Brownian Dynamics simulation of reversible polymeric networks

Non-l~zmmi~ Fluid Medmuics ELSEVIER

J. Non-Newtonian Fluid Mech., 60 (1995) 303 334

Brownian Dynamics simulation of reversible polymeric networks B.H.A.A. van den Brule *, P.J. Hoogerbrugge Shell Research B.V., P.O. Box 60, 2280AB-Rijswi/k, The Netherland*

Received 30 May 1995

Abstract In this paper we describe the construction of a Brownian Dynamics simulation of reversibly cross-linked networks. The simulation differs from existing analytical approaches and computer simulations of networks in the sense that we also take the topology of the network into account. The motion of the junction points between different molecules is not prescribed but is calculated from a force balance. This makes it possible to measure the effect of network reorganisations on the stress relaxation. The response of networks to shear flow is measured and analysed in terms of transient network theory and within the framework of linear viscoelasticity. It is shown that the average motion of the junction is affine but that there is a long time diffusive process around the affine path. It was found that, even in systems with Gaussian chains and fixed association and disociation rates, a shear thickening of the viscosity and primary normal stress coefficient can occur. The reason was found to be that dangling segments are recaptured by the network before they had the opportunity to fully relax to the equilibrium state where the probability of reattachment to the network increases linearly with the length of a segment. Due to this mechanism the fraction of long segments present in the network is increased. This explanation Of shear thickening seems to be consistent with experimental findings. Keywords: Brownian Dynamics; Network theory; Non-affine motion; Polymer gels; Rheol-

ogy; Shear thickening

* Corresponding author.

SSD1 0377-0257(95)01378-4

304

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Flukl Mech. 60 (1995) 303 334

1. Introduction

In recent years there has been a growing interest in reversibly cross-linked polymer gels, partly because of their intriguing theological properties but also because they find wide application in the food, paint and oil and gas industries. In industry, associative polymers are often used as thickeners to increase the viscosity of a liquid, e.g. to reduce the settling velocity of suspended particles, whereas in other applications it might be important to create a fluid with a shear rate dependent viscosity. There are many ways of generating reversible networks. Often, the polymer chains contain specific groups which may e.g. interact by hydrogen bonding or ionic association, e.g. borate cross-linked hydroxypropyl-guar a system widely used in the oil industry. Also, there are systems where the molecules may interact due to the formation of double or triple helical junction zones (e.g. Xanthan and Shellflo-S) or by partial crystallisation. Other interesting systems are based upon a solution of a multi-block polymer in a solvent with a very different quality for the different blocks; under suitable conditions of temperature and concentration the virtually insoluble parts of the polymer will aggregate into micelles which are connected by the soluble parts. From a rheological point of view reversible gels differ markedly form "ordinary" polymer solutions without cross-links. Quite often it has been observed that there is a significant amount of shear thickening, a failure of the Cox-Merz rule, and a high value of the dynamic moduli which extends to very low frequencies. Due to the presence of the cross-links between the macromolecules it is attractive to analyse the behaviour of reversible gels within the context of network theory. The first network theories for polymeric materials were developed some fifty years ago to model the behaviour of permanently cross-linked polymeric materials such as vulcanised natural rubber. The essence of network theories is that the polymer polymer interactions are assumed to be localised in so-called junctions. In rubber, a junction is a point at which two macromolecules are chemically cross-linked. It is commonly assumed that the junction points between the chains move as if they are embedded in the macroscopic continuum (affine motion). The macromolecular material connecting two junctions is called a network segment. The network is described in terms of the lengths and orientations of the individual network segments and in doing so the actual spatial topology, i.e. the way the segments are connected to each other, is not taken into account. Using an entropy argument it can be shown that a network segment can be represented by an elastic spring which for small extensions can be taken to be linear. A comprehensive review of the network theory for rubber can be found in the excellent monograph by Treloar [1]. Later, the theory for rubber elasticity was extended by Green and Tobolski [2] to allow for the stress relaxation observed in reversible networks. In this transient network theory the junctions are not of permanent nature like a chemical cross-link but of temporary nature, so junction points be created and destroyed. Yamamoto [3] and Lodge [4] further explored the concept of transient network models in order to find a description of the behaviour of entangled polymer melts. A review of

B.H,A.A. lan den Brule, P.J. Hoogerbrugge

J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

305

network theories and their associated constitutive equations for polymeric liquids can be found in the books by Bird et al. [5] and by Larson [6]. Most transient network theories have in c o m m o n that they assume the rate of creation of junctions to be constant. The rationale is that the creation of junctions is driven by diffusion. The theories can be discriminated in the way they account ti~r destruction of junctions. In models that lead to closed form constitutive equations for the stress tensor the probability of destruction of a junction is either constant or it depends on global quantities like the average extension of segments. An example of a theory with a constant rate of destruction of junctions is the Lodge rubber-like liquid model [4] and an example with length dependent destruction is the model of Phan-Thien and Tanner [7]. From a physical point of view it is more natural to assume that the probability of a segment detaching from the network depends on its own length instead of being dependent on some average quantity as lbr instance the average length of all the segments. To investigate this problem Petruccione and Biller [8] performed stochastic simulations of a network where they followed the evolution of an ensemble of segments under a prescribed deformation. It should be noted however that in these simulations the topology of the network is still not taken into account. For this reason, creation or destruction of a segment does not lead to a further reorganisation of the network and a subsequent relaxation of the stress. The drawback of network models for polymer melts is that there is no physical modal available to account for the formation and destruction of network segments. Due to this problem the network models for melts have to be classified as phenomenological rather than molecular. In modern theories for polymer melts the topological interactions between molecules are taken into account by means of the reptation concept; a molecule can effectively only relax its configuration by a diffusive motion along its backbone since lateral motions are hindered by surrounding molecules. Recently~ a series of papers appeared which focussed on the modelling of reversibly cross-linked polymer gels. These systems can be divided into two classes. First, we have systems where the molecular weight between cross-links is large compared to the molecular weight between entanglements. These entanglements may persist on a time scale long compared to the average lifetime of a junction. For this reason the dynamics of such a system is dominated by the presence of these long living entanglements. The molecular motion in an entangled network can be studied using the concept of "sticky reptation". A recent contribution in this area is the work by Leibler et al. [9]. Also, there are systems with a molecular weight between cross-links that is small compared to the entanglement molecular weight. In these systems the number of cross-links exceeds the number of entangelements and the dynamics is governed by the formation and destruction of the cross-links or junctions. In this paper we want to focus on these materials. The dynamic behaviour of unentangled, reversibly cross-linked networks has been recently investigated by T a n a k a and Edwards [10]. In their model, the polymer chains are all of equal molecular weight and have sticky functional groups attached at both ends of every chain. These sticky end groups associate the chains to junctions where it is assumed that there are no limitations on the functionality of

306

B.H,A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

a junction. It is assumed that the chains obey Gaussian statistics and hence that the force transmitted through a chain increase in proportion to its length. Two types of chains are discerned: active chains in which both ends are connected to separate junctions and dangling chains which are only connected to the network at one of their ends. When an active segment becomes dangling it is assumed that it relaxes into the equilibrium state on a time scale which is of the order of the Rouse relaxation time of a single chain. This time scale is very short compared to the average lifetime of a junction. Under equilibrium conditions the detachment rate of chains from a junction is governed by thermal fluctuations but under flowing conditions the chain breakage probability depends on the extension of the chain segment under consideration. With regard to the motion of the junctions the classical assumption of network theory, that of affine motion, is adopted. To explain the shear thickening behaviour which is often observed in polymer gels a number of modifications of the T a n a k a - E d w a r d s theory have been put forward. Wang [11] explained shear thickening by a shear induced coagulation of free chains to the network thus increasing the number of active chains. Marrucci et al. [12] proposed two modifications of the T a n a k a - E d w a r d s model. First, they took the finite extensibility into account and showed that only a tiny amount of shear thickening results from this model. Next, they introduced the concept of partial relaxation of dangling segments into the theory. When a chain dissociates from a junction the just released sticky end starts to move through the network in order to relax the elastic stress. It will however not succeed in relaxing completely since the spacing between the junctions is typically small compared to the length of a stretched chain and hence the dangling end is soon recaptured by the network. It was shown that this model leads to a significant shear thickening and possibly to the correct scaling of the flow curves. It should be mentioned that there are also other explanations for the shear thickening behaviour of cross-linked gels. Some years ago Witten and Cohen [13] proposed the idea that with increasing shear rate the probability of formation of an intermolecular association increase over that of forming an intramolecular one, thus effectively increasing the molecular weight and the viscosity. Recently, Ahn and Osaki [14] studied shear thickening behaviour using a network model with strain dependent creation and loss rates. Here we report on a simulation of a reversible network similar to the one proposed by Tanaka and Edwards [10], where we take the topology of the network into account. We can either prescribe the motion of the junctions but we also have the possibility of calculating the actual junction motion from a force balance. Furthermore, processes like association and dissociation and the resulting relaxation of the entire network can be considered in detail.

2. Description of the simulation The networks we consider in this paper consist of a large number of macromolecules with uniform molecular weight and it is assumed that the molecular weight between the junctions is small compared to the entanglement molecular

B.H.A.A. lan den Brule, P.J. Hoogerbrugge / J, Non-Newtonian Fluid Mech. 60 (1995) 303 334

307

weight. For this reason all the interactions between molecules are concentrated in the junctions and there are no topological constraints on the motion of the chains. To further reduce the amount of computational work, a polymer chain is modelled as an elastic dumbbell. The beads of the dumbbells are sticky which causes the gelation of the system. For the physical gels we want to study here~ the bond energy between two beads sticking together is typically in the same range as the thermal energy so that bonds may be created and destroyed. In this way a temporary network of elastic dumbbells, modelling a reversible gel, is created. In order to simulate an infinite network in shear flow, periodic boxes are used that are allowed to slide over each other. It turned out to be non-trivial to find a way to correctly implement the association/dissociation process in the simulation. The natural way to implement association is by making the beads attractive; neighbouring beads will then attract each other and form a junction. Now the problem arises that the attractive force from a junction consisting of two beads on a third bead will be doubled etc. The result is that a few clusters will rapidly grow in size and eventually the system will collapse. The natural solution to this problem is to use beads of a size larger than, or comparable to, the range of the attractive force. This solution comes at the expense of very small time steps in the simulation and is thus not considered practical. Another option is to introduce a critical radius d in the model; two beads separated by a distance less than a certain critical distance are put together in a junction. This would indeed yield a simple way to associate the beads but now the dissociation process beomes problematic: after dissociation of a junction there will be an immediate reassociation since all the beads just released will be within the critical radius. A way to solve these problems is by the introduction of a virtual particle we called a "sticky node". "File beads associate to these sticky nodes rather than to each other. The way the model works is as follows. Say we start with a free node. A free node jumps through the simulation space in a random fashion. It continues doing this until it jumps to a position within the critical distance e/from a bead. When this occurs it will j u m p to the position of the bead and attach itself to the bead. Next, the node will move with the bead until they find another bead within the critical distance. When this happens the node will glue the two beads together to form a junction at the position of their centre-of-mass. This process continues until the number of beads glued together by a node equals the maximum functionality Ji,~,~ of the node. Every time step cSt the node has a certain probability P to detach itself from the surrounding beads and it will j u m p to a random new position in the simulation space after which the whole process starts again. In this paper the probability of detachment is taken to be constant, P = 1 - c x p ( - ~ t r) where r is the average lifetime of a node sticking to a bead. A detailed discussion of the resulting composition of the network in terms of the model parameters can be found in the Appendices. The actual simulation of the motion of the beads is based upon the method of Brownian Dynamics. The dumbbells are suspended in a continuous solvent with a viscosity /L and a homogeneous velocity field given by v~(r)= vo + K ' r , where K denotes the transpose of the velocity gradient, ~,-,/= ~t~,,c~.vi. It is assumed that the

308

B.H.A.A. van den Brule, P.J. Hoogerbrugge /J. Non-Newtonian FluM Mech. 60 (1995) 303 334

beads do not disturb the velocity field of solvent which means that hydrodynamic interaction effects are not accounted for in this model. The beads experience a viscous drag force when they move with respect to the solvent. This viscous drag force is given by F "'~= - ~ ( v - v s ) where ( is the friction coefficient. Besides the (systematic) viscous drag force there will be a stochastic component in the interaction between the bead and the solvent. Regarding the stochastic force F ~ we adopt the usual assumption that the stochastic force is uncorrelated on the time scale of the motion of the particle and define [15] ( U s~) = 0

(2. i )

(F(S)(t)F~S)(t + r ) ) = 2(kTd(r)l,

(2.2)

and

where k is the Boltzmann constant and T is the temperature. The angular brackets indicate an expectation value which for our numerical simulation equals the mean taken over the beads in the system. The elastic connector between two beads gives rise to a force F (c) = HQ where H is the spring constant which might be a function of the length of the connector and where Q is the vector connecting the two beads of the dumbbell. In this paper, however, we restrict ourselves to systems with a linear spring force-law. In order to keep the system homogeneous, i.e. to avoid the gradual build-up of a concentration gradient in bead positions, it turned out to be necessary to introduce a repulsive force between the beads. To limit the number of interactions of a bead with surrounding beads an interaction range a is introduced and the repulsion force exerted by a neighbouring bead at a relative position s is given by Fro(s)= _ F 0 ( 1 s ~ s _ a/ \

for

s
F(r'(s)

=

0

for s > a .

(2.3)

s

The interaction range is t y p ~ e l e c t e d in the order of 0.5 times the equilibrium length of the dumbbells (x/3kT/H) and for the densities used in our simulations this means that the number of neighbours in this range is between 10 and 30. In the simulations the strength of the repulsive force is tuned in such a way that the contribution of the repulsive force to the pressure is approximately nkT, thus balancing the contribution of the connector forces to the pressure which is known to be -nkT. It turned out that indeed the network stayed homogeneous in this way (probably it is possible to obtain homogeneous networks with smaller repulsive forces since the Brownian motion also assists in keeping the network open. This point needs further investigation). The contribution of the repulsive forces to the stress tensor remained isotropic under all circumstances investigated (the contribution of the repulsive forces to the deviatoric part of the stress tensor was always less than 1%). The repulsive force also had no effect on the association/dissociation statistics. Neglecting inertial effects, the forces mentioned above should balance and the following Langevin equation for the velocity of a free bead, i.e. a bead not connected to a node, can be formulated:

B.H.A.A. van den Brule, P.J. Hoogerbrugge , J. Non-Newtonian Fluid Mech. 60 (1995) 303 dd4

v = co + r • r + ( U ~1 + r F ~) + F ~ ) / ( .

309

(2.4)

To obtain the velocity of a junction, i.e. a node plus the f beads connected to it, a similar equation can be used when we simply sum up all the forces working on the beads attached to the node and replace the friction coefficient ~ byJ~L It can be shown that with this procedure the internal energy of the system does not change due to association or dissociation of the the polymers to the network. A straightforward method to obtain a numerical algorithm to simulate a Langevin equation is the one proposed by E r m a k and M c C a m m o n [16]. In this approach the force balance is integrated over a numerical timestep ~t where the time increment is small enough to assume that the configuration of the system does not change appreciably during a timestep. For this reason the systematic forces can be considered to be constant during the integration. The random force, however, goes through many fluctuations. In this way the following algorithm which generates a new configuration from an existing configuration is obtained: r(t-+- ~t) = r ( t ) + vomit + K • r(t),~t + 1 [U,,~(t) + XF'"'(t)]bt + ~h',

(2.5)

where the stochastic step dr is given by

(~r ==-1 ~ U~>(t ') dt'.

(2.6)

J t

Using the properties of the stochastic force it can be shown that the random step can be sampled from a Gaussian distribution with zero mean and variance c;2= ( 2 k T (),~tl. It can be shown that without loss of accuracy the random step can also be sampled from a uniform distribution on the interval [ - , 3 o -2, , / 3 a 2] [17]. In the simulations the time is made dimensionless with the relaxation time 2 = (/(4H) of a free dumbbell and the lengths are scaled with a measure of the equilibrium size of a dumbbell , , / k T / H . In order to compare the actual motion of a network to the prescribed deformations which are commonly assumed in the analytical network models (for an overview see, e.g. Larson [6]) we also build in the possibility of defining the motion of the junctions. The dumbbells with both beads connected to a junction now act as a network segment with a prescribed motion as in standard network theory and the dangling ends, with only one end connected to the network, will relax in a Rouse-like manner. This model can be analysed analytically and will be considered in more detail in Section 4.

3. Simulation results in shear flow It is assumed that the extra stress tensor, i.e. the deviation of the stress tensor from the equilibrium value, can be written as a sum of a Newtonian solvent contribution and a polymer contribution, r = r~ + %. In the remainder of this paper we will only consider the polymer contribution. The Kramers expression for the polymer contribution to the stress tensor reads [5]

310

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303-334

rp = - n (F~C)Q)

(3.1)

where n is the number density of dumbbells and where we ignore isotropic contributions. In the simulation only a limited number of chains is considered and the averaging is performed over the dumbbells present in the simulation domain. In shear flow, the velocity field is given by vx = ~y, v.v = 0, and v= = 0, where ;~ is the shear rate and may be a function of time. In start-up of shear flow the fluid is initially at equilibrium. For t > 0 a constant shear rate ~ is generated, and the stresses grow until they eventually reach their steady-state values. The material functions of interest in this flow field are the viscosity 4, the primary normal-stress coefficient q~ and the secondary normal-stress coefficent q~2. The material functions are time dependent and are defined by tAx(t) = - v/+(~,t)~,, r,.~(t)

r.~.x(t)=

-

-

q ~l+t T", t ") 72,

(3.2)

r y x ( t ) - r=(t) = - qJ+ (~),t)~2. Since the stochastic processes in the y-direction and the z-direction are identical and since we use a linear force law, the secondary normal-stress coefficient will be equal to zero for the model considered here. The steady-state values of the viscosity and the primary normal stress coefficient are obtained from the limiting behaviour of the transient material functions at long times: lim ~/+(~,t),

~(~) = - r , . , . / f , =

t~

(3.3)

z¢

qJl(7) = - ( r x , - r ~ , ) / f '2= lim q~/-(~,t).

(3.4)

t~JG

For small or slow deformations, the response of a material can be described within the context of the theory of linear viscoelasticity. The shear stress of a linear viscoelastic material is related to the history of deformation by t

v,.,.(t) = .f G ( t -

t')~(t') dt',

(3.5)

where the function G ( t ) is the relaxation modulus. One way to measure the relaxation modulus is in a step strain experiment. At t = 0 the material is suddenly sheared by an amount 7. The shear stress measured in such an experiment is given by r,..,.(t)/7 = G(t).

(3.6)

Often, the relaxation modulus is expressed in terms of a discrete spectrum {At, Gi}, G ( t ) = y , G e -'/;''. i

The transient viscosity is given by

(3.7)

B.H.A.A. van den Brule, P.J. Hoogerbrugge

J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

q+(t) = i G(s) ds = ~ Gi2/(1 - e -t.~,),

311

(3.8)

0

and the steady-state viscosity is thus found to be

tl = i G(s) ds = ~ Gs2t~.

(3.9)

o

For the primary normal stress coefficient no expression can be given within the theory of linear viscoelasticity since the first normal-stress difference grows quadratically with shear rate. Using the principle of frame invariance, Larson [6] showed that the lowest order term of the expansion for the first normal-stress difference is related to the relaxation modulus by l

• ~(t)=2

sG(s) ds=2

G~Jt 1 - = e

- e - ' ~' .

(3.10)

0

The steady-state value of the first normal-stress coefficient is thus found to be equal to

u?,.,,=2

f

sG(s)ds=2~G,22,

(3.11)

i

0

which is a well known result commonly derived from simple fluid theory. The storage modulus and the loss modulus, G'(~o) and G"((~), are related to the relaxation modulus by the Fourier transforms:

G'(e)) = so

G(t) sin((ot)dt

~-

GiAiO-

(3.12)

2 ,' 1 +2t(o-

(3.13)

• 1+

zTo)-

0

G"(~o)=a)

G(t)cos(e)t)dt= 0

3. 1. Simulated material functions A series of simulations was carried out to assess the influence of several parameters in our model on the stress response. To this end, a periodic domain of (20 x 5 × 5) x / - ~ / H was filled with 2000 dumbbells with random position and orientation and equilibrated. Then, a number of nodes, between 1000 and 2000, was inserted at random positions and the system was again allowed to equilibrate for 150 to 1500 2, depending on the dissociation time constant r. Starting from such a well-equilibrated state, a shearing deformation was suddenly started using shear

312

B.H.A.A. van den Brule, P.J. Hoogerbrugge /J. Non-Newtonian Fluid Mech. 60 (1995) 303-334

rates o f 2;7'-----0.01, 0.03, 0.1, 0.3, 1, 3 and continued for 150 to 1500 ~,. The shear stress and the normal stresses were recorded as described in the previous section. The typical time step size used in the simulations was 0.01 2. In Table 1 an overview is presented o f the parameter settings o f the simulated systems. F o r a system o f free H o o k e a n dumbbells, an excellent agreement was f o u n d between the simulated stress response and the theoretical prediction. Variation o f the time step had no effect on the results for free dumbbells. In an associating system, the free nodes move to a r a n d o m new position every time step. A change o f the size o f the time step thus affects the mobility o f the free nodes. In order to keep the statistics o f the system invariant under a change o f the time step it is therefore necessary to keep the frequency o f association attempts constant under a change o f the time step. This was indeed confirmed with a simulation where we reduced the time step by a factor ten, and an association attempt was made only every tenth step. F o r all systems were no shear was applied, x x , y y and z z c o m p o n e n t s o f the polymer contribution to the stress tensor were equal to n k T within statistical error, independent o f the association/dissociation processes taking place. Fig. 1 shows the transient dimensionless viscosity vs. time for different shear rates for system 4. In this particular case the viscosity is independent o f shear rate, because the rate of association was low e n o u g h to allow dangling segments to fully relax to the equilibrium state. In Fig. 1 it can also be seen that the simulation results for the viscosity o f the non-affine network are much lower than for the affinely moving network. In the non-affine network, network reorganisations can take place a r o u n d dissociating and associating segments which results in a significant relaxation o f the stress. In Fig. 2 the same behaviour as described above for the viscosity is observed for the first normal stress coefficient. The y y and z z c o m p o n e n t s o f the stress tensor are not affected by the shear, and are therefore always equal to n k T within statistical error. Assuming that our simulated systems are dominated by a single, slow, relaxation mechanism with a measured relaxation time 2m and a modulus G m w e can derive the following relations: Table l Systems l-10: parameter settings and resulting fraction of active segments at low shear rates No.

~

d

l 2 3 4 5 6 8 9 10

l0 l0 10 30 100 10 10 t0 10

0.1 0.175 0.3 0.175 0.1 0.175 0.175 0.175 0.175

Jn,~

Nnoo~ 4 4 4 4 4 6 10 4 4

4D<,theory

q$<,affine

0.542 0.672 0.750 0.794 0.826 0.812 0.856 0.873 0.929

0.500 0.671 0.790 0.791 0.801 0.816 0.865 0.880 0.937

1000 1000 1000 1000 1000 1000 1000 1500 2000

System 7 never existed. Error in ~b, is 2 x 10 3.

~b,, non-affine 0.493 0.688 0.810 0.803 0.802 0.816 0.843 0.889 0.941

B.H.A.A. ran den Brule, P.J. Hoogerbrugge ,' J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

25

313

i i

o

20 i

c"

.. I

o

I

m

/t

a,ine

/ -

non-affine

-5'

o

2'0

4'0

60

80

lOO

time

Fig. 1. Growth of the viscosity upon inception of shear flow tbr atline and non-affine motion. Shear rates are ,~;= 0.1 and 0.3. The full line represents the analytical result. Simulation system 4.

Jl~ = 2,,,

G,,1

and

%t'~'° - 2m.

(3.14)

2r/o

In Fig. 3, tl",.0/2r/o is plotted aginst the measured relaxation time. It can be seen that indeed the response seems to be dominated by a single relaxation time which becomes shorter when we move from an affinely moving system to a non-affine system. Since the system contains a total o f n segments in a unit volume, we expect 1200, i

0

affine

lOOO~ i

.J--

800

j .-//

*~

600 ~

o

400

=

2o0 ~

=

0

c

" "/

-~-~--~

~~-"i2~-

I -200

' 0

2'0

4'0

60

80

1 O0

time

Fig. 2. Growth of the first normal-stress coefficient upon inception of shear flow for affine and non-affine motion. Shear rates are 2;:,=0.1 and 0.3. The full line represents the analytical result. Simulation system 4.

314

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303-334

l

z a~

i

_ e-

~0

./

t

c 0 c)

E

I0~

J i

t

nl

i /J 10 Psi1 / 2 visc Fig. 3. Relation between the viscosity and the first normal-stress co-efficient for affine and non-affine systems. The dashed line represents the relation 9°1,o= 2qo2~. al refers to an affine simulation of system 1 and nl to a non-affine simulation etc. that, in a n a l o g y to elastic d u m b b e l l t h e o r y o r n e t w o r k theory, G(0) = n k T . In Fig. 4, we p l o t t e d % / n k T vs. the r e l a x a t i o n time. N o w we observe that the expected b e h a v i o u r is o n l y f o u n d for the affine systems. In the non-affine systems the viscosity is r e d u c e d b y a larger f a c t o r t h a n the r e l a x a t i o n time. This suggests t h a t in these systems the effective m o d u l u s is lower t h a n n k T . T h e r e l a x a t i o n m o d u l u s G ( t ) is o b t a i n e d f r o m a step strain simulation. Despite the fact that, in o r d e r to o b t a i n a c c u r a t e results, we h a d to r e p e a t the e x p e r i m e n t s a n u m b e r o f times, this p r o c e d u r e t u r n e d out to be m o r e efficient t h a n use o f a a5 J

n5

~ /~I 0 J

/ J

E n4

6 E

/-

I0~

n1~9

/./

./

//

!nl

/

/

/~

J

.J

10 normalised viscosity Fig. 4. Relation between the viscosity and the measured time constant for affine and non-affine systems. The dashed line represents the relation % / n k T = 2 m.

B.H.A.A. van den Brule, P.J. Hoogerbrugge ,, J. Non-Newtonian Fluid Mech. 60 (1995)303 334

315

affine: simulation affine: two-mode theory - non-affine: simulation non-affine: multi-mode - - \

"\ "\

"\

0.1~ i

o

\

"-~.

,i

"'X

'

\.

0.01

0

20

4'0

6'0

I

I 100

80

time Fig. 5. Relaxation modulus G(t)/nkT of a typical system for both affine and non-afline motion.

G r e e n - K u b o relation. In Fig. 5, the resulting relaxation modulus for a typical system is shown, both for the affine and the non-affine case. It can be seen that, in agreement with expectations, G(0) = nkT. From the final slopes o f the two curves it is clear that the dominant relaxation time in the non-affinely moving system is slightly shorter• The value of the modulus associated with this slow process is significantly lower than nkT as can be seen by an extrapolation of the final part of the curve to t = 0. The major difference between the two situations appears to be due to the relatively fast initial relaxation in the non-affine system.

5 0

E

#,

0.1q i

0 c-(13

G' G"

O~

o

0.01

0,01

0.1 Omega

1

Fig. 6. Dynamic moduli for system 4 in non-affine motion.

--

316

B.H.A.A. van den Brule, P.J. Hoogerbrugge /J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

E~J

viscosity normal stress

(.J

~E

+k_

O) 0

100

E 0 c

c

10~

f i .

.

0.01

.

.

.

.

.

,

.

.

0.1

.

.

.

.

.

.

.

.

1

.

.

.

.

.

10

shear rate

Fig. 7. Shear thickening behaviour of the viscosity and the first normal-stress coefficient at high shear stress. The dynamic moduli are obtained from the measured relaxation functions. In Fig. 6 the moduli for a non-aNne simulation of system 4 are shown. In a number of systems a significant shear thickening was observed, in both the viscosity and the first normal-stress coefficient. This in spite of the fact that we used linear springs and a fixed probability for segments to dissociate from the network. An example is shown in Fig. 7. Shear thickening was found to occur in those systems where a dangling segment has a high probability of reassociating to the network again before it has fully relaxed to the equilibrium state. In agreement with experimental findings the onset of shear thickening occurs at a critical value of the shear stress rather than the shear rate (see Fig. 16 below). 3.2. Actual motion o f the network compared to affine motion

In a non-affinely moving system, the motion of an active segment depends not only on its own length and orientation but it is also affected by the motion of the rest of the network which is subjected to continuous rearrangements due to creation and destruction. To eventually include non-affine motion in an analytical model, it is necessary to find a way to take the influence of the topology of the network into account in some average sense. In order to gain more insight into the dependence of the non-affine motion on global properties several additional simulations of networks with non-affinely moving junction points have been performed. First, the relation between the motion of a junction and its age is studied. A junction is " b o r n " at t = to, when a node associates to its first bead and "dies" when the node disconnects from the surrounding beads. The age of a junction is represented by 0 = t - t o . Since in the simulation the time progresses in discrete steps we introduce Oi = iSt where i = 0 at t = to. The difference s between the affine

B . H . A . A . z an den Brule, P.J. Hoogerbrugge ' J. N o n - N e w t o n i a n Fluid :14ech. 60 (1995) 303

334

317

path and the actual path of a junction is calculated at each time step, and these differences are accumulated over the lifetime of the junction:

s(O,,,)= ~ [~i-v(r,)]bt, i

(3.15)

O

where r~ is the position of the junction at time O, = to + i,~t. Changes in position that take place at the moment that additional beads associate to the junction are not taken into account. Over the lifetime of a junction this happens only a few times, and therefore the influence will be very small. In a homogeneous flow field, due to symmetry, the relative motion of a junction with respect to the macroscopic continuum is just as likely to be in one direction as in the direction exactly opposite to this. For this reason, the relative displacement averaged over all junctions will be zero, irrespective of age:
(3.16)

So, the average motion of a junction is still affine. This, however, does not imply that the network segments deform in an affine way as well. Junctions at opposite ends of an active segment have a general tendency to move towards each other, thereby allowing the segment to relax its stress. (This is similar to the motion of an elastic dumbbell; the average motion of the beads is affine but the motion of the connectors between the beads is certainly non-affine.) To obtain non-trivial information regarding the junction motion, it is necessary to record the mean-squared relative displacement of junctions as a function of junction age, or rather the dyadic product . A typical case is shown in Fig. 8. At low ages, the junction motion is dominated by Brownian forces on a single bead, and the diffusion is isotropic. At high ages, the reduction in the slope is much 9 i

8~

E E o

/-

C m

E

i

j/

7i!

/

61

.

YY ZZ XY

.ffl "o

0"

.f"

a4

/./

/

//

/

/

i

27

ii/, / 0/ 0

//

2'0

4'0 60 age of junction

8'0

1 O0

Fig. 8. Mean-squared relative displacement of junctions as function of their age. Simulation system 4, shear rate 0.1.

318

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

Table 2 Comparison between stress tensor and long-time diffusion tensor components for various systems and shear rates System

4 4 4 4 8 10

)~

0 0.01 0.1 1.0 0.1 0.1

Stress

Diffusion

Diff/stress

xx

yy,zz

xy

xx

yy,zz

xy

xx

yy,zz

xy

1 1.03 3.65 363 2.18 6.69

I 1 1 1 1 1

0 0.094 0.92 10.8 0.69 1.40

0.017 0.017 0.075 6.0 0.03 0.3

0.017 0.017 0.017 0.017 0.013 0.03

0 0.0008 0.008 0.08 0.0042 0.01

0.017 0.017 0.021 0.017 0.014 0.045

0.017 0.017 0.017 0.017 0.013 0.03

--0.0085 0.0087 0.0074 0.0061 0.007

more than can be expected from the fact that the junction will have a higher functionality. Nevertheless, the slope becomes constant again for all components, which clearly indicates that the junciton motion is again diffusive. At these high ages, however, the diffusive motion is dominated by elastic forces, which fluctuate due to the continuous reorganisations taking place in the network. Because the magnitude of the elastic forces on the junction depend linearly on the size of associated segments, it is expected that the long-time diffusion tensor, determined from the slopes at high ages in Fig. 8 is proportional to (QQ) and therefore to the stress tensor. As can be seen in Table 2, the diagonal elements of the long-time diffusion tensor and the stress tensor have a very similar proportionality constant. The xy-component of the tensor shows a different proportionality, about a factor 2 smaller, the cause of which is currently not well understood. The motion of segments has also been studied in a more direct way. For this, the average dependence of Q and Q on Qx has been recorded for active, dangling and free segments. Also, the number density and the association/dissociation probabilities have been recorded as a function of Qx. As expected, there is no dependence of , on Qx at equilibrium. A typical result for
B.H.A.A. ran den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

319

r

-0.2

'~ -

-~

j'

A

.

acti,ve

-0.4 o

"o v

-0.6 -

-0.8 -

. j,

dangling

free -11 0

0.5

i

115

2

2.5 Qx

:3

3'5

4

4.5

5

Fig. 9. Relation between Q, and (Q,) for active, dangling and free segments. Simulation system 4. system at rest. the segment size Qx- However, especially for the active segments, the curve has become very non-linear, in contrast to the equilibrium case. Small active segments now even show a tendency to stretch rather than to relax. For large segments, the tendency to relax becomes less and less Q,.-dependent. This is probably due to the fact that large active segments are likely to have been active a long time, and much of their possible relaxation has already taken place. The only way such segments can relax further is by dissociation processes of other segments in the surrounding network. This will yield relaxations in Q proportional to the size of the released

cO

Q1

\\\

~

t~

active

' \ \

~langlin&

"cl

O.Ol c

-,\

E

f~ee,

(1)

,,

o.ool 7

o.oool

, 0

05

1

is

2

25

3

35

4

4~

Qx

Fig. 10. Distribution of segment lengths in the x-direction. System at rest.

320

B.H.A.A. van den Brule, P.J. Hoogerbrugge /J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

2.5!

2~

active ' ~,,.v

I

/~"

'~'~

15 ! A

]

'

.

0v 0.5"

,, ~//'/free

o

2

,;,

~

8

1'o

1'2

4

Qx Fig. l l. Relation between (Q,.) and Q~ for active, dangling and free segments. Simulation system 4, shear rate 0.1.

other segments, rather than its own size. Such a size-age correlation does not exist at equilibrium. For the region 0.6 < Q,, < 2, the relaxation of dangling segments can still be described accurately as half the active plus half the free relaxation.

-0.2 4I A

::., ,,/,~¢j,

i!

-0,4 4

II

i

O "O V

qtanglipg i, !

If i " '?

!

i

i', I /

(

-0.6 1

free

~

-0.8 !

iI -1! 0

2

~

6

Qx

S

1'o

1'2

14

Fig. 12. Relation between (Qv) and Q~ for active, dangling and free segments. Simulation system 4, shear rate 0.1.

B.H.A.A. ~an den Brule, P.J. Hoogerbrugge

J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

321

[ o~ -0.2

active

,

/

I \

A

-0.4 ~

/

,~\

o V

-0.6 -

i

~"~ "~

~l

i

" - dangling

osJ

,

free -1 i 0

, 2

, ", 4

4

, 6

8

1'0

1'2

14

Qx Fig. 13. Relation between
4. Analysis of the results

4.1. Analytical model In the analytical model two types of segments are distinguished: active segments which are connected to a node (junction) at both ends and dangling segments which are only connected to one node. The number of free segments is usually small and their behaviour is similar to the dangling segments. Variables which refer to the active segments are indicated by a subscript ' a ' and variables related to the dangling segments are indicated by a 'd'. The average life time of a node is given by the time constant r. The probability that an active segment becomes dangling between t and t + dt can therefore be written as h,.,l dt where h,,j = 2/T. Likewise, the probability for a dangling segment to become active is given by h,~., dt. To specify the orientation of dangling and active segments the following distribution functions are introduced with i = a,d:

tFi(Q,t)d3 Q = number of segments of type i with a configuration in the range d3Q about Q at time t in a unit volume.

(4.1)

A balance in configuration space yields the following equation for the evolution of the distribution function:

~t u?i = ___Q . (Qudi) + hj. tFj_ h u ~ i

(ij = a,d and i C j).

(4.2)

The first term on the right-hand side of these equations accounts for the convection of segments while the last two terms account for the creation and destruction of segments. Integration of this equation over configuration space yields the following relation:

322

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

hd,~na = h,,an~,

(4.3)

where na and na denote the number density of dangling and active segments, respectively. In order to obtain a constitutive equation, the motion of the segments needs to be further specified. Based upon the expression for the average velocity of a free elastic dumbbell with orientation Q the following generic expression is postulated:

kT Q = ~c. Q - ~t-~- ~

in ud~ - ½fl~Q

(i = a,d)

(4.4)

which reduces to the usual equation for an elastic dumbbell if we choose c~ = 2 and fl~ = 4H/(. If we select ~i = fli = 0 we will have affine convection of the segments, as in standard network theory, and if we select c~ = 1 and fit = 2 H / ( the segment will behave as a dumbbell with one bead moving affinely and one bead subjected to Brownian motion. Going through the usual procedure of standard network or dumbbell theory the following equation for the second moment of the distribution function can be derived: d -~t ( Q Q ) t -

K . (QQ),-

( Q Q ) , . K*

2o~intk T - ~ 1 - fl~(QQ)t + h j . ~ ( Q Q ) j - h ~ j ( Q Q ) ,

( i j = a,d and i C j). (4.5)

In this equation the notation ( . . . ) ~ refers to the weighted average over the distribution u?~. Using the Kramers expression for the polymer contribution to the stress tensor, and ignoring isotropic contributions, we have Tp,i =

-H(QQ)~

(id= a,d).

(4.6)

Elimination of ( Q Q ) i between Eqs. (4.5) and (4.6) yields the following equation for the stress tensor:

flt)z 1 iJ 1 + h~./] p.t + --h~j--6t Vp.t =

o~,n,kT ni 22ht¢ 1 + ny--zpj

( i j = a,d and i -¢j),

(4.7)

where 6/c~t denotes the upper-convected time derivative,

6/6t(..) = D/Dt(..)-

Vv*. ( . . ) -

( . . ) - Vv,

(4.8)

and 2 = ( / 4 H is the relaxation time of a free dumbbell.

4.1.1. Affinely moving network First, we will analyse a model with affinely moving active segments (~, = flu = 0) and dangling segments with one end moving affinely and the other end subjected to Brownian motion (c~d = 1 and fld = 2 H / ( = 1/(22)). Using ha.d= 2/r and hd.a = 2na/ (rnd) we obtain for the active segments z ~ na Vp., + -~ -~t r°'~ = - - Vp.a, nd

(4.9)

B.H.A.A. ran den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 -334

323

and for the dangling segments 42n,~

5

42

1 + rna /]'rP"l -}- 22 ~ '~'p.d= - n a k T l + - - %,,. T

(4.10)

The viscosity and the primary normal-stress coefficient for this model can be found to be

(4.11)

tl = nkT2c~,A(1 + e) and k~JI =

2nkT.;t2gb,,A:(1 + 2~b,, + s2q~,,),

(4.12)

where we introduced ~b~ = %/n, ~,, = na/n, A = z/(2Z) and ~:= 2/(Aq~,q~a). The linear viscoelastic properties can be written as a spectrum {).,,G~} consisting of two modes given by 21.2 = c_+ x f ~ - 2A,

< : - nkT[-2 [_ / 1 + X/ 1

(4.13) s(c272A)2A -1 ,

(414/

where c = A/2 + 1/~ba. In Fig. 6 the calculated relaxation modulus is compared with the simulation results. The agreement is excellent. Also, the calculated response to inception to shear flow compares well with the affine simulations as can be seen in Figs. 2 and 3.

4.1.2. A non-affine model The measurements on the junction motion suggest that the long-time diffusive motion of the junctions is approximately proportional to the stress tensor. To incorporate the relaxation due to this mechanism into the theory we can put fl,, # 0 in the evolution equation for the active segments. Eq. (4.5) and obtain

~t ~QQ),, = -ilu (QQ)~, + h,l.,, ( Q Q ) ~ - hj.a~QQ),,

(4.15)

where we have taken % = 0. To simplify the results we furthermore assume that the dangling segments relax immediately to the equilibrium state:

~QQ)a=

nakT H

1.

(4.16)

Using the Kramers expression for the stress tensor, Eq. (4.6), we arrive at the following constitutive equation for the network: 2 r5 %.,,4 2+fl,,r2c$t%.~,-

2

2+fl,,rn,kT1.

(4.17)

Comparing this expression with the result of the corresponding affine model which is obtained by putting ft, = 0, shows that the time constant is reduced by a factor p = 2/(2 + flS). The viscosity is reduced by a factor p2 and the primary normal-

324

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303-334

stress coefficient by a factor p3. As can be seen in Fig. 14 this behaviour is indeed observed in many simulations. A problem in this non-aNne model is that is fails to describe the contribution of the elastic forces to the pressure. It is found in the simulations that the elastic forces generate a contribution - n , k T t%the pressure, whereas in the model this is reduced to - p n , k T . This is probably caused by the fact that we only modelled the long-time behaviour of the network and ignored the fast processes that also contribute to the modulus. This point will be discussed in detail in the next section. 4.2. Description in terms o f the relaxation spectrum

For the affinely moving system we expect that there are two relaxation times in the spectrum: a short relaxation time due to the relaxation of dangling segments and a long relaxation time related to the average lifetime of a segment in the network. From the analytical theory it followed that coupling of the relaxation times in fact causes a decrease of the short relaxation time and an increase of the long relaxation time. The analytical results are shown as the solid line in Fig. 5. The modulus of the system G = G(0) is given by G = E i Gi = n k T . In the measurements of the relaxation modulus we found, in accordance with expectations, that the modulus did not change when we moved from an a n n e to a non-aNne system. In a non-aNne system however the junctions are free to move around their affine path. This means that the number of modes in the spectrum increases dramatically since now there are many modes which are due to relaxations of the network structure. Especially, in analogy to the Rouse spectrum for a linear spring, we expect to find many fast modes related to local reorganisations of the network. Since the sum of

'I

5

0

2-

60

0 t"0 r-

.-"""

I 0.14

......

i

8

..••""9 ...'"

E if) co e~

10

E 0.01 0.1

viscosity ratio

Fig. 14. Ratio between affine and non-affine value for the time constant (solid line with slope 1/'2), zero-shear viscosityand zero-shear first normal-stress coefficient (dashed line with slope 3,.'2)for systems 1-10.

B.H.A.A. tart den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

325

2 ,

16! g

L,---".- ---

1.4i

--

o8i

....

06~

="--

04

lO0-node random network i

0.2 ~ 0 (

0

',1~,~

©

500-node random network 4096-node cubic network

011 012 ola 014 ols 016 07 sequence number i/N

018 019

Fig. 15. Spectrum of relaxation times for a random and a cubic network.

the moduli still has to add up to n k T we arrive at the conclusion that the value of the modulus associated with the slowest mode has to decrease which would result in significant lowering of the stress. To illustrate this phenomenon we analysed the relaxation spectra of a permanent random network and of a permanent cubic network. In analogy to the Rouse model the spectrum is obtained from the eigenvalues of the equivalent of the Rouse matrix. The resulting spectra are given in Fig. 15. If we assume that for the last modes of the simulated networks it is not important whether the network is permanent or not, we can use this spectrum to describe the fast modes in the network. In addition there will be a slow mode related to the destruction process. In the simulation, the separation between the time scales is limited for computational reasons so there will be a slight coupling between the slow mode and the fast modes. In a real gel, the time scales can be separated much more and we expect the effect of coupling to be negligible in practice. In Fig. 5 we used this type of spectrum and it can be seen that an accurate description of the relaxation modulus can be obtained this way. There is a free parameter however which determines the fraction of the modulus related to the fast modes compared to the contribution to the modulus of the association/dissociation process. In Fig. 1 and 2 the predictions for the transient viscosity and normal stress difference are compared with simulation results and again a good agreement is obtained.

5. Shear thickening As mentioned in the introduction, shear thickening phenomena are often encountered in associative polymer solutions. In the literature it is either attributed to

326

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303-334

stretching of elastic polymeric segments into the non-linear elastic regime, or to a change from intramolecular to intermolecular associations due to the unfolding of a polymer chain at higher shear rates [13]. Also, it has been suggested that shear thickening is the result of shear induced coagulation of free chains to the network [11]. In our simulations, neither of these explanations can be valid, because a linear spring force law is used, and the number of segments that are loop-connected to themselves is insignificant; typically in the order of 1% of the total number of dumbbells. Also the change of the composition of the system, i.e. a change of the ratio between active and dangling segments can only partially explain the amount of shear thickening observed in our simulations. Significant shear thickening is found to occur if a dangling segment has a high probability of reassociating to the network before it is fully relaxed to the equilibrium state. In Fig. 16 it can be seen that the onset of shear thickening is found to occur at a constant value of the shear stress, independent of model parameters like the rate of dissociation, and functionality or trapping radius of the nodes. The explanation for the shear thickening in our simulations is as follows. When a long active segment becomes dangling it will retract its free end faster than a short segment will do, since the elastic force in a long segment is higher. In a given time interval, a long retracting free end will therefore explore a larger volume of the simulation domain than a short one. This leads to a relative increase of the association probability for longer segments. In Fig. 17 it can indeed be seen that the association probability of a segment increases in proportion to its length. It may be possible to discriminate experimentally between shear thickening due to the above effect and shear thickening due to the stretching of segments into the non-linear regime. To this end we will first briefly describe the line of reasoning of 1O0

o

"',.

I~§

10 " o

1

.

0.01

.

.

.

o

.

.

.

.

i

",

~

.

.

.

~

n2

",

.

.

.

0.1

,

1

.

.

.

.

.

.

.

.

10

shear rate Fig. 16. Viscosity vs. the shear rate for a number of shear thickening systems. The dashed line indicates a constant stress level.

B . H . A . A . r,an den Brule, P.J. Hoogerbrugge

0.7,

J. N o n - N e w t o n i a n Fluid Mech. 60 (1995) 303

! :

0.5

:,J"

~ 1

.,,: ~:J

i ]

/x/.,v,'~~x~'

0.4

to t~

0.3

o

327

I

0.6

o

334

I,

0.2] 0.14

I I

01 0

1'0

2'0

3'0 Qx

4~0

5'0

60

Fig. 17. The association probability of a dangling segment vs. its length.

Marrucci et at. [12]. They assumed that a segment that detaches from an aggregate can only partially relax its extended conformation since it is soon captured by the network again. This is indeed confirmed in our simulations. The consequence of this partial relaxation mechanism, as opposed to the immediate relaxation which is commmonly assumed in network theory, is that above st certain critical shear rate the majority of the chains will be stretched close to their maximum extension. At these shear rates the stress becomes constant. The critical shear rate at which the segments become fully stretched is estimated using the following result from elastic dumbbell theory: ( x -~) = (xg)(1 + 2Zm)'-),

(5.1)

where ,(,, is the effective relaxation time of a chain which is significantly longer than the Rouse time since a relaxing chain is continuously being recaptured by the network. ( x 2) is the mean square projection of the chain dimension in the flow direction and (Xo) = Nb2/3 is the equilibrium value. The parameters N and h refer to the number of links and their length in a Kramers model for the polymer molecule. The contour length of a chain equals Nb so in order for the chain to become fully stretched the shear rate should be of the order of :',,i,= v/N/2m •

5.2)

The viscosity curve has a maximum at the critical shear rate and is estimated to be higher by a factor E/kT than the zero shear rate viscosity, where E is the energy of the barrier for dissociation. In our simulation we expect the onset of shear thickening as soon as the end-to-end distance of a chain becomes significantly larger than the equilibrium value. When this is the case, the elastic retraction of the free end of a dangling segment becomes

328

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

significant to the Brownian motion and the association probability will increase. Using Eq. (5.1) the critical shear rate for this mechanism to become important is estimated to be ~"cri,= 1/;tm,

(5.3)

which can be considerably lower than the prediction of Marrucci et al. [12]. An experiment on H E U R - A T by Annable et al. [18] shows onset of shear thickening at a shear rate of 2 m ? = 0 . 3 which seems to support our predictions. Another argument which seems to support the mechanism proposed here is that in neutron scattering experiments on labelled chains by Pedley et al. [19] no appreciable stretching of the chains was observed at the onset of shear thickening. It can be expected that in many practical systems both mechanisms discussed above are active and careful experiments with varying molecular weights between the end groups are clearly needed. In agreement with existing theories and experiments, shear thickening sets in at a fixed value of the shear stress rather than the shear rate. This critical shear stress is also the same for affine and non-affine simulations. This can be understood by realising that the shear stress is a direct measure of the deviation of the configuration distribution from equilibrium.

6. Conclusion A Brownian Dynamics simulation of a complete network of end-linked chains enabled us to investigate the influence of the topology of the network on stress relaxation after creation or destruction of a network segment. For the affine case, where the motion of the junction is prescribed, an exact analytical theory could be derived. Also, the composition of the network in terms of simulation parameters is well understood. Measurement of the motion of the junctions showed that the average motion is affine, but that there is a diffusion of the junctions around the affine path on a long time scale. In a first approximation, the diffusion tensor describing this process appears to be proportional to the stress tensor, but the oil-diagonal component in shear flow does not follow this proportionality. In order to arrive at a correct constitutive equation for the non-affine model this process needs to be investigated further. A comparison of the response of affine and non-affine network indicated that the modulus associated with the long relaxation times in the spectrum is significantly lower in the non-affine case. This can be understood from the fact that the sum of the moduli has to add up to nkT. In the affine case there is only one fast mode in the spectrum which is related to the relaxation of dangling segments. In the non-affine case, however, there are many fast modes due to small scale motions of the network which are suppressed in the affine model. Since these fast modes all contribute to the total modulus, the modulus associated with the dominant slow modes is reduced. It turned out to be possible to get an accurate description of the

B.H.A.A. i'an den Brule, P.J. Hoogerhrugge

J. Non-Newtonian HuM Mech. 60 (1995) 303 334

329

relaxation modulus by a superposition of the spectrum of a permanent network (for the fast modes) and the relaxation process due to dissociation of the network (the slow mode). Simulation of a polymer gel consisting of linearly elastic chains showed a marked shear thickening effect which can not be explained by the current theories t\~r the rheological behaviour of gels. The shear thickening in these simulations is caused by the fact that the association probability of a dangling segment increases with its length. The predictions based upon this concept appear lo be consistent with currently available experimental results. To explore this area further it is worthwhile to perform experiments with end-linked polymers with varying molecular weight. In this paper we limited ourselves to gels with Gaussian chains and fixed dissociation rates. The simulator, however, is capable of handling finitely extensible chains and dissociation rates that are dependent oil the forces transmitted through the chain segments. This will be the subject of future investigations.

Appendix A: Association/dissociation statistics of nodes

The association/dissociation statistics of beads and nodes can be rather well predicted using some basic assumptions. Let us first recall the basics of these processes. A bead, positioned at one end of an elastic dumbbell, can associate to one node and can have two stales: ,fi'ee and connected. A node can connect to one or more beads, up to some m a x i m u m f ..... . An empty node is repositioned to a random new position in the simulation domain every time step. When a node is connected to a bead it moves with the bead (or beads) and every time step it has a fixed probability of dissociating from the attached beads. It is assumed that the distribution of free beads in the simulation domain is homogeneous and uncorrelated with the distribution of nodes. Furthermore, it is assumed that the motion of the free beads is dominated by Brownian motion. The average number of free beads captured by an empty node when it jumps to a new position is given by the concentration of free beads C~ times the capturing volume 4/37rd 3, where d is the trapping radius. The transition rate of empty to partially filled nodes thus becomes C I 4 7rd3"

(A1 }

<~~}.l - ($t 3

The capturing of free beads by partially filled nodes is dominated by the diffusive motion of the free beads (this is shown in Appendix B). The rate of association is approximately <~1;I

- ~

=

4rcdDCl,

(A2)

330

B.H.A.A. van den Brule, P.J. Hoogerbrugge /J. Non-Newtonian Fluid Mech. 60 (1995) 303 -334

where D = k T / ( is the diffusion coefficient of free beads. Full nodes can not capture more beads: gfm~,~, ..... , = 0.

(A3)

The dissociation rate of any node is fixed, independent of the number of beads attached to it, g~:o = I/r,

(A4)

where T is the average life time of a node. In steady state, the number of empty, partially filled and full nodes is stationary and therefore the following balances should hold: go.l No = ~ gl:0(N - No),

(A5)

f

g/ 1,1N~ 1 = (g1:o + g~;l+ I)Nr,

(A6)

where N~ is the number of nodes with f beads connected and N is the total number of nodes in the simulation domain. For the number of nodes we have N=

~ Nr l-

(A7)

0

For the concentration of free beads we can derive the following equation - t ~ 0 f N j ],

(18)

where V is the volume of the simulation domain and Nb is the total number of beads present in V. The above set of equations can be solved numerically to find Nj and C / f r o m the simulation parameters. The number density of active, dangling and free segments, n,,, n~ and nr, can be calculated from the fraction of free beads Ot = C j V / N b. If the probability for a bead to be free is qS/, then the probability for a segment to have two free beads equals ~bt~ and the number density of free segments thus becomes n r = q~n.

(A9)

Similarly, the probability for a segment to have two beads connected to a node is given by (1 - q~t)2 and the number density of active segments equals n,, = (1 - ~t)2n.

(A 10)

For the number density of dangling segments we may write n,t = 2~br(1 - ~br)n.

(A 11 )

In Table 1 the predicted results are compared with measurements in a number of simulations. It can be seen that there is good agreement.

B.H.A.A. t:an den Brule, P.J. Hoogerbrug,ge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

331

Appendix B: Association of beads to partially filled nodes A partially filled node can meet a free bead through three different mechanisms that play a role in our simulations: by diffusion, by convection and by retraction of a free end through the fluid. We will discuss each of these processes to assess their relative importance.

D(ffusion It is assumed that the motion of a partially filled node, which is connected to the network, is negligible compared to the motion of a free bead. Following the classical Smoluchowski theory for coagulation we furthermore assume that the concentration profile C(r) of free beads around a partially filled node has become time independent in an averaged sense. The number of particles reaching the surface of the static spherical trap of radius d can then be found by solving Fick's diffusion equation and the mass balance, subject to the boundary conditions C(d) = 0 and C ( o c ) = Q ,

j= -DVC,

(B1)

~C --=V ~t

(B2)

.j=0

where j is the flux of free beads per unit area, C is the local concentration of tree beads and D = kT/( is the diffusion coefficient. The solution for the flux becomes

j = - DC~d/r 2

(B3)

and the number of free beads entering the spherical trap by diffusion is given by

4rcdDQ..

g'dm- = ,]', 47~d2 =

(B4)

A simple correction for the motion of the nodes would be to introduce a diffusion coefficient for the motion of a node with functionality f: D ' = D + Dnod~ = D(1 + l/J).

(B5)

Usually, the effect of this correction turns out be small.

Convection When the fluid is sheared free beads will be convected into the spherical trap around a node. This number of beads entering the trap in unit time is given by

g~he~r= ct f

f

v . n dA =4 d3~,Ct.

(B6,

(v-n
Retraction A bead at the free end of an extended segment will be dragged through the fluid

332

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

with a velocity that can be derived from a force balance, where we now ignore Brownian motion: F ~c) + F ~v) = - H Q - (v = 0,

(B7)

where v is the relative motion of the bead with respect to the surrounding fluid. The number of partially filled nodes encountered by a free bead in unit time therefore becomes h = v;~d2rtpf = 7td2H(Q ~lnpf/(,

(Bg)

npf is the number density of partially filled nodes. The association rate of free beads due to retraction to a particular node can now be expressed as a function of the length of the connector vectors averaged over all free beads where


gretr = =d2H< Q ) , C t / ( .

(B9)

In dimensionless form, where we use 2 to dimensionalise the time and ( x / ~ / H ) for the length scale, the following results for the various contributions to the association rate are obtained: oadi~r= 4~3/5Q, ^

4 A3 "*

(B 10)

gshe,r=~d ~Q,

(Bll)

A ! gretr = g /'C~2( 0 ~ f ~ .

(BI2)

With the largest values for the trapping radius and the shear rate, d - - 0 . 3 and = 3, the association rates become gdiff ~, ~ Cf,

(B13)

g~h~a~~ 0.01 ~I,

(B14)

gretr ~' 0 . 0 7 ( Q ) r G ,

(B15)

where we used /5 = 1/4. This comparison reveals that the diffusion process is dominant, except when { Q ) I becomes very large. (In equilibrium ( Q ) l = x f l S . ) The probability that a free bead associates to a node in a unit time interval can be written as P = gnpf/C I.

(B16)

The association probabilty has been measured directly in simulation system 4 where d = 0 . 1 7 5 and ~ = 1. The observed value for the concentration of free beads is C'j= 0.774 and the number density of partially filled nodes was measured to be h p f = 0.34. According to the above theory the association probability in a dimensionless time unit can be written as P ~ npf{4rrd/5 + ¼~rd2(O ~1}

B.H.A.A. I~an den Brule, P,J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

0.t9 + 0.008(0)~

333

(B17)

This is in reasonable agreement with the simulation result (see Fig. 17) of P = 0.27 + 0.005
(B18)

References [1] L.R.G. Treloar, The Physics of Rubber Elasticity, Oxford University Press, Oxford, 1975. [2] M.S. Green and A.V. Tobolski, A new approach to the theory of relaxing polymeric medial, J, Chem. Phys., 14 (1946) 8 0 9 2 . [3] M. Yamamoto, The viscoelastic properties of network structure-l. General formulism. J. Phys. Soc. Jpn., 11 (1956) 413 421; The viscoelastic properties of network structure-2. Structural viscosity, J. Phys, Soc. Jpn., 12 (1957) 1148 1158: ]'he viscoelastic properties of network structure. 3. Normal stress effect (Weissenberg effect), J. Phys. Soc. Jpn., 13 (1958) 12(10 1211. [4] A.S. Lodge, A network theory of flow birefringence and stress in concentrated polymer solutions, Trans. Faraday Soc., 52 (1956) 120 130: Constitutive equations from molecular network theories lk)r polymer solutions, Rheol. Acta, 7 (1968) 379 392. [5] R.B. Bird, C.F. Curtiss, R.C. Armstrong and O. Hassager. Dynamics of Polymeric Liquids. Vol. 2. Kinetic Theory, J. Wiley & Sons, New York, 1987. [6] R.G. Larson, Constitutive Equations for Polymer Melts and Solutions, Butterworth Publishers, Stoneham, 1988. [7] N. Phan-Thien and R.I. Tanner, A new constitutive equation derived from network theory. J. Non-Newtonian Fluid Mech., 2 (1977) 353 365. N. Phan-Thien, J. Rheol., 22 (1978) 259 283. R.[. Tanner, Some useful constitutive equations with kinematic slip hypothesis, J. Non-Newtonian Fhfid Mech.. 5 (1979) 103 112. [8] F. Petruccione and Biller P., A numerical stochastic approach to network theories of polymeric fluids, J. Chem. Phys., 89 (1988) 577 582. W. Hermann and F. Petruccione, Quantitative rheological predictions of a transient network model of Lodge Yamamoto type: Simple and multiaxial elongational flow, J. Rheol., 36 (1992) ~461 1476. [9] L. Leibler, M. Rubinstein and R.H. Colby, Dynamics of reversible networks, Macromolecules, 24 ( 1991 ) 4701 4707. [10] F. Tanaka and S.F. Edwards, Viscoelastic properties of physically cross-linked networks. Transient network theory, Macromolecules, 25 (1992) 1516 1523. [11] S.-Q Wang, Transient network theory for shear thickening fluids and physically cross-linked systems, Macromolecules, 25 (1993) 7003 7010. [12] G. lVlarrucci, S. Bhargava and S.L. Cooper, Models of shear thickening behaviour in physically cross-linked networks, Macromolecules, 26 (1993) 6483 6488. [13] T.A. Witten and M.H. Cohen, Cross-linking in shear thickening ionomers, Macromolecules, 18 (1985) 1915 1918. [14] K.H. Ahn and H. Osaki, A network model for predicting the shear thickening behaviour of a poly(vinyl alcohol)/sodium borate aqueous solution, J. Non-Newtonian Fluid Mech.. 55 (1994) 215 227: Mechanism of shear thickening investigated by a network model, ]. Non-Newtonian Fluid Mech., 56 (1995) 267 288. [15] E.J. Hinch, Application of the Langevin equation to fluid suspensions, J. Fluid Mech.. 72 (1975). [16] D.L. Ermak and J.A. McCammon, Brownian Dynamics with hydro-dynamic interactions..1. Chem. Phys,, 69 (1978) 1352-1360. [17] A. Greiner, W. Strittmatter and J. Honerkamp, Numerical integration o|" stochastic differential equations, J. Star. Phys., 51 (1988) 95 108. [18] T. Annable, R. Buscall, R. Ettelaie and D. Whittlestone, The theology of solutions of associatir~g polymers: Comparison of experimental behavior with transient network theory, J. Rheol., 37 (19931 695 726.

334

B.H.A.A. van den Brule, P.J. Hoogerbrugge / J. Non-Newtonian Fluid Mech. 60 (1995) 303 334

[19] A.M. Pedley, J.S. Higgins, D.G. Peiffer and E. Staples, Single chain dimensions in ionomer solutions under quiescent and shear conditions as determined by small angle neutron scattering, Polym. Commun., 30 (1989) 162 165.