- Email: [email protected]
Numerical study of the rheology of rigid fillers suspended in long-chain branched polymer under planar extensional flow
J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Contents lists available at ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm
Numerical study of the rheology of rigid fillers suspended in long-chain branched polymer under planar extensional flow M. Ahamadi ∗ , O.G Harlen Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, United Kingdom
a r t i c l e
i n f o
Article history: Received 1 September 2008 Received in revised form 3 December 2009 Accepted 7 January 2010
Keywords: Finite elements (65M60) Fluid mechanics Suspensions rheology Planar extension Periodic boundary conditions Viscoelastic fluid
a b s t r a c t We report a detailed numerical study of the rheology of two-dimensional rigid fillers suspended in branched polymer melt under planar extensional flow. The polymer melt is modelled using the pom-pom constitutive equation. The numerical method uses a finite element solution of the flow in a unit cell within the self-replicating lattice for planar extensional flow identified by Kraynik and Reinelt [Int. J. Multiphase flow 18 (1992) 1045]. This is implemented using a quotient space representation that maps all points in space onto points within the unit cell. We show that the Kraynik and Reinelt cell allows simulations of suspensions under planar extensional flow to be conducted to large strains in a truly periodic cell. The influence of both the pom-pom parameters and the particle area fraction on the rheology of the suspension are investigated. We find a reduction in the degree of extension-rate thickening with the increase of both particles concentration and Weissenberg numbers in agreement with published experimental and numerical findings on other viscoelastic models. This reduction is found to be due to flow disturbance created by the particles which disrupts the alignment of backbone tube segments with the extensional axis. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Understanding the extensional rheology of filled polymer melts is important for the prediction of material performance in many commercial applications. One of the most important classes of such fillers are solid particles, typically glass or carbon particles, suspended in a polymeric matrix. The rheology of polymer blends under extensional flow has been studied by many investigators [22]. For long-chain branched polymers the effective viscosity of the polymer during the start-up (defined as the ratio of the stress difference to the imposed strain rate) grows above the level predicted by linear viscoelasticity. This phenomenon, known as strain hardening, is associated with the stretching of sections of the molecular chains, which occur at strain rates that exceed the effective Rouse or stretch relaxation time of the chain. In entangled long-chain polymer melts this relaxation time is controlled by the constraints on the motion of the branch points due to entanglements. Strain hardening is a desirable feature that improves melt processability in blow molding, foaming, and melt spinning [21]. Thus the effective extensional viscosity and the degree to which it shows strain hardening at high strain rates is an important rheological measurement. For filled polymers it has been observed experimentally [15,16,20] and numerically [6] that the degree of strain hardening
∗ Corresponding author. E-mail address: [email protected] (M. Ahamadi). 0377-0257/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2010.01.002
generally decreases with increasing particle concentration. A precise explanation of this phenomenon is still lacking, although a number of hypotheses have been suggested. Kobayashi et al. [16] argued that the presence of the particles perturbs the deformation of matrix polymer chains around them and thus causes the suppression of strain hardening. In their numerical study of concentrated particle suspensions, D’Avino et al. [6] suggested that the reduction in strain hardening with increasing Weissenberg number is caused by the shear flow between the particles. The effect of adding fillers on the extensional properties depends upon the type of polymer matrix used. Le Meins et al. [20] measured the elongational viscosity of monodisperse polystyrene spheres suspended in a viscoelastic polyisobutene matrix. They observed an increase in the extensional viscosity with increasing filler volume fraction and strain rate. At large strain rate, strain hardening occurs in the pure polymer but is suppressed by adding particles. Kobayashi et al. [16] studied the effect of glass beads on the elongation flow of high density polyethylene (HDPE). Their experimental results showed weaker strain hardening with increasing strain rate in both the pure polymer, as well for the suspension. The suppression of the strain hardening by the beads is stronger at high strain rate where the strain hardening of the pure polymer is weak. In this paper we present a detailed numerical study of the rheology of fillers under planar extensional flow. We focus our attention on the case of hard and non-Brownian circular particles suspended in long-chain branched polymer melts modelled using the pompom constitutive model. This class of polymer was chosen because
282
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
long-chain branched melts exhibit strain hardening in extensional flow, whereas linear polymers are generally strain softening in both shear and extensional flow. The pom-pom model developed by McLeish and Larson [19] contains the key features of branched polymers, i.e., it models chain segments whose reptation is inhibited by the presence of branch points. The original single mode model is able to correct qualitative behaviour of branched polymer in both shear and extension (and multimode versions are able to capture this rheology quantitatively [13,4,27]). In addition, the parameters in the models are directly linked to the microstructure of the polymer melt. The idealised pom-pom molecule has a single backbone with multiple branches emerging from each end. In extensional flow the backbone becomes stretched if the extensionrate is faster than the stretch relaxation time, s of the backbone. However, since the tension forces at the branch-point must balance there is an upper-bound to this extension, q, equal to the number of arms connected to the branch-point. This upper-bond on extension limits the degree of extension hardening. In planar extensional flow, as shown in Fig. 1, fluid particles separate exponentially in time in the direction of extension. Thus for a complex fluid any finite computational domain will usually become very highly distorted within the characteristic time for the fluid to reach steady state. Consequently, in order to be able to reach a steady state, it is desirable to use a self-replicating lattice that recovers its original geometry after a period of strain. This is not straightforward to simulate since the dynamical motion of extensional flow involves a change in the shape of the unit cell within the computational domain. Kraynik and Reinelt [17] found a set of self-replicating lattice for planar extension flow, which we will use to impose the periodic structure of a two-dimensional computational domain ˝(t). The numerical method is based upon the Lagrangian–Eulerian approach of Harlen et al. [9]. In this application this approach provides a natural way to implement the particle motion within the Kraynik and Reinelt lattice, as well as enabling the constitutive equation to be solved in the natural co-deforming frame of the fluid. 2. Computational domain We consider a two-dimensional unbounded domain, ˝∞ (t), containing non-Brownian hard circular particles suspended in viscoelastic fluid under a planar extensional flow. In order to make the computations tractable we approximate the fluid microstructure by a spatially periodic fluid microstructure. Therefore we solve the flow in an initially rectangular unit cell of size H × contai-
ning N particles, ˝(t) that when periodically replicated covers the infinite domain. Under the action of the velocity gradient, the unit cell within the lattice deforms and in order to be able to continue to large strains it is desirable to use a self-replicating lattice that recovers its original geometry after a period of strain. In shear flow, this can be achieved using the Lees–Edwards [12,18] lattice in which cells slide across one another in layers parallel to the flow direction. However, in a planar extensional flow (see Fig. 1), a rectangular lattice aligned along the principle flow axes will not periodically replicate its structure. A lattice of this type, suggested by Heyes [10] is used in the simulation method developed by Hwang and Hulsen [11] for planar elongational flow of suspensions in a Newtonian fluid. However, because the aspect ratio of the unit cell increases exponentially in time their simulations are limited to moderate strains. D’Avino et al. [5] have recently developed an alternative fixed grid method for simulating extensional flow. This does not use periodic boundary conditions but uses a three layer domain in which particles leaving the domain are reintroduced at random positions on the inlet boundaries. Kraynik and Reinelt [17] found a set of self-replicating lattice for planar extension flow by using the condition of reproducibility of a two-dimensional square lattice L(t) = n1 p(t) + n2 q(t), where p(t), q(t) are linearly independent basis vectors, n1 , n2 integers. They demonstrated that the square lattice with the minimum period for replication was one where the√sides are aligned at an angle of approximately = tan−1 (2/(1 + 5)) = 31.7◦ to the principal axes √ of the flow, see Fig. 1b, with the Hencky strain p = log((3 + 5)/2) = 0.9624. The important outcome of the Kraynik and Reinelt’s findings is that for this initial orientation the square lattice will reproduce itself after a strain of p in the sense that vertices of the original lattice move to other vertex points of this lattice. The period of reproducibility of the Kraynik and Reinelt cell depends only on the Hencky strain and so this lattice may be used for a variable strain-rate simulation provided that the principle axes of strain remain the same. In our previous paper [1], we detailed how to implement the Kraynik and Reinelt lattice structure for suspension in Newtonian fluid by using an appropriate image system that map positions in the infinite domain ˝∞ (t) to positions in the computational domain, ˝(t). In this paper we show how our method can be extended to viscoelastic fluids. 3. Governing equations The flow is assumed to be incompressible, isothermal and inertialess, so that the equations of conservation of mass and
Fig. 1. (a) Particle paths in planar extensional flow. (b) The unit cells in the critical lattice.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
momentum are given by,
∇ · =0
(1)
∇ ·u=0
where u is the fluid velocity and is the stress tensor. The constitutive model we will consider is the pom-pom model of McLeish and Larson [19]. This model is obtained by considering a melt of pom-pom molecules formed by connecting two q-armed star polymers with a backbone chain. The presence of the branch points at the ends of the backbone chain inhibits its motion along its tube, so that stretch and orientation relaxation times of the backbone, s and , respectively, are controlled by the relaxation of the star arms. We shall consider a single pom-pom mode describing the stretch and orientation of the backbone together with a Newtonian stress term for the fast relaxing modes associated with the arms so that the stress is of the form = −pI + 2s D + P ,
(2) P
where p is the pressure, 2s D is the Newtonian stress and is the backbone stress contribution. In this study, we will adopt the differential version of the original pom-pom model, but with the modification introduced by Blackwell et al. [4]. P = 3G2
A . Tr(A)
(3)
Here the auxillary tensor A satisfies an upper-convected Maxwell equation ∇ 1 A = − (A − I),
(4)
283
around the surface of the particle l as a Lagrange multiplier to force the fluid inside particle l to behave as a rigid solid. Here n is the outward unit normal to the boundary of particle Pl . The angular velocity ωl and the translational velocity Ul are found from the conditions of no net force and torque on each particle,
fl ds = 0,
l = 1, . . . , N
(8)
∂Pl
(x − Xl ) × fl ds = 0,
l = 1, . . . , N
where fl is the traction force density vector on particle l. The motions of particles are described by the kinematic equations dXl = Ul , dt
Xl |t=0 = Xl,0 ,
l = 1, . . . , N
1 ij = V
D e2(−1)/q = A : ∇u − ( − 1) Dt TrA s
where Vf is the volume of the fluid phase.
up to the maximum stretch, q. For this model we can define four dimensionless groups: the non-dimensional elastic modulus c = (3G/s ), the number of arms q, the Weissenberg number based on the orientation relaxation time We = ˙ and the Weissenberg number based upon the stretch relaxation time Wes = Wes /. Note that c is independent of ˙ and so remains fixed for a given material. In this paper we will show results for c = 9 and Wes = (1/2)We and vary q and We. The ratio We/Wes is proportional to the entanglement length of the backbone and a value of 2 would correspond to pom-pom molecules where the backbone length is equal to 5 entanglements and makes up 37% of the molecular weight of the polymer. Varying this ratio does not qualitatively affect the results. The extensional flow is imposed through the boundary conditions of the computational domain, which will be discussed in Section 4.2. 3.1. Particles The filler particles are assumed to be rigid and subject to zero net force and couple. Assuming that there is no slip between the filler particles and the fluid matrix, the fluid velocity at the surface of particle l is given by u = Ul + ωl × (x − Xl ),
l = 1, . . . , N
(6)
where Ul and ωl = ωl zˆ are the unknown velocity and angular velocity of particle l (with zˆ the unit vector in the zˆ direction). To enforce these boundary conditions we impose a surface force density fl (s) = n,
(7)
(10)
d l = ωl , l |t=0 = l,0 , l = 1, . . . , N (11) dt where l is the angular rotation of particle l. Note that the kinematic Eqs. (10) and (11) are decoupled from Eqs. (6), (8) and (9). Hence these equations of motions can be solved separately by using a standard ordinary differential equation method as we will describe later in Section 4.2. The quantity measured in rheological experiments is the spaced averaged stress of the suspension. We shall refer to this as the “bulk” stress to distinguish it from time-averaged quantities. Batchelor [2] derived an expression for the bulk stress in a suspension as sum of the stress in the fluid phase plus a contribution from the surface forces on the particles, in the form
where is the orientation relaxation time, while the backbone stretch, , is given by (5)
(9)
∂Pl
Vf
1 ij dVf + 2
[(xi − Xi )fj + fi (xj − Xj )ds]
(12)
∂Pl
4. Computational method 4.1. Finite element formulation We shall use the numerical method based on the split Eulerian–Lagrangian technique of Harlen et al. [9]. This uses a decoupled method in which the constitutive Eqs. (4) and (5) are solved as ordinary differential equations in the co-deforming (Lagrangian frame) within a separate step from the calculation of the fluid velocity. The fluid velocity and pressure are calculated for the current values of microstructural variables A and and positions of the particles. This is achieved by solving the Stokes Eq. (1) with Eqs. (8) and (9), and the no-slip boundary condition (6) via a finite element method. We then step forward in time to find the new values of A, and particle positions by solving the constitutive Eqs. (4) and (5). To render the Eqs. (1), (6), (8) and (9) in finite element form we shall adopt the combined equation of motion of Glowinski et al. [8] to derive the weak form of the equation for the solid–liquid mixture. The natural combined velocity space for the fluid and particle velocities is given by 2
V = {((v, Vl , l )|v ∈ H1 (˝)) , Vl ∈ R2 , l ∈ R, v = Vl + l zˆ × (x − Xl ) in Pl (t), and v biperiodic on ∂˝}.
(13)
In the distributed Lagrange multiplier method the extended weak form for whole domain can be obtained by removing the
284
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
constraint Eq. (6) from the velocity space and enforcing it weakly as a side constraint. This is done by introducing a Lagrange multiplier (fl (s) in our case), which can be interpreted as the traction force required to maintain the rigid-body motion of particle Pl (t). In our methodology, the constraint Eqs. (8) and (9) are used to determine the unknown Ul and ωl and therefore they must be incorporated into the final weak form.
where p(t), q(t) are the base vectors of the Kraynik and Reinelt cell and W = {(a, b) ∈ R × R such
Q : R2 →
−p∇ · w dV +
˝
∇ u : ∇ w dV − ˝
N
l=1
fl · w dS
(14)
˝
q∇ · u dV = 0,
−
(15)
˝
l · (u − (Ul + ωl × (x − Xl ))) dS = 0,
(16)
∂Pl (t)
fl dS = 0,
(17)
fl × (x − Xl ) dS = 0,
(18)
∂Pl (t)
∂Pl (t)
where l = 1, . . . , N and w, q, l are respectively the variations for the fluid velocity, pressure, and traction force. We use triangular P1 − P1 elements, which have linear continuous pressures and velocities, and a continuous linear interpolation for the force density fl around the particles to discretise the weak form [Eq. (14)–(18)]. The viscoelastic stress P is piecewise constant, P0 , over elements. The element combination P1 − P1 for the velocity and pressure is unstable to spurious pressure modes. These are removed by introducing a pressure stabilisation in which the weak form of continuity Eq. (15) is replaced by
q∇ · u dx − ˇ
− ˝
is the quotient map allowing us to identify all images. From the general form of the quotient map, Eq. (21), one can derive an explicit relation between a point in ˝(t) and its images in ˝∞ and between the velocities and other quantities held at these points. Thus after a strain of = ˙ t in planar extensional flow the images of a point (xlattice , ylattice ) are given by ximag = xlattice + me p01 + ne q01 H
−
(21)
ximag → Q(ximag ) = [xlattice ]
∂Pl (t)
P : ∇ w dV,
=−
R2 W
h ∇ q · ∇ p dx = 0, 2
(19)
˝
where h2 is twice the area of the triangle and ˇ is a positive constant. Following Silvester and Wathen [26], we choose ˇ = 0.025 as this gives the optimal convergence rate for the Stokes problem. For a given finite element mesh, the discretisation of Eqs. (14)–(18) leads to a symmetric indefinite linear system of algebraic equations. These are solved using a preconditioned conjugate residuals method with block preconditioner of the form suggested by Silvester and Wathen [26]. 4.2. Implementation of the Kraynik and Reinelt lattice The extensional flow is imposed through the periodically selfreplicating Kraynik and Reinelt lattice using the method described in ref. [1]. We generate an initial biperiodic mesh for the computational domain, ˝(t0 ), in which the nodes on opposite edges are identified as being equivalent points. This is done using the “Triangle” grid generation program [24] to generate an initial grid in the unit cell and then relabelling boundary points so that equivalent points have the same label. This equivalence can be extended further by noting that every point in the computational domain has an equivalence class of image nodes in the infinite domain, ˝∞ . Let ∼ be an equivalence relation on V = R2 defined by x∼r ⇔ ∃ m ∈ Z, n ∈ Z such
(20)
yimag = ylattice + me− p02 + ne− q02 H,
(22)
where the integer pair (m, n) denote the shift in the periodic image. Since the finite element matrices depend only upon the relative displacement of nodes and not their absolute position we define the displacement dxlattice of nodes on an element to be the displacement between images that lies in the space dxlattice = dpp(t) + dqq(t) −
< dp ≤ , 2 2
−
with
H H < dq ≤ . 2 2
(23)
The external extensional flow is imposed by the condition that the fluid velocity uimag of the (m, n) image is given by the time derivative of the quotient map given in Eq. (21). Consequently, the fluid velocity, uimag is given by uimag = ulattice + ˙ (Tu , Tv ), me˙ t p01 + ne˙ t q01 H
(24) −me−˙ t p02 − ne−˙ t q02 H.
and Tv = where Tu = All other variables, including pressure and the microstructure variables A and have the same value at each image point. A major drawback of using a Lagrangian mesh is that severe deformation-induced mesh distortions develop due to velocity gradients within the fluid, which will degrade the accuracy of the finite element solution. A partial solution to this problem is to retain the nodes as material points, but reconnect them in the way that produces a Delaunay triangulation [7]. The Delaunay triangulation ensures global equiangular triangulation in the sense that it maximises the minimum angle in any triangle [25]. The Lagrangian mesh movement combined with Delaunay reconnection and the use of the equivalence class to relate points with their images automatically reproduces the self-replicating lattice of Kraynik and Reinelt. This is illustrated in Fig. 2 where the evolution of a very simple mesh consisting of eight triangular elements connecting four replicated points and their appropriate images under planar extension is shown. The initial mesh is shown in Fig. 2a with the square unit cell (indicated by the bold lines) oriented at an angle = 31.7◦ to the x-axis. Fig. 2b and c show the lattice after a strain of (1/2)p . The original unit cell (shown with solid bold lines in Fig. 2b) is deformed into a parallelogram, however, we can reform an equivalent square unit cell (shown with dashed lines), by replacing p with p + q. Between Fig. 2b and c we have performed a Delaunay reconnection. Although we have recovered a square lattice it is not the original lattice as it is now oriented at 58.3◦ to the x-axis. Fig. 2d shows the lattice after a strain of p = 0.962424. We can now recover the original square unit cell by choosing appropriate images. A further Delaunay reconnection will reproduce the original grid shown in Fig. 2d.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
285
Fig. 2. A simple example showing how the mesh movement and reconnection replicate the unit cell after a strain of p = 0.9624. (a) The initial lattice with a unit cell made up of eight triangular elements. (b) The lattice after a strain of (1/2)p . The initial unit cell is shown with the solid line, while the dashed line shows an equivalent square unit cell formed by replacing two of the elements with their images. (c) The lattice after a strain of (1/2)p following a Delaunay reconnection. (d) The lattice after a strain of p = 0.9624. The original unit cell and lattice can be recovered by a further replacement of two elements with their images and a further Delaunay reconnection.
4.3. Lagrangian calculation of the viscoelastic stress
we employ the methods of grid improvement detailed in [1] during the calculation.
The value of the viscoelastic stress within each element is calculated from the values of microstructure variables A and , which are updated by integrating Eqs. (4) and (5) forward in time in a frame deforming with the fluid [3]. This is done by moving the vertices of the triangles with their current velocity. Each triangle is, therefore, a material volume (since the velocity varies linearly over the triangle) and the vectors along any two sides of the triangle form a set of base vectors for the frame that moves and deforms with the fluid. The time derivative in this frame is therefore the upper-convected derivative, so that A may be found by solving the ordinary differential equation ˜ dA 1 ˜ ˜ A−I =− dt
(25)
˜ and ˜I are the components of the A and the identity tensor, I where A in triangle coordinates [9]. Similarly, the Lagrangian time derivative (D/Dt) becomes (d/dt) in this frame. The resulting ordinary differential equations are solved by forward time-stepping and then A is transformed back into Cartesian coordinates. The grid is then reconnected where necessary to give a Delaunay triangulation. For the simulations containing particles the Delaunay reconnection is not allowed to reconnect an edge on a particle boundary. However edges on the boundary of the unit cell are reconnected so that the outline of the original cell is lost. As noted earlier, due to the equivalence class it does not matter which image of a triangle is used in the assembly. Once particles are introduced reconnection alone is not sufficient to maintain mesh quality and so
5. Results and discussion In planar extensional flow the bulk extensional viscosity is defined as =
− xx yy ˙
.
For a purely linear viscoelastic fluid this is equal to 4 times the transient shear viscosity. However, in a strain-hardening material it will be a function of both t and the applied strain rate ˙ . For the pom-pom model given in Eq. (3) the extensional viscosity is given by = 4s +
3G 2 c 2 (Sxx − Syy ) (Sxx − Syy ) = s 4 + We ˙
where S = A/Tr(A) is the backbone tube orientation distribution. For large strain rates (Wes > 1) both and Sxx − Syy approach their limiting values of q and unity, respectively, so that /s will grow from 4 to a limiting value of 4 + (c/We)q2 at high strain. This is shown in Fig. 3a for the case of q = 2 and Wes = (1/2)We = 1.75 where it can be seen that steady state is achieved after a Hencky strain of 3. At intermediate strain rates (We > 0.5, Wes < 1) the backbone tubes become fully aligned, but do not stretch so that remains close to unity, while for strain rates for which We 1 we recover the linear viscoelastic limit.
286
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Fig. 3. Transient extensional viscosity for a pom-pom fluid with q = 2 and the stretch Weissenberg number Wes = 1.75: (a) pure fluid, (b) single particle of radius 0.3 (˚ = 0.2827) for the three different meshes given in Table 1.
Table 1 Details of the four finite element meshes used to calculate the extensional viscosity at t = 0 for a circular particle of area fraction ˚ = 0.2827. Mesh
Nb
Np
Area min
Elements, Ne
Average of bulk viscosity
Mesh 1 Mesh 2 Mesh 3
32 44 62
140 196 277
0.0025 0.001 5 × 10−4
1150 2210 4054
14.2439592866843 14.2284206433883 14.2163790630068
The addition of filler particles affects the flow in number of ways. First, since the particles do not deform the average strain rate in the −1 fluid phase increases by a factor (1 − ˚) , where ˚ is the area fraction of the particles and so increases the effective values of the Weissenberg numbers. However, the particles also perturb the flow away from pure planar extension which will rotate the tube segments away from perfect alignment. 5.1. Single particle per cell We first consider the case of a single particle at the centre of computational cell. In Fig. 3b we plot the bulk extensional viscosity for a particle of radius 0.3 (giving a particle area fraction ˚ = 0.2827) suspended in a pom-pom fluid with q = 2 and Wes = 1.75. In order to test the spatial convergence of our numerical method we compare the results from three different uniform meshes whose details are given in Table 1. Here, Nb is the number of nodes along each edge of the unit cell, Np the number of nodes around the particle, and Area min is the threshold value for the area of each triangle
used by “Triangle”. Mesh 2 is shown in Fig. 4. In our previous paper [1] we found that the spatial convergence for Newtonian fluids was of order h2 , and it can be seen that there is very little difference between the results for the two finest meshes. In subsequent calculations we will use meshes similar to mesh 3. Comparing Fig. 3 a and b we see that the bulk viscosity of the filled systems grows in a similar manner to that of the pure fluid, but with a small oscillation superimposed on it. The frequency of this oscillation is equal to twice that of the self-replicating lattice and is due to hydrodynamic interactions of the particle with its images. The maximum extensional viscosity occurs when the particles are arranged in a square lattice, which occurs twice in each period p . The lattice configuration after a strain of (1/2)p is square, but is rotated relative to its initial configuration. By averaging over the period of this oscillation at large strains we can obtain the average bulk viscosity and compare this with the steady state value for the pure fluid. In Fig. 5a we compare the average bulk viscosity for ˚ = 0.237 with the steady viscosity of the pure fluid at a stretch Weissenberg number Wes = 2.25 for varying numbers of arms q. As expected, increasing the number of arms increases the extensional viscosity for both the filled and unfilled materials. For Newtonian fluids the addition of rigid filler particles always increases the extensional viscosity, however, at this Weissenberg number adding particles reduces the bulk extensional viscosity. As noted earlier, since Wes is greater than unity, the backbone segments become fully aligned and stretched in the pure fluid so that they contribute their maximum possible value to the extensional stress. Adding particles disrupts the alignment and the consequent reduction in the poly-
Fig. 4. The computational domain for the single particle problem. (a) The unit cell at time t = 0 in relation to its periodic images. (b) The finite element Mesh 2 used in the calculations described in Table 1.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
287
Fig. 5. The effect of the varying the arm number q on the extensional viscosity of rheology of a suspension at ˚ = 0.237. (a) Plot of the average bulk extensional viscosity against arm number for Wes = 2.25. (b) Plot of the ratio of the suspension average viscosity to that of the pure fluid for Wes = 1, 1.25, 1.75 and 2.25.
meric stress is sufficient to overcome the addition to the solvent stress resulting from the addition of particles. By plotting the ratio of average bulk extensional viscosity of the filled system to the steady extensional viscosity of the suspending fluid we can see more easily the effect that adding particles has on the extensional viscosity. In Fig. 5b we plot this ratio for stretch Weissenberg numbers ranging from 1 to 2.25. The ratio decreases with increasing stretch Weissenberg number showing that the reduction in extensional viscosity occurs at high strain rates where the viscoelastic stress contribution is at its maximum in the unfilled system. For the two highest stretch Weissenberg numbers the ratio decreases with increasing arm number, however at a stretch Weissenberg number of unity the ratio increases with increasing arm number. This Weissenberg number corresponds to the transition between stretched and unstretched backbones, and so the backbone segments are not fully extended in the pure material. This means that in regions of the flow where the local extension rates are higher the backbone segments can become more highly extended than in the pure material. In Fig. 6 we compare the average backbone extension in suspensions of area fractions ˚ = 0.126 and 0.237 to the value of in the pure melt for a 3-armed pom-pom. At stretch Weissenberg numbers below 0.3 (Fig. 6a) the average extension increases with filler concentration. Because the particles do not strain, the average extension-rate in melt phase of a suspension is higher than in −1 the pure material by a factor (1 − ˚) . Consequently, the effective values of the Weissenberg number in the melt phase of a suspension is higher than the value calculated from the bulk deformation rate. Infact for small Weissenberg numbers the results can
be superimposed by shifting the Weissenberg number by this factor. However, at stretch Weissenberg number above 0.3 the effect is reversed (Fig. 6b) with the average backbone extension decreasing with filler concentration. In all three materials the degree of extension increases with We, but in the pure melt, achieves its maximum value for Wes > 2.5, whereas the velocity perturbations in the suspensions prevent every from reaching this value. By examining the ratio of the average backbone extension to the value of in the pure melt in Fig. 7a we can identify a number of different regimes. For very small Weissenberg numbers where is close to unity and the orientation is close to being isotropic, grows faster than due to the increase in the average extension-rate, so that the ratio is greater than unity. As can be seen in Fig. 8 at these low Weissenberg numbers the degree of tube alignment increases with ˚. However at around We = 0.5 (Wes = 0.25) the ˚ dependence of Sxx − Syy is reversed. This Weissenberg number marks the transition from random to oriented tube segments. At higher Weissenberg numbers tube segments in the pure melt become fully aligned, whereas in the suspension the velocity perturbations caused by the particles keep the average order parameter Sxx − Syy away from full alignment, even at large Weissenberg numbers. The effect of this change in Sxx − Syy is to reduce the ratio of to . The ratio decreases further as increases and then drops steeply to a minimum at the Weissenberg number at which reaches its maximum value of q in the pure material. As shown in Fig. 7b the main effect of the parameter q is to vary the position of this minimum.
Fig. 6. Comparison of the average backbone extension as a function of the stretch Weissenberg number Wes for ˚ = 0, 0.126 and 0.237 with q = 3. (a) For small Wes the average extension in the suspension is higher than in the pure melt. (b) However, at higher Wes the hydrodynamic disturbance created by the particles reduces the average extension.
288
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Fig. 7. Graphs of the ratio of the average backbone extension in suspension, , to for the unfilled melt: (a) q = 3 for ˚ = 0.126 and 0.237, (b) ˚ = 0.237 for q = 2 and 3.
Fig. 8. Comparison of the backbone tube orientation Sxx − Syy as a function of Weissenberg number We between ˚ = 0.126 and 0.237 and the pure melt, ˚ = 0.
These results show that while the viscoelastic stress increases with particle area fraction at very low Weissenberg numbers, at higher Weissenberg numbers the polymer stress contribution is reduced by adding particles. This reduction is due to a combination of two factors both associated with the velocity disturbance caused by the particles. First, the backbone tube segments do not become fully aligned, reducing the value of Sxx − Syy . Second, the polymers
Fig. 10. Ratio of the average of bulk viscosity against stretch Weissenberg numbers for different number of particles per unit for ˚ = 0.1257. The hydrodynamic interactions between particles create an artificial increase of the extensional viscosity.
do not become fully extended everywhere so that remains less than q. 5.2. Multiple particles per cell The previous calculations contain a single particle per unit cell, so that the interparticle spacing is set by the lattice. In these simu-
Fig. 9. The effect of increasing the stretch Weissenberg number to the behaviour of the viscosity (a). Transient viscosity against strain for Wes = 0.25. The viscosity of the filled system is higher than that of the unfilled system (b). Transient viscosity against strain for Wes = 3. In this case the viscosity of the filled system is lower than that of the unfilled system. At higher stretch Weissenberg number, the presence of the particles lowers the viscosity of the suspension.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
289
Fig. 11. The effect of increasing the area fraction of the particle, to the rheology of the suspension (a). The average bulk viscosity of the suspension against area fraction ˚ (b). The strain hardening parameter against area fraction.
lations the hydrodynamic interactions between particles produces a periodic disturbance in the extensional viscosity which is an artifact of the constraint imposed by the lattice structure. In order to remove the symmetry imposed by the lattice system we simulate the motion of multiple particles per cell, whose initial positions are chosen at random (subject to them not overlapping). For Newtonian fluids, we found [1] that while the single particle system gave the correct qualitative behaviour for the extensional viscosity, it underpredicted the nonlinear increase in extensional viscosity at higher area fractions. Thus we need to ensure that the reduction in strain hardening is not an artifact of the lattice structure, but also occurs when the particles are not constrained to a lattice. Multiparticle simulations also allow us to investigate whether particles form aggregates, or develop an anisotropic pair distribution function. In shear flows particles in viscoelastic fluids aggregate into long chains [23]. However, there have been no reports of similar structures in extensional flows. We begin by considering how the number of particles in the unit cell affects the calculated bulk viscosity. In Fig. 9 we compare the results for 1, 5 and 10 particles per cell at Wes = 0.25 and Wes = 3 for ˚ = 0.1257 and q = 10. At the lower stretch Weissenberg number the viscosity of the suspension is higher than that of the unfilled melt for all three simulations. However, the 5 and 10 particle simulations give a slightly higher prediction of the suspension viscosity than the single particle simulation. At the high Weissenberg number, shown in Fig. 9b, the suspension viscosity is lower than that of the pure melt in all simulations, but again the multiparticle simulations predict a higher extensional viscosity than the single particle simulations. In their simulations with the Giesekus model D’Avino et al. found that average extensional viscosity always increases with ˚, though the rate of increase decreases with increasing Weissenberg number. However, their results were limited to Weissenberg numbers below one and the Giesekus model does not have the same asymptotic behaviour as the pom-pom model at high Weissenberg numbers. In the multiparticle simulations the extensional viscosity is no longer periodic and contains fluctuations in the viscous stress contribution due to close interactions between particles. However we can still obtain average values by averaging over a long simulation. In Fig. 10 we compare the ratio of the average bulk viscosity of suspension to that of the pure fluid for the simulations with 1, 5 or 10 particles at an area fraction of ˚ = 0.1257. The ratio decreases with increasing Weissenberg number for all simulations. The only difference between the results for different numbers of particles is a vertical shift, which can be accounted for by the differences in the calculated value of the bulk extensional viscosity in a Newtonian fluid between single and multiparticle simulations [1]. This suggests that although the single particle simulations underpredict the
suspension viscosity, the prediction that adding fillers reduces the extensional viscosity of a suspension at high Weissenberg number remains valid. We focus now on the properties of the suspension at high Weissenberg number Wes > 1 and large arm number q (=10) where the unfilled melt is strongly strain hardening. The simulations contain 25 particle per cell and were performed with a mesh containing approximately 12,000 elements. The grid was adaptively refined in the regions between surrounding the particles and the gaps between close neighbouring particles were carefully resolved in order to prevent the particles from overlapping. This also required redu−1 cing the time step to a minimum of tn = 0.000125˙ when two particles become close to one another. In Fig. 11a we plot the average bulk viscosity of the suspension as a function of ˚ for Wes = 2 and 2.25. At this range of Weissenberg numbers the bulk extensional viscosity decreases with increasing particle area fraction. At the lower Weissenberg number Wes = 2 the decrease is small, with only 1% decrease in extensional viscosity for a suspension with ˚ = 0.16 compared to the unfilled polymer. However, increasing the Weissenberg number to 2.25 produces a larger decrease. Since adding particles increases the extensional viscosity of the suspension at small Weissenberg numbers, to obtain a quantitative measure of the degree of strain-rate thickening in a suspension, we define the strain-rate thickening parameter [14,20,6] as the ratio of the bulk extensional viscosity to its value in the limit of zero Weissenberg number, (Wes , ˚, t) =
(Wes , ˚, t) . (0, ˚, t)
(26)
At steady state, (0, ˚, t) is equal to the extensional viscosity when the suspending fluid is a Newtonian fluid of viscosity s (1 + (1/3)c) and so we can make use of the results calculated in reference [1]. The values of for Wes = 2 and 2.25 are plotted in Fig. 11b. The strain-rate thickening parameter decreases with increasing ˚ in qualitative agreement with both experimental measurements [16,20] and the simulations of D’Avino et al. [6]. In Fig. 12,1 we show snapshots of the backbone stretch distribution in the unit cell together with three of its images for the case of ˚ = 0.16 at a stretch Weissenberg number Wes = 2.25 and q = 10. In these pictures regions of higher backbone extension are denoted by darker shading. In Fig. 12a, taken after a strain of one, it can be seen that highest values of form thin lines connecting particles that have been compressed together in the y direction and then stretched apart in the x-direction. In their simulations D’Avino et
1
The motion is illustrated in the file named movieSuspension.mpg.
290
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Fig. 12. Snapshot showing the backbone stretch , for a simulation of 25 particles at ˚ = 0.16 displayed together with 3 of its images, under planar extensional flow at (a) = 1.0; (b) = 2.5; (c) = 2.75 and (d) = 3.0.
al. [6] showed that these regions are where the highest velocity gradients are found. Fig. 12b–d show the stretch distribution in the suspension for strains between 2.5 and 3 and show how the shape of the flow domain changes during a lattice period. At these strains the unfilled material has reached steady state, but the microstructure in the suspension continues to evolve due to the hydrodynamic disturbances caused by the particles. Although the distribution of is more uniform than found at earlier times in Fig. 12a, it is still possible to see the presence of multiple dark horizontal lines resulting from previous pair interactions. In addition we also observe regions of lower stretch at the sides of the particles where the particles screen the extensional flow allowing the stretch to relax. Finally, we examine whether viscoelasticity causes particles to aggregate in extensional flow. We define the local density of par-
Fig. 13. Plot of the local particle density, l , as a function of Hencky strain for one particular initial configuration for Wes = 2 and 3. The dashed line indicates the value l = 0.78 for non-overlapping randomly distributed particles.
ticles, l to be the average number of other particles whose centres lie within a region of 3 radii around a given particle, so that l is the integral of the pair distribution function over a circle whose radius is 3 times that of the particles. Fig. 13 shows the evolution of l for one particular initial configurations at stretch Weissenberg numbers of 2 and 3. The values of l fluctuate around l = 0.78 the mean value for randomly distributed non-overlapping particles. We have run a number of different initial configuration and although there are variations in l due to the finite size of the system we do not find any evidence of aggregation. Nor do we find any difference between the average distances between neighbouring particles in the compressional and extensional directions. This was also the case with our simulations of suspensions in Newtonian fluids [1]. 6. Conclusion In this paper we report on a numerical study of the effect of filler particles on the extensional rheology of a branched polymer melt. In Newtonian fluids the addition of filler particles always increases the suspension viscosity. However, we find that at high Weissenberg numbers fillers not only reduce the degree of strain hardening of the extensional viscosity (as measured by the ratio of the extensional viscosity to the linear envelope), but they can reduce the extensional viscosity even below that of the unfilled melt. This is a stronger result than reported by D’Avino et al. [6] who found a reduction in strain hardening, but the viscosity of the suspension still increased with area fraction. Although the fillers increase the average extension-rate in the suspension, which would be expected to increase the molecular extension, the localised velocity disturbances caused by the particles disrupt the tube orientation distribution causing the stretch to relax. The effect is strongest at high Weissenberg numbers for fluids that are strongly strain hardening and provides a qualitative explanation for the dramatic reduction in strain hardening observed in experiments. Although we have only reported results for the pom-pom model here, we have also studied other constitutive models including the FENE-CR model and find qualitatively similar behaviour. For the FENE-CR model we find that for an extensibility L = 10 the exten-
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
sional viscosity in a suspension is lower than that of the pure fluid for Weissenberg numbers greater than 2, but that for a lower extensibility, L = 5, the suspension viscosity is always higher than that of the pure fluid, though the degree of strain hardening is reduced. In the Giesekus model used by D’Avino et al. [6] the high Weissenberg extensional viscosity is controlled by the parameter ˛ which is roughly equivalent to 1/L2 in the FENE-CR model. Thus the value of ˛ = 0.1 used in the simulations corresponds to the low values of L where we find that the extensional viscosity remains higher than that of the pure fluid. The Lagrangian numerical method is well suited to this type of simulation. It allows the simulations to be performed on a boundary-free implementation of the self-replicating Kraynik and Reinelt lattice, so that simulations can be continued to high strains without the unit cell becoming highly distorted [1]. While the Lagrangian property of the mesh allows the constitutive equation to be solved as ordinary differential equations in the frame of the deforming fluid. Although we have only considered rigid circular particles in this paper it is possible to extend this to both non-spherical rigid particles and also deformable particles such as droplets or bubbles. The simulations in this paper are two-dimensional, whereas the real system is three-dimensional. The method can be extended to three-dimensional suspensions by replacing triangular elements with tetrahedral elements. However, the additional complexity of reconnecting a three-dimensional compared to a two-dimensional grid makes the Lagrangian scheme computationally expensive compared to other methods that produce less mesh distortion. Another restriction is that our method relies on the existence of a periodically self-replicating lattice and so cannot be used for uniaxial extension for which no such lattice exists. Acknowledgements We are grateful to the EPSRC for financial support through the Microscale Polymer Processing project GR/T11807/01. We would also like to thank Jonathan Shewchuk for the “Triangle” grid generation program and Tim Nicholson for the “FlowDis” visualisation software. References [1] M. Ahamadi, O.G. Harlen, A Lagrangian finite element method for simulation of a suspension under planar extensional flow, J. Comp. Phys. 227 (2008) 7543. [2] G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. 41 (1970) 545. [3] G.B. Bishko, O.G. Harlen, T.C.B. McLeish, T.M. Nicholson, Numerical simulation of the transient flow of branched polymer melts through a planar contraction using the pom-pom model, J. Non-Newtonian Fluid Mech. 82 (1999) 255.
291
[4] R. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag–strain coupling in branched polymer melts, J. Rheol. 44 (2000) 121. [5] G. D’Avino, P.L. Maffettone, M.A. Hulsen, G.M.W. Peters, A numerical method for simulating concentrated rigid particle suspensions in an elongational flow using a fixed grid, J. Comp. Phys. 226 (2007) 688. [6] G. D’Avino, P.L. Maffettone, M.A. Hulsen, G.M.W. Peters, Numerical simulation of planar elongational flow of concentrated rigid particle suspensions in a viscoelastic fluid, J. Non-Newtonian Fluid Mech. 150 (2008) 65. [7] B.N. Delaunay, Sur la sphere vide, Bull. Acad. Sci. USSR VII: Class. Sci. Math. (1934) 793. [8] R. Glowinski, T.-W. Pan, T.I. Hesla, D. Joseph, A distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiphase Flow 25 (1999) 755. [9] O.G. Harlen, J.M. Rallison, P. Szabo, A split Lagrangian–Eulerian method for simulating transient viscoelastic flows, J. Non-Newtonian Fluid Mech. 60 (1995) 81. [10] D.M. Heyes, Molecular dynamics simulations of extensional, sheet and unidirectional flow, Chem. Phys. 98 (1985) 15. [11] W.R. Hwang, M.A. Hulsen, Direct numerical simulations of hard particle suspensions in planar elongational flow, J. Non-Newtonian Fluid Mech. 136 (2006) 167. [12] W.R. Hwang, M.A. Hulsen, H.E.H. Meijer, Direct simulation of particle suspensions in sliding bi-periodic frames, J. Comput. Phys. 194 (2004) 742. [13] N.J. Inkson, T.C.B. McLeish, O.G. Harlen, D.J. Groves, Predicting low density polyethylene melt rheology in elongational and shear flows with pom-pom constitutive equations, J. Rheol. 43 (4) (1999) 873. [14] O. Ishizuka, K. Koyama, Elongational viscosity at a constant elongational strain rate of polypropylene melt, Polymer 21 (1980) 164. [15] M. Kobayashi, T. Takahashi, J. Takimoto, K. Koyama, Flow-induced whisker orientation and viscosity for molten composite systems in a uniaxial elongational flow-field, Polymer 36 (1995) 3927. [16] M. Kobayashi, T. Takahashi, J. Takimoto, K. Koyama, Influence of glass beads on the elongational viscosity of polyethylene with anomalous strain rate dependence of the strain-hardening, Polymer 37 (1996) 3745. [17] A.M. Kraynik, D.A. Reinelt, Extensional motions of spatially periodic lattices, Int. J. Multiphase Flow 18 (1992) 1045. [18] A.W. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C 5 (1921). [19] T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the pom-pom polymer, J. Rheol. 42 (1) (1998) 81. [20] J.F. Le Meins, P. Moldenaers, J. Mewis, Suspension of monodisperse spheres in polymer melts: particle size effects in extensional flow, Rheol. Acta 42 (2003) 184. [21] K. Ogura, M. Takahashi, Uniaxial extensional flow and thermoformability of PMMA melts with very high molecular weight component, Nihon Reoroji Gakkaishi [J. Soc. Rheol. Japan] 28 (2000) 99. [22] C.J.S. Petrie, One hundred years of extensional flow, J. Non-Newtonian Fluid Mech. 137 (2006) 1. [23] R. Scirocco, J. Vermant, J. Mewis, Effect of the viscoelasticity of the suspending fluid on structure formation in suspensions, J. Non-Newtonian Fluid Mech. 117 (2004) 183. [24] J.R. Shewchuk, Triangle: engineering a 2D quality mesh generator and delaunay triangulator, Lect. Notes Comp. Sci. 1148 (1996) 203. [25] R. Sibson, Locally equiangular triangulation, Comput. F 21 (1977) 243. [26] D. Silvester, A. Wathen, Fast iterative solution of stabilised Stokes systems. Part 2: Using general block preconditioners, SIAM J. Numer. Anal. 30 (1994) 630. [27] W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Differential constitutive equations for polymer melts: the extended Pom-Pom model, J. Rheol. 45 (4) (2001) 823.
Contents lists available at ScienceDirect
Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm
Numerical study of the rheology of rigid fillers suspended in long-chain branched polymer under planar extensional flow M. Ahamadi ∗ , O.G Harlen Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, United Kingdom
a r t i c l e
i n f o
Article history: Received 1 September 2008 Received in revised form 3 December 2009 Accepted 7 January 2010
Keywords: Finite elements (65M60) Fluid mechanics Suspensions rheology Planar extension Periodic boundary conditions Viscoelastic fluid
a b s t r a c t We report a detailed numerical study of the rheology of two-dimensional rigid fillers suspended in branched polymer melt under planar extensional flow. The polymer melt is modelled using the pom-pom constitutive equation. The numerical method uses a finite element solution of the flow in a unit cell within the self-replicating lattice for planar extensional flow identified by Kraynik and Reinelt [Int. J. Multiphase flow 18 (1992) 1045]. This is implemented using a quotient space representation that maps all points in space onto points within the unit cell. We show that the Kraynik and Reinelt cell allows simulations of suspensions under planar extensional flow to be conducted to large strains in a truly periodic cell. The influence of both the pom-pom parameters and the particle area fraction on the rheology of the suspension are investigated. We find a reduction in the degree of extension-rate thickening with the increase of both particles concentration and Weissenberg numbers in agreement with published experimental and numerical findings on other viscoelastic models. This reduction is found to be due to flow disturbance created by the particles which disrupts the alignment of backbone tube segments with the extensional axis. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Understanding the extensional rheology of filled polymer melts is important for the prediction of material performance in many commercial applications. One of the most important classes of such fillers are solid particles, typically glass or carbon particles, suspended in a polymeric matrix. The rheology of polymer blends under extensional flow has been studied by many investigators [22]. For long-chain branched polymers the effective viscosity of the polymer during the start-up (defined as the ratio of the stress difference to the imposed strain rate) grows above the level predicted by linear viscoelasticity. This phenomenon, known as strain hardening, is associated with the stretching of sections of the molecular chains, which occur at strain rates that exceed the effective Rouse or stretch relaxation time of the chain. In entangled long-chain polymer melts this relaxation time is controlled by the constraints on the motion of the branch points due to entanglements. Strain hardening is a desirable feature that improves melt processability in blow molding, foaming, and melt spinning [21]. Thus the effective extensional viscosity and the degree to which it shows strain hardening at high strain rates is an important rheological measurement. For filled polymers it has been observed experimentally [15,16,20] and numerically [6] that the degree of strain hardening
∗ Corresponding author. E-mail address: [email protected] (M. Ahamadi). 0377-0257/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2010.01.002
generally decreases with increasing particle concentration. A precise explanation of this phenomenon is still lacking, although a number of hypotheses have been suggested. Kobayashi et al. [16] argued that the presence of the particles perturbs the deformation of matrix polymer chains around them and thus causes the suppression of strain hardening. In their numerical study of concentrated particle suspensions, D’Avino et al. [6] suggested that the reduction in strain hardening with increasing Weissenberg number is caused by the shear flow between the particles. The effect of adding fillers on the extensional properties depends upon the type of polymer matrix used. Le Meins et al. [20] measured the elongational viscosity of monodisperse polystyrene spheres suspended in a viscoelastic polyisobutene matrix. They observed an increase in the extensional viscosity with increasing filler volume fraction and strain rate. At large strain rate, strain hardening occurs in the pure polymer but is suppressed by adding particles. Kobayashi et al. [16] studied the effect of glass beads on the elongation flow of high density polyethylene (HDPE). Their experimental results showed weaker strain hardening with increasing strain rate in both the pure polymer, as well for the suspension. The suppression of the strain hardening by the beads is stronger at high strain rate where the strain hardening of the pure polymer is weak. In this paper we present a detailed numerical study of the rheology of fillers under planar extensional flow. We focus our attention on the case of hard and non-Brownian circular particles suspended in long-chain branched polymer melts modelled using the pompom constitutive model. This class of polymer was chosen because
282
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
long-chain branched melts exhibit strain hardening in extensional flow, whereas linear polymers are generally strain softening in both shear and extensional flow. The pom-pom model developed by McLeish and Larson [19] contains the key features of branched polymers, i.e., it models chain segments whose reptation is inhibited by the presence of branch points. The original single mode model is able to correct qualitative behaviour of branched polymer in both shear and extension (and multimode versions are able to capture this rheology quantitatively [13,4,27]). In addition, the parameters in the models are directly linked to the microstructure of the polymer melt. The idealised pom-pom molecule has a single backbone with multiple branches emerging from each end. In extensional flow the backbone becomes stretched if the extensionrate is faster than the stretch relaxation time, s of the backbone. However, since the tension forces at the branch-point must balance there is an upper-bound to this extension, q, equal to the number of arms connected to the branch-point. This upper-bond on extension limits the degree of extension hardening. In planar extensional flow, as shown in Fig. 1, fluid particles separate exponentially in time in the direction of extension. Thus for a complex fluid any finite computational domain will usually become very highly distorted within the characteristic time for the fluid to reach steady state. Consequently, in order to be able to reach a steady state, it is desirable to use a self-replicating lattice that recovers its original geometry after a period of strain. This is not straightforward to simulate since the dynamical motion of extensional flow involves a change in the shape of the unit cell within the computational domain. Kraynik and Reinelt [17] found a set of self-replicating lattice for planar extension flow, which we will use to impose the periodic structure of a two-dimensional computational domain ˝(t). The numerical method is based upon the Lagrangian–Eulerian approach of Harlen et al. [9]. In this application this approach provides a natural way to implement the particle motion within the Kraynik and Reinelt lattice, as well as enabling the constitutive equation to be solved in the natural co-deforming frame of the fluid. 2. Computational domain We consider a two-dimensional unbounded domain, ˝∞ (t), containing non-Brownian hard circular particles suspended in viscoelastic fluid under a planar extensional flow. In order to make the computations tractable we approximate the fluid microstructure by a spatially periodic fluid microstructure. Therefore we solve the flow in an initially rectangular unit cell of size H × contai-
ning N particles, ˝(t) that when periodically replicated covers the infinite domain. Under the action of the velocity gradient, the unit cell within the lattice deforms and in order to be able to continue to large strains it is desirable to use a self-replicating lattice that recovers its original geometry after a period of strain. In shear flow, this can be achieved using the Lees–Edwards [12,18] lattice in which cells slide across one another in layers parallel to the flow direction. However, in a planar extensional flow (see Fig. 1), a rectangular lattice aligned along the principle flow axes will not periodically replicate its structure. A lattice of this type, suggested by Heyes [10] is used in the simulation method developed by Hwang and Hulsen [11] for planar elongational flow of suspensions in a Newtonian fluid. However, because the aspect ratio of the unit cell increases exponentially in time their simulations are limited to moderate strains. D’Avino et al. [5] have recently developed an alternative fixed grid method for simulating extensional flow. This does not use periodic boundary conditions but uses a three layer domain in which particles leaving the domain are reintroduced at random positions on the inlet boundaries. Kraynik and Reinelt [17] found a set of self-replicating lattice for planar extension flow by using the condition of reproducibility of a two-dimensional square lattice L(t) = n1 p(t) + n2 q(t), where p(t), q(t) are linearly independent basis vectors, n1 , n2 integers. They demonstrated that the square lattice with the minimum period for replication was one where the√sides are aligned at an angle of approximately = tan−1 (2/(1 + 5)) = 31.7◦ to the principal axes √ of the flow, see Fig. 1b, with the Hencky strain p = log((3 + 5)/2) = 0.9624. The important outcome of the Kraynik and Reinelt’s findings is that for this initial orientation the square lattice will reproduce itself after a strain of p in the sense that vertices of the original lattice move to other vertex points of this lattice. The period of reproducibility of the Kraynik and Reinelt cell depends only on the Hencky strain and so this lattice may be used for a variable strain-rate simulation provided that the principle axes of strain remain the same. In our previous paper [1], we detailed how to implement the Kraynik and Reinelt lattice structure for suspension in Newtonian fluid by using an appropriate image system that map positions in the infinite domain ˝∞ (t) to positions in the computational domain, ˝(t). In this paper we show how our method can be extended to viscoelastic fluids. 3. Governing equations The flow is assumed to be incompressible, isothermal and inertialess, so that the equations of conservation of mass and
Fig. 1. (a) Particle paths in planar extensional flow. (b) The unit cells in the critical lattice.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
momentum are given by,
∇ · =0
(1)
∇ ·u=0
where u is the fluid velocity and is the stress tensor. The constitutive model we will consider is the pom-pom model of McLeish and Larson [19]. This model is obtained by considering a melt of pom-pom molecules formed by connecting two q-armed star polymers with a backbone chain. The presence of the branch points at the ends of the backbone chain inhibits its motion along its tube, so that stretch and orientation relaxation times of the backbone, s and , respectively, are controlled by the relaxation of the star arms. We shall consider a single pom-pom mode describing the stretch and orientation of the backbone together with a Newtonian stress term for the fast relaxing modes associated with the arms so that the stress is of the form = −pI + 2s D + P ,
(2) P
where p is the pressure, 2s D is the Newtonian stress and is the backbone stress contribution. In this study, we will adopt the differential version of the original pom-pom model, but with the modification introduced by Blackwell et al. [4]. P = 3G2
A . Tr(A)
(3)
Here the auxillary tensor A satisfies an upper-convected Maxwell equation ∇ 1 A = − (A − I),
(4)
283
around the surface of the particle l as a Lagrange multiplier to force the fluid inside particle l to behave as a rigid solid. Here n is the outward unit normal to the boundary of particle Pl . The angular velocity ωl and the translational velocity Ul are found from the conditions of no net force and torque on each particle,
fl ds = 0,
l = 1, . . . , N
(8)
∂Pl
(x − Xl ) × fl ds = 0,
l = 1, . . . , N
where fl is the traction force density vector on particle l. The motions of particles are described by the kinematic equations dXl = Ul , dt
Xl |t=0 = Xl,0 ,
l = 1, . . . , N
1 ij = V
D e2(−1)/q = A : ∇u − ( − 1) Dt TrA s
where Vf is the volume of the fluid phase.
up to the maximum stretch, q. For this model we can define four dimensionless groups: the non-dimensional elastic modulus c = (3G/s ), the number of arms q, the Weissenberg number based on the orientation relaxation time We = ˙ and the Weissenberg number based upon the stretch relaxation time Wes = Wes /. Note that c is independent of ˙ and so remains fixed for a given material. In this paper we will show results for c = 9 and Wes = (1/2)We and vary q and We. The ratio We/Wes is proportional to the entanglement length of the backbone and a value of 2 would correspond to pom-pom molecules where the backbone length is equal to 5 entanglements and makes up 37% of the molecular weight of the polymer. Varying this ratio does not qualitatively affect the results. The extensional flow is imposed through the boundary conditions of the computational domain, which will be discussed in Section 4.2. 3.1. Particles The filler particles are assumed to be rigid and subject to zero net force and couple. Assuming that there is no slip between the filler particles and the fluid matrix, the fluid velocity at the surface of particle l is given by u = Ul + ωl × (x − Xl ),
l = 1, . . . , N
(6)
where Ul and ωl = ωl zˆ are the unknown velocity and angular velocity of particle l (with zˆ the unit vector in the zˆ direction). To enforce these boundary conditions we impose a surface force density fl (s) = n,
(7)
(10)
d l = ωl , l |t=0 = l,0 , l = 1, . . . , N (11) dt where l is the angular rotation of particle l. Note that the kinematic Eqs. (10) and (11) are decoupled from Eqs. (6), (8) and (9). Hence these equations of motions can be solved separately by using a standard ordinary differential equation method as we will describe later in Section 4.2. The quantity measured in rheological experiments is the spaced averaged stress of the suspension. We shall refer to this as the “bulk” stress to distinguish it from time-averaged quantities. Batchelor [2] derived an expression for the bulk stress in a suspension as sum of the stress in the fluid phase plus a contribution from the surface forces on the particles, in the form
where is the orientation relaxation time, while the backbone stretch, , is given by (5)
(9)
∂Pl
Vf
1 ij dVf + 2
[(xi − Xi )fj + fi (xj − Xj )ds]
(12)
∂Pl
4. Computational method 4.1. Finite element formulation We shall use the numerical method based on the split Eulerian–Lagrangian technique of Harlen et al. [9]. This uses a decoupled method in which the constitutive Eqs. (4) and (5) are solved as ordinary differential equations in the co-deforming (Lagrangian frame) within a separate step from the calculation of the fluid velocity. The fluid velocity and pressure are calculated for the current values of microstructural variables A and and positions of the particles. This is achieved by solving the Stokes Eq. (1) with Eqs. (8) and (9), and the no-slip boundary condition (6) via a finite element method. We then step forward in time to find the new values of A, and particle positions by solving the constitutive Eqs. (4) and (5). To render the Eqs. (1), (6), (8) and (9) in finite element form we shall adopt the combined equation of motion of Glowinski et al. [8] to derive the weak form of the equation for the solid–liquid mixture. The natural combined velocity space for the fluid and particle velocities is given by 2
V = {((v, Vl , l )|v ∈ H1 (˝)) , Vl ∈ R2 , l ∈ R, v = Vl + l zˆ × (x − Xl ) in Pl (t), and v biperiodic on ∂˝}.
(13)
In the distributed Lagrange multiplier method the extended weak form for whole domain can be obtained by removing the
284
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
constraint Eq. (6) from the velocity space and enforcing it weakly as a side constraint. This is done by introducing a Lagrange multiplier (fl (s) in our case), which can be interpreted as the traction force required to maintain the rigid-body motion of particle Pl (t). In our methodology, the constraint Eqs. (8) and (9) are used to determine the unknown Ul and ωl and therefore they must be incorporated into the final weak form.
where p(t), q(t) are the base vectors of the Kraynik and Reinelt cell and W = {(a, b) ∈ R × R such
Q : R2 →
−p∇ · w dV +
˝
∇ u : ∇ w dV − ˝
N
l=1
fl · w dS
(14)
˝
q∇ · u dV = 0,
−
(15)
˝
l · (u − (Ul + ωl × (x − Xl ))) dS = 0,
(16)
∂Pl (t)
fl dS = 0,
(17)
fl × (x − Xl ) dS = 0,
(18)
∂Pl (t)
∂Pl (t)
where l = 1, . . . , N and w, q, l are respectively the variations for the fluid velocity, pressure, and traction force. We use triangular P1 − P1 elements, which have linear continuous pressures and velocities, and a continuous linear interpolation for the force density fl around the particles to discretise the weak form [Eq. (14)–(18)]. The viscoelastic stress P is piecewise constant, P0 , over elements. The element combination P1 − P1 for the velocity and pressure is unstable to spurious pressure modes. These are removed by introducing a pressure stabilisation in which the weak form of continuity Eq. (15) is replaced by
q∇ · u dx − ˇ
− ˝
is the quotient map allowing us to identify all images. From the general form of the quotient map, Eq. (21), one can derive an explicit relation between a point in ˝(t) and its images in ˝∞ and between the velocities and other quantities held at these points. Thus after a strain of = ˙ t in planar extensional flow the images of a point (xlattice , ylattice ) are given by ximag = xlattice + me p01 + ne q01 H
−
(21)
ximag → Q(ximag ) = [xlattice ]
∂Pl (t)
P : ∇ w dV,
=−
R2 W
h ∇ q · ∇ p dx = 0, 2
(19)
˝
where h2 is twice the area of the triangle and ˇ is a positive constant. Following Silvester and Wathen [26], we choose ˇ = 0.025 as this gives the optimal convergence rate for the Stokes problem. For a given finite element mesh, the discretisation of Eqs. (14)–(18) leads to a symmetric indefinite linear system of algebraic equations. These are solved using a preconditioned conjugate residuals method with block preconditioner of the form suggested by Silvester and Wathen [26]. 4.2. Implementation of the Kraynik and Reinelt lattice The extensional flow is imposed through the periodically selfreplicating Kraynik and Reinelt lattice using the method described in ref. [1]. We generate an initial biperiodic mesh for the computational domain, ˝(t0 ), in which the nodes on opposite edges are identified as being equivalent points. This is done using the “Triangle” grid generation program [24] to generate an initial grid in the unit cell and then relabelling boundary points so that equivalent points have the same label. This equivalence can be extended further by noting that every point in the computational domain has an equivalence class of image nodes in the infinite domain, ˝∞ . Let ∼ be an equivalence relation on V = R2 defined by x∼r ⇔ ∃ m ∈ Z, n ∈ Z such
(20)
yimag = ylattice + me− p02 + ne− q02 H,
(22)
where the integer pair (m, n) denote the shift in the periodic image. Since the finite element matrices depend only upon the relative displacement of nodes and not their absolute position we define the displacement dxlattice of nodes on an element to be the displacement between images that lies in the space dxlattice = dpp(t) + dqq(t) −
< dp ≤ , 2 2
−
with
H H < dq ≤ . 2 2
(23)
The external extensional flow is imposed by the condition that the fluid velocity uimag of the (m, n) image is given by the time derivative of the quotient map given in Eq. (21). Consequently, the fluid velocity, uimag is given by uimag = ulattice + ˙ (Tu , Tv ), me˙ t p01 + ne˙ t q01 H
(24) −me−˙ t p02 − ne−˙ t q02 H.
and Tv = where Tu = All other variables, including pressure and the microstructure variables A and have the same value at each image point. A major drawback of using a Lagrangian mesh is that severe deformation-induced mesh distortions develop due to velocity gradients within the fluid, which will degrade the accuracy of the finite element solution. A partial solution to this problem is to retain the nodes as material points, but reconnect them in the way that produces a Delaunay triangulation [7]. The Delaunay triangulation ensures global equiangular triangulation in the sense that it maximises the minimum angle in any triangle [25]. The Lagrangian mesh movement combined with Delaunay reconnection and the use of the equivalence class to relate points with their images automatically reproduces the self-replicating lattice of Kraynik and Reinelt. This is illustrated in Fig. 2 where the evolution of a very simple mesh consisting of eight triangular elements connecting four replicated points and their appropriate images under planar extension is shown. The initial mesh is shown in Fig. 2a with the square unit cell (indicated by the bold lines) oriented at an angle = 31.7◦ to the x-axis. Fig. 2b and c show the lattice after a strain of (1/2)p . The original unit cell (shown with solid bold lines in Fig. 2b) is deformed into a parallelogram, however, we can reform an equivalent square unit cell (shown with dashed lines), by replacing p with p + q. Between Fig. 2b and c we have performed a Delaunay reconnection. Although we have recovered a square lattice it is not the original lattice as it is now oriented at 58.3◦ to the x-axis. Fig. 2d shows the lattice after a strain of p = 0.962424. We can now recover the original square unit cell by choosing appropriate images. A further Delaunay reconnection will reproduce the original grid shown in Fig. 2d.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
285
Fig. 2. A simple example showing how the mesh movement and reconnection replicate the unit cell after a strain of p = 0.9624. (a) The initial lattice with a unit cell made up of eight triangular elements. (b) The lattice after a strain of (1/2)p . The initial unit cell is shown with the solid line, while the dashed line shows an equivalent square unit cell formed by replacing two of the elements with their images. (c) The lattice after a strain of (1/2)p following a Delaunay reconnection. (d) The lattice after a strain of p = 0.9624. The original unit cell and lattice can be recovered by a further replacement of two elements with their images and a further Delaunay reconnection.
4.3. Lagrangian calculation of the viscoelastic stress
we employ the methods of grid improvement detailed in [1] during the calculation.
The value of the viscoelastic stress within each element is calculated from the values of microstructure variables A and , which are updated by integrating Eqs. (4) and (5) forward in time in a frame deforming with the fluid [3]. This is done by moving the vertices of the triangles with their current velocity. Each triangle is, therefore, a material volume (since the velocity varies linearly over the triangle) and the vectors along any two sides of the triangle form a set of base vectors for the frame that moves and deforms with the fluid. The time derivative in this frame is therefore the upper-convected derivative, so that A may be found by solving the ordinary differential equation ˜ dA 1 ˜ ˜ A−I =− dt
(25)
˜ and ˜I are the components of the A and the identity tensor, I where A in triangle coordinates [9]. Similarly, the Lagrangian time derivative (D/Dt) becomes (d/dt) in this frame. The resulting ordinary differential equations are solved by forward time-stepping and then A is transformed back into Cartesian coordinates. The grid is then reconnected where necessary to give a Delaunay triangulation. For the simulations containing particles the Delaunay reconnection is not allowed to reconnect an edge on a particle boundary. However edges on the boundary of the unit cell are reconnected so that the outline of the original cell is lost. As noted earlier, due to the equivalence class it does not matter which image of a triangle is used in the assembly. Once particles are introduced reconnection alone is not sufficient to maintain mesh quality and so
5. Results and discussion In planar extensional flow the bulk extensional viscosity is defined as =
− xx yy ˙
.
For a purely linear viscoelastic fluid this is equal to 4 times the transient shear viscosity. However, in a strain-hardening material it will be a function of both t and the applied strain rate ˙ . For the pom-pom model given in Eq. (3) the extensional viscosity is given by = 4s +
3G 2 c 2 (Sxx − Syy ) (Sxx − Syy ) = s 4 + We ˙
where S = A/Tr(A) is the backbone tube orientation distribution. For large strain rates (Wes > 1) both and Sxx − Syy approach their limiting values of q and unity, respectively, so that /s will grow from 4 to a limiting value of 4 + (c/We)q2 at high strain. This is shown in Fig. 3a for the case of q = 2 and Wes = (1/2)We = 1.75 where it can be seen that steady state is achieved after a Hencky strain of 3. At intermediate strain rates (We > 0.5, Wes < 1) the backbone tubes become fully aligned, but do not stretch so that remains close to unity, while for strain rates for which We 1 we recover the linear viscoelastic limit.
286
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Fig. 3. Transient extensional viscosity for a pom-pom fluid with q = 2 and the stretch Weissenberg number Wes = 1.75: (a) pure fluid, (b) single particle of radius 0.3 (˚ = 0.2827) for the three different meshes given in Table 1.
Table 1 Details of the four finite element meshes used to calculate the extensional viscosity at t = 0 for a circular particle of area fraction ˚ = 0.2827. Mesh
Nb
Np
Area min
Elements, Ne
Average of bulk viscosity
Mesh 1 Mesh 2 Mesh 3
32 44 62
140 196 277
0.0025 0.001 5 × 10−4
1150 2210 4054
14.2439592866843 14.2284206433883 14.2163790630068
The addition of filler particles affects the flow in number of ways. First, since the particles do not deform the average strain rate in the −1 fluid phase increases by a factor (1 − ˚) , where ˚ is the area fraction of the particles and so increases the effective values of the Weissenberg numbers. However, the particles also perturb the flow away from pure planar extension which will rotate the tube segments away from perfect alignment. 5.1. Single particle per cell We first consider the case of a single particle at the centre of computational cell. In Fig. 3b we plot the bulk extensional viscosity for a particle of radius 0.3 (giving a particle area fraction ˚ = 0.2827) suspended in a pom-pom fluid with q = 2 and Wes = 1.75. In order to test the spatial convergence of our numerical method we compare the results from three different uniform meshes whose details are given in Table 1. Here, Nb is the number of nodes along each edge of the unit cell, Np the number of nodes around the particle, and Area min is the threshold value for the area of each triangle
used by “Triangle”. Mesh 2 is shown in Fig. 4. In our previous paper [1] we found that the spatial convergence for Newtonian fluids was of order h2 , and it can be seen that there is very little difference between the results for the two finest meshes. In subsequent calculations we will use meshes similar to mesh 3. Comparing Fig. 3 a and b we see that the bulk viscosity of the filled systems grows in a similar manner to that of the pure fluid, but with a small oscillation superimposed on it. The frequency of this oscillation is equal to twice that of the self-replicating lattice and is due to hydrodynamic interactions of the particle with its images. The maximum extensional viscosity occurs when the particles are arranged in a square lattice, which occurs twice in each period p . The lattice configuration after a strain of (1/2)p is square, but is rotated relative to its initial configuration. By averaging over the period of this oscillation at large strains we can obtain the average bulk viscosity and compare this with the steady state value for the pure fluid. In Fig. 5a we compare the average bulk viscosity for ˚ = 0.237 with the steady viscosity of the pure fluid at a stretch Weissenberg number Wes = 2.25 for varying numbers of arms q. As expected, increasing the number of arms increases the extensional viscosity for both the filled and unfilled materials. For Newtonian fluids the addition of rigid filler particles always increases the extensional viscosity, however, at this Weissenberg number adding particles reduces the bulk extensional viscosity. As noted earlier, since Wes is greater than unity, the backbone segments become fully aligned and stretched in the pure fluid so that they contribute their maximum possible value to the extensional stress. Adding particles disrupts the alignment and the consequent reduction in the poly-
Fig. 4. The computational domain for the single particle problem. (a) The unit cell at time t = 0 in relation to its periodic images. (b) The finite element Mesh 2 used in the calculations described in Table 1.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
287
Fig. 5. The effect of the varying the arm number q on the extensional viscosity of rheology of a suspension at ˚ = 0.237. (a) Plot of the average bulk extensional viscosity against arm number for Wes = 2.25. (b) Plot of the ratio of the suspension average viscosity to that of the pure fluid for Wes = 1, 1.25, 1.75 and 2.25.
meric stress is sufficient to overcome the addition to the solvent stress resulting from the addition of particles. By plotting the ratio of average bulk extensional viscosity of the filled system to the steady extensional viscosity of the suspending fluid we can see more easily the effect that adding particles has on the extensional viscosity. In Fig. 5b we plot this ratio for stretch Weissenberg numbers ranging from 1 to 2.25. The ratio decreases with increasing stretch Weissenberg number showing that the reduction in extensional viscosity occurs at high strain rates where the viscoelastic stress contribution is at its maximum in the unfilled system. For the two highest stretch Weissenberg numbers the ratio decreases with increasing arm number, however at a stretch Weissenberg number of unity the ratio increases with increasing arm number. This Weissenberg number corresponds to the transition between stretched and unstretched backbones, and so the backbone segments are not fully extended in the pure material. This means that in regions of the flow where the local extension rates are higher the backbone segments can become more highly extended than in the pure material. In Fig. 6 we compare the average backbone extension in suspensions of area fractions ˚ = 0.126 and 0.237 to the value of in the pure melt for a 3-armed pom-pom. At stretch Weissenberg numbers below 0.3 (Fig. 6a) the average extension increases with filler concentration. Because the particles do not strain, the average extension-rate in melt phase of a suspension is higher than in −1 the pure material by a factor (1 − ˚) . Consequently, the effective values of the Weissenberg number in the melt phase of a suspension is higher than the value calculated from the bulk deformation rate. Infact for small Weissenberg numbers the results can
be superimposed by shifting the Weissenberg number by this factor. However, at stretch Weissenberg number above 0.3 the effect is reversed (Fig. 6b) with the average backbone extension decreasing with filler concentration. In all three materials the degree of extension increases with We, but in the pure melt, achieves its maximum value for Wes > 2.5, whereas the velocity perturbations in the suspensions prevent every from reaching this value. By examining the ratio of the average backbone extension to the value of in the pure melt in Fig. 7a we can identify a number of different regimes. For very small Weissenberg numbers where is close to unity and the orientation is close to being isotropic, grows faster than due to the increase in the average extension-rate, so that the ratio is greater than unity. As can be seen in Fig. 8 at these low Weissenberg numbers the degree of tube alignment increases with ˚. However at around We = 0.5 (Wes = 0.25) the ˚ dependence of Sxx − Syy is reversed. This Weissenberg number marks the transition from random to oriented tube segments. At higher Weissenberg numbers tube segments in the pure melt become fully aligned, whereas in the suspension the velocity perturbations caused by the particles keep the average order parameter Sxx − Syy away from full alignment, even at large Weissenberg numbers. The effect of this change in Sxx − Syy is to reduce the ratio of to . The ratio decreases further as increases and then drops steeply to a minimum at the Weissenberg number at which reaches its maximum value of q in the pure material. As shown in Fig. 7b the main effect of the parameter q is to vary the position of this minimum.
Fig. 6. Comparison of the average backbone extension as a function of the stretch Weissenberg number Wes for ˚ = 0, 0.126 and 0.237 with q = 3. (a) For small Wes the average extension in the suspension is higher than in the pure melt. (b) However, at higher Wes the hydrodynamic disturbance created by the particles reduces the average extension.
288
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Fig. 7. Graphs of the ratio of the average backbone extension in suspension, , to for the unfilled melt: (a) q = 3 for ˚ = 0.126 and 0.237, (b) ˚ = 0.237 for q = 2 and 3.
Fig. 8. Comparison of the backbone tube orientation Sxx − Syy as a function of Weissenberg number We between ˚ = 0.126 and 0.237 and the pure melt, ˚ = 0.
These results show that while the viscoelastic stress increases with particle area fraction at very low Weissenberg numbers, at higher Weissenberg numbers the polymer stress contribution is reduced by adding particles. This reduction is due to a combination of two factors both associated with the velocity disturbance caused by the particles. First, the backbone tube segments do not become fully aligned, reducing the value of Sxx − Syy . Second, the polymers
Fig. 10. Ratio of the average of bulk viscosity against stretch Weissenberg numbers for different number of particles per unit for ˚ = 0.1257. The hydrodynamic interactions between particles create an artificial increase of the extensional viscosity.
do not become fully extended everywhere so that remains less than q. 5.2. Multiple particles per cell The previous calculations contain a single particle per unit cell, so that the interparticle spacing is set by the lattice. In these simu-
Fig. 9. The effect of increasing the stretch Weissenberg number to the behaviour of the viscosity (a). Transient viscosity against strain for Wes = 0.25. The viscosity of the filled system is higher than that of the unfilled system (b). Transient viscosity against strain for Wes = 3. In this case the viscosity of the filled system is lower than that of the unfilled system. At higher stretch Weissenberg number, the presence of the particles lowers the viscosity of the suspension.
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
289
Fig. 11. The effect of increasing the area fraction of the particle, to the rheology of the suspension (a). The average bulk viscosity of the suspension against area fraction ˚ (b). The strain hardening parameter against area fraction.
lations the hydrodynamic interactions between particles produces a periodic disturbance in the extensional viscosity which is an artifact of the constraint imposed by the lattice structure. In order to remove the symmetry imposed by the lattice system we simulate the motion of multiple particles per cell, whose initial positions are chosen at random (subject to them not overlapping). For Newtonian fluids, we found [1] that while the single particle system gave the correct qualitative behaviour for the extensional viscosity, it underpredicted the nonlinear increase in extensional viscosity at higher area fractions. Thus we need to ensure that the reduction in strain hardening is not an artifact of the lattice structure, but also occurs when the particles are not constrained to a lattice. Multiparticle simulations also allow us to investigate whether particles form aggregates, or develop an anisotropic pair distribution function. In shear flows particles in viscoelastic fluids aggregate into long chains [23]. However, there have been no reports of similar structures in extensional flows. We begin by considering how the number of particles in the unit cell affects the calculated bulk viscosity. In Fig. 9 we compare the results for 1, 5 and 10 particles per cell at Wes = 0.25 and Wes = 3 for ˚ = 0.1257 and q = 10. At the lower stretch Weissenberg number the viscosity of the suspension is higher than that of the unfilled melt for all three simulations. However, the 5 and 10 particle simulations give a slightly higher prediction of the suspension viscosity than the single particle simulation. At the high Weissenberg number, shown in Fig. 9b, the suspension viscosity is lower than that of the pure melt in all simulations, but again the multiparticle simulations predict a higher extensional viscosity than the single particle simulations. In their simulations with the Giesekus model D’Avino et al. found that average extensional viscosity always increases with ˚, though the rate of increase decreases with increasing Weissenberg number. However, their results were limited to Weissenberg numbers below one and the Giesekus model does not have the same asymptotic behaviour as the pom-pom model at high Weissenberg numbers. In the multiparticle simulations the extensional viscosity is no longer periodic and contains fluctuations in the viscous stress contribution due to close interactions between particles. However we can still obtain average values by averaging over a long simulation. In Fig. 10 we compare the ratio of the average bulk viscosity of suspension to that of the pure fluid for the simulations with 1, 5 or 10 particles at an area fraction of ˚ = 0.1257. The ratio decreases with increasing Weissenberg number for all simulations. The only difference between the results for different numbers of particles is a vertical shift, which can be accounted for by the differences in the calculated value of the bulk extensional viscosity in a Newtonian fluid between single and multiparticle simulations [1]. This suggests that although the single particle simulations underpredict the
suspension viscosity, the prediction that adding fillers reduces the extensional viscosity of a suspension at high Weissenberg number remains valid. We focus now on the properties of the suspension at high Weissenberg number Wes > 1 and large arm number q (=10) where the unfilled melt is strongly strain hardening. The simulations contain 25 particle per cell and were performed with a mesh containing approximately 12,000 elements. The grid was adaptively refined in the regions between surrounding the particles and the gaps between close neighbouring particles were carefully resolved in order to prevent the particles from overlapping. This also required redu−1 cing the time step to a minimum of tn = 0.000125˙ when two particles become close to one another. In Fig. 11a we plot the average bulk viscosity of the suspension as a function of ˚ for Wes = 2 and 2.25. At this range of Weissenberg numbers the bulk extensional viscosity decreases with increasing particle area fraction. At the lower Weissenberg number Wes = 2 the decrease is small, with only 1% decrease in extensional viscosity for a suspension with ˚ = 0.16 compared to the unfilled polymer. However, increasing the Weissenberg number to 2.25 produces a larger decrease. Since adding particles increases the extensional viscosity of the suspension at small Weissenberg numbers, to obtain a quantitative measure of the degree of strain-rate thickening in a suspension, we define the strain-rate thickening parameter [14,20,6] as the ratio of the bulk extensional viscosity to its value in the limit of zero Weissenberg number, (Wes , ˚, t) =
(Wes , ˚, t) . (0, ˚, t)
(26)
At steady state, (0, ˚, t) is equal to the extensional viscosity when the suspending fluid is a Newtonian fluid of viscosity s (1 + (1/3)c) and so we can make use of the results calculated in reference [1]. The values of for Wes = 2 and 2.25 are plotted in Fig. 11b. The strain-rate thickening parameter decreases with increasing ˚ in qualitative agreement with both experimental measurements [16,20] and the simulations of D’Avino et al. [6]. In Fig. 12,1 we show snapshots of the backbone stretch distribution in the unit cell together with three of its images for the case of ˚ = 0.16 at a stretch Weissenberg number Wes = 2.25 and q = 10. In these pictures regions of higher backbone extension are denoted by darker shading. In Fig. 12a, taken after a strain of one, it can be seen that highest values of form thin lines connecting particles that have been compressed together in the y direction and then stretched apart in the x-direction. In their simulations D’Avino et
1
The motion is illustrated in the file named movieSuspension.mpg.
290
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
Fig. 12. Snapshot showing the backbone stretch , for a simulation of 25 particles at ˚ = 0.16 displayed together with 3 of its images, under planar extensional flow at (a) = 1.0; (b) = 2.5; (c) = 2.75 and (d) = 3.0.
al. [6] showed that these regions are where the highest velocity gradients are found. Fig. 12b–d show the stretch distribution in the suspension for strains between 2.5 and 3 and show how the shape of the flow domain changes during a lattice period. At these strains the unfilled material has reached steady state, but the microstructure in the suspension continues to evolve due to the hydrodynamic disturbances caused by the particles. Although the distribution of is more uniform than found at earlier times in Fig. 12a, it is still possible to see the presence of multiple dark horizontal lines resulting from previous pair interactions. In addition we also observe regions of lower stretch at the sides of the particles where the particles screen the extensional flow allowing the stretch to relax. Finally, we examine whether viscoelasticity causes particles to aggregate in extensional flow. We define the local density of par-
Fig. 13. Plot of the local particle density, l , as a function of Hencky strain for one particular initial configuration for Wes = 2 and 3. The dashed line indicates the value l = 0.78 for non-overlapping randomly distributed particles.
ticles, l to be the average number of other particles whose centres lie within a region of 3 radii around a given particle, so that l is the integral of the pair distribution function over a circle whose radius is 3 times that of the particles. Fig. 13 shows the evolution of l for one particular initial configurations at stretch Weissenberg numbers of 2 and 3. The values of l fluctuate around l = 0.78 the mean value for randomly distributed non-overlapping particles. We have run a number of different initial configuration and although there are variations in l due to the finite size of the system we do not find any evidence of aggregation. Nor do we find any difference between the average distances between neighbouring particles in the compressional and extensional directions. This was also the case with our simulations of suspensions in Newtonian fluids [1]. 6. Conclusion In this paper we report on a numerical study of the effect of filler particles on the extensional rheology of a branched polymer melt. In Newtonian fluids the addition of filler particles always increases the suspension viscosity. However, we find that at high Weissenberg numbers fillers not only reduce the degree of strain hardening of the extensional viscosity (as measured by the ratio of the extensional viscosity to the linear envelope), but they can reduce the extensional viscosity even below that of the unfilled melt. This is a stronger result than reported by D’Avino et al. [6] who found a reduction in strain hardening, but the viscosity of the suspension still increased with area fraction. Although the fillers increase the average extension-rate in the suspension, which would be expected to increase the molecular extension, the localised velocity disturbances caused by the particles disrupt the tube orientation distribution causing the stretch to relax. The effect is strongest at high Weissenberg numbers for fluids that are strongly strain hardening and provides a qualitative explanation for the dramatic reduction in strain hardening observed in experiments. Although we have only reported results for the pom-pom model here, we have also studied other constitutive models including the FENE-CR model and find qualitatively similar behaviour. For the FENE-CR model we find that for an extensibility L = 10 the exten-
M. Ahamadi, O.G Harlen / J. Non-Newtonian Fluid Mech. 165 (2010) 281–291
sional viscosity in a suspension is lower than that of the pure fluid for Weissenberg numbers greater than 2, but that for a lower extensibility, L = 5, the suspension viscosity is always higher than that of the pure fluid, though the degree of strain hardening is reduced. In the Giesekus model used by D’Avino et al. [6] the high Weissenberg extensional viscosity is controlled by the parameter ˛ which is roughly equivalent to 1/L2 in the FENE-CR model. Thus the value of ˛ = 0.1 used in the simulations corresponds to the low values of L where we find that the extensional viscosity remains higher than that of the pure fluid. The Lagrangian numerical method is well suited to this type of simulation. It allows the simulations to be performed on a boundary-free implementation of the self-replicating Kraynik and Reinelt lattice, so that simulations can be continued to high strains without the unit cell becoming highly distorted [1]. While the Lagrangian property of the mesh allows the constitutive equation to be solved as ordinary differential equations in the frame of the deforming fluid. Although we have only considered rigid circular particles in this paper it is possible to extend this to both non-spherical rigid particles and also deformable particles such as droplets or bubbles. The simulations in this paper are two-dimensional, whereas the real system is three-dimensional. The method can be extended to three-dimensional suspensions by replacing triangular elements with tetrahedral elements. However, the additional complexity of reconnecting a three-dimensional compared to a two-dimensional grid makes the Lagrangian scheme computationally expensive compared to other methods that produce less mesh distortion. Another restriction is that our method relies on the existence of a periodically self-replicating lattice and so cannot be used for uniaxial extension for which no such lattice exists. Acknowledgements We are grateful to the EPSRC for financial support through the Microscale Polymer Processing project GR/T11807/01. We would also like to thank Jonathan Shewchuk for the “Triangle” grid generation program and Tim Nicholson for the “FlowDis” visualisation software. References [1] M. Ahamadi, O.G. Harlen, A Lagrangian finite element method for simulation of a suspension under planar extensional flow, J. Comp. Phys. 227 (2008) 7543. [2] G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. 41 (1970) 545. [3] G.B. Bishko, O.G. Harlen, T.C.B. McLeish, T.M. Nicholson, Numerical simulation of the transient flow of branched polymer melts through a planar contraction using the pom-pom model, J. Non-Newtonian Fluid Mech. 82 (1999) 255.
291
[4] R. Blackwell, T.C.B. McLeish, O.G. Harlen, Molecular drag–strain coupling in branched polymer melts, J. Rheol. 44 (2000) 121. [5] G. D’Avino, P.L. Maffettone, M.A. Hulsen, G.M.W. Peters, A numerical method for simulating concentrated rigid particle suspensions in an elongational flow using a fixed grid, J. Comp. Phys. 226 (2007) 688. [6] G. D’Avino, P.L. Maffettone, M.A. Hulsen, G.M.W. Peters, Numerical simulation of planar elongational flow of concentrated rigid particle suspensions in a viscoelastic fluid, J. Non-Newtonian Fluid Mech. 150 (2008) 65. [7] B.N. Delaunay, Sur la sphere vide, Bull. Acad. Sci. USSR VII: Class. Sci. Math. (1934) 793. [8] R. Glowinski, T.-W. Pan, T.I. Hesla, D. Joseph, A distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiphase Flow 25 (1999) 755. [9] O.G. Harlen, J.M. Rallison, P. Szabo, A split Lagrangian–Eulerian method for simulating transient viscoelastic flows, J. Non-Newtonian Fluid Mech. 60 (1995) 81. [10] D.M. Heyes, Molecular dynamics simulations of extensional, sheet and unidirectional flow, Chem. Phys. 98 (1985) 15. [11] W.R. Hwang, M.A. Hulsen, Direct numerical simulations of hard particle suspensions in planar elongational flow, J. Non-Newtonian Fluid Mech. 136 (2006) 167. [12] W.R. Hwang, M.A. Hulsen, H.E.H. Meijer, Direct simulation of particle suspensions in sliding bi-periodic frames, J. Comput. Phys. 194 (2004) 742. [13] N.J. Inkson, T.C.B. McLeish, O.G. Harlen, D.J. Groves, Predicting low density polyethylene melt rheology in elongational and shear flows with pom-pom constitutive equations, J. Rheol. 43 (4) (1999) 873. [14] O. Ishizuka, K. Koyama, Elongational viscosity at a constant elongational strain rate of polypropylene melt, Polymer 21 (1980) 164. [15] M. Kobayashi, T. Takahashi, J. Takimoto, K. Koyama, Flow-induced whisker orientation and viscosity for molten composite systems in a uniaxial elongational flow-field, Polymer 36 (1995) 3927. [16] M. Kobayashi, T. Takahashi, J. Takimoto, K. Koyama, Influence of glass beads on the elongational viscosity of polyethylene with anomalous strain rate dependence of the strain-hardening, Polymer 37 (1996) 3745. [17] A.M. Kraynik, D.A. Reinelt, Extensional motions of spatially periodic lattices, Int. J. Multiphase Flow 18 (1992) 1045. [18] A.W. Lees, S.F. Edwards, The computer study of transport processes under extreme conditions, J. Phys. C 5 (1921). [19] T.C.B. McLeish, R.G. Larson, Molecular constitutive equations for a class of branched polymers: the pom-pom polymer, J. Rheol. 42 (1) (1998) 81. [20] J.F. Le Meins, P. Moldenaers, J. Mewis, Suspension of monodisperse spheres in polymer melts: particle size effects in extensional flow, Rheol. Acta 42 (2003) 184. [21] K. Ogura, M. Takahashi, Uniaxial extensional flow and thermoformability of PMMA melts with very high molecular weight component, Nihon Reoroji Gakkaishi [J. Soc. Rheol. Japan] 28 (2000) 99. [22] C.J.S. Petrie, One hundred years of extensional flow, J. Non-Newtonian Fluid Mech. 137 (2006) 1. [23] R. Scirocco, J. Vermant, J. Mewis, Effect of the viscoelasticity of the suspending fluid on structure formation in suspensions, J. Non-Newtonian Fluid Mech. 117 (2004) 183. [24] J.R. Shewchuk, Triangle: engineering a 2D quality mesh generator and delaunay triangulator, Lect. Notes Comp. Sci. 1148 (1996) 203. [25] R. Sibson, Locally equiangular triangulation, Comput. F 21 (1977) 243. [26] D. Silvester, A. Wathen, Fast iterative solution of stabilised Stokes systems. Part 2: Using general block preconditioners, SIAM J. Numer. Anal. 30 (1994) 630. [27] W.M.H. Verbeeten, G.W.M. Peters, F.P.T. Baaijens, Differential constitutive equations for polymer melts: the extended Pom-Pom model, J. Rheol. 45 (4) (2001) 823.