Viscosity of curved fibers in suspension

J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Viscosity of curved fibers in suspension C.G. Joung∗ , N. Phan-Thien, X.J. Fan Department of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney, NSW 2006, Australia Received 28 November 2000; received in revised form 12 August 2001

Abstract In this paper, the relationship between fiber shape and relative viscosity of a fiber suspension is explored. A numerical simulation has been used to model non-Brownian curved rigid fibers in suspension under the influence of Newtonian shear flow. Curvature in the simulated fibers was taken to represent general deformities of real fibers in suspension. The simulation method was previously used by Joung et al. [J. Non-Newtonian Fluid Mech. 99 (2001) 1] to determine suspension viscosity for flexible fibers in suspension. When compared to the equivalent straight rigid fiber suspension, fiber curvature was found to contribute to a large increase in suspension viscosity. For typical semi-concentrated to concentrated suspensions, curved fibers were observed to produce viscosity increases of the order twice that of straight fiber suspensions. Results indicate that even a small bend in the fibers may cause a large bulk viscosity increase. Suspension viscosity is therefore highly dependent on the quality control measures taken during sample preparation. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Direct simulation; Flexible fibers; Suspension viscosity; Bent fibers; Curved fibers

1. Introduction The inclusion of short fibers in the melt during the injection molding process can significantly improve the mechanical performance of the finished molded product [2]. Therefore, this type of injection molding is a popular method for the enhancement of molded plastics. There are many contributing factors to the ultimate mechanical performance of the product. In particular, understanding the microstructure during the filling stage is advised if one is to take full advantage of the benefits of fiber injection molding. It is known that the flow field within the mould during the injection and cooling stages directly determines the final position and orientation of the embedded fibers, and hence the mechanical properties of the finished composite. If, for example, flow patterns cause the fibers to become aligned in one common orientation, the molded part is typically stiffer in that direction, and more compliant in transverse directions [3]. The embedded fiber position and orientation ∗

Corresponding author. E-mail address: [email protected] (C.G. Joung). 0377-0257/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 0 1 ) 0 0 1 6 3 - X

2

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Fig. 1. Photographs of rayon fibers used by Milliken et al. (Fig. 2 of [31]; reproduced with permission.)

may also drastically affect the warpage behavior of the molded part as it cools [4,5]. The ability to predict and control the fiber motion and orientation within an injection mould is therefore very important to the quality of injection molded products. The demand for this ability warrants the current efforts of researchers in the field of fiber suspension rheology and it is from this basis that the current work is undertaken. In this work, as in [1], the main focus is on the viscosity of fiber suspensions since viscosity is important to the flow field development of molten plastics. As in [1] we use a direct numerical simulation of fibers in suspension to observe the microstructural development in a shear flow and calculate suspension viscosity. In [1], however, the simulation fibers were flexible and deforming with time and flow conditions. In this work, we now force the fibers to be rigid but curved. The enforced rigidity eliminates all fiber flexibility effects which was the particular focus of [1]. Curvature in the fiber is now introduced as a metaphor for general shape variations and deformities observed in real fibers—as shown in Fig. 1. Deformities of this kind are presumably produced in the laboratory by the method used to cut fiber samples. Deformations may also arise from fiber degradation over time, thermal variations (warpage), abrasions and other factors not usually addressed in theoretical/numerical works. 1.1. Shape effects on fiber suspension viscosity One of the common approaches in fiber suspension theory is to assume that fibers are rigid straight cylindrical rods. This assumption greatly simplifies the theory since a rigid straight cylinder is an easily defined geometry and complications arising from shape changes during flow may be avoided. It has been shown, however, that fiber shape is an important factor in the resultant suspension viscosity [6,7]. Nawab and Mason [8] observed experimentally that Nylon fiber suspensions could produce a measured viscosity three times the order of theoretical predictions. They speculated that this difference might be a result of fiber curvature. Forgacs and Mason [9] showed that a slight curvature in the fiber would greatly alter fiber rotations. Kitano et al. [10] presented an empirical relationship between relative viscosity and volume fraction for various fillers such as glass and carbon fibers. Their results suggested that different materials would produce different suspension viscosities. It is reasonable to suppose that the varying filler stiffness (and hence shape deformation) was the cause for the viscosity variation. Likewise Goto et al. [11,12] looked experimentally at the relationship between fiber flexibility and suspension viscosity. They found that more flexible fibers produced relatively higher suspension viscosity than the equivalent stiffer fiber suspensions. Blakeney [13] found that slightly curved fibers (of the order

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

3

O(5◦ ) end-to-end fiber axis tangent angle) could increase suspension viscosity by the order 10–15% above that of a comparable straight fiber suspension. While we do not claim that all straight rigid rod theories will generally underestimate viscosity, the evidence does seem to imply that fiber deformity does contribute to a significant suspension viscosity increase over the equivalent straight fiber suspension. With this work we use a theoretical/numerical basis to reproduce and describe this viscosity–curvature relationship. 2. Numerical simulation method 2.1. Background: numerical methods for fiber suspensions Jeffrey [14] provides the basis on which most fiber simulations build. Jeffrey showed that the centroid of a lone ellipsoid suspended in shear flow will translate with a velocity equal to the region of fluid of which the centroid displaces. If the orientation of the ellipsoid’s main axis of revolution is described by a unit vector P, then P will change in a repeating cyclic ‘orbit’ commonly called the Jeffrey Orbit. If we consider the cylindrical fiber shape to be an approximation to the ellipsoid, the Jeffrey equation may be used to describe the motion of the fiber. This motion is described by Eq. (1), a2 − 1 P˙ = −ω · P + e2 [D · P − D : PPP], ae + 1

(1)

where ar is the fiber aspect ratio and ae ≈ 0.7ar is the effective aspect ratio, (ae2 − 1)/(ae2 + 1) is a shape correction factor for cylindrical rods [14], ω = (1/2)(⵱u − ⵱uT ) is the vorticity tensor, D = (1/2)γ˙ = (1/2)(⵱u + ⵱uT ), and u is the flow velocity vector. The distinguishing characteristic of this orbit is that the vector P spends a great majority of time aligned in the direction of the shear flow (x axis) broken by regular transition periods. During these transition periods the fiber rapidly flips over along the orbital path until the opposite fiber end is aligned to the shear flow. The transitions occur in regular time intervals described by Eq. (2),
 2π 1 Torbit = ae + . (2) γ˙ ae  The generalized shear rate used here is γ˙ = (1/2) Tr γ˙ 2 . The orbital path of vector P remains constant and unchanging for ‘rigid rods’ in purely Newtonian sheared fluids. While Jeffrey’s Eq. (1) is adequate for a lone rigid fiber, it is not suitable for the more industrially useful case of multiple interacting fibers in a suspension. Eq. (1) does not account for the influence of inter-fiber interactions which can drastically alter fiber motion from the Jeffrey orbit. A fiber in suspension with peers is disturbed by hydrodynamic interaction forces to the extent that the Jeffrey orbit is often totally overwhelmed. Several strategies at this point may be taken to account for interactions. The motion of an interacting fiber within a suspension may now be modeled using Jeffrey’s Eq. (1), with an additional stochastic component to model the effect of interaction forces and flow velocity disturbances from neighboring fibers. Folgar and Tucker [3] extended Jeffrey’s equation resulting in an equation similar to Eq. (3), P˙ = −ω · P + λ[D · P − D : PPP] + (δ − PP) · F(t).

(3)

4

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

The last term is an additional Gaussian random rotational diffusion force F(t) (δ is the unit tensor) which has zero mean and delta correlation function with magnitude Ci γ˙ . Ci is termed the interaction coefficient which at this time can only be determined empirically. Experimental work by Bay [15] suggests 10−3 < Ci < 10−2 , while numerical predictions by Joung et al. [1] (10−3 < Ci < 10−2 ), Yamane et al. [16] (10−5 < Ci < 10−4 ), and Fan et al. [17] (1.625 × 10−3 ) agree to varying degrees. Eq. (3) can be used to produce evolution equations for the average configurational structure tensors A2 and A4 as shown by Advani and Tucker [18] (see also [19]). These structure tensors are defined by   A2 = PPψ(P) dP = PP, A4 = PPPPψ(P) dP = PPPP, (4) where ψ is the orientation distribution function for P. Explicit knowledge of the function ψ is not required due to the ensemble averaging process in Eq. (4). Advani and Tucker’s equation for the evolution of A2 reads DA2 3 = −{ω · A2 − A2 · ω} + λ{D · A2 + A2 · D − 2D : A4 } + 2Ci γ˙ {δ − 3A2 }, (5) Dt where D/Dt denotes the material time derivative. The second-order tensor A2 is the average orientation tensor and its eigenvectors reveal the average suspension orientation and the degree of variation from that average. The second term of Eq. (5) contains the fourth-order A4 moment tensor and its presence produces the problem of closure for this equation. Several closure approximation methods have been developed that perform somewhat satisfactorily [15,20,21], however, they are deficient as none perform well in all flow regimes. Bay [15] has identified the closure approximation as the most likely cause of inaccuracy in Eq. (5). As an alternative to the Advani–Tucker Eq. (5), one could base a simulation on Eq. (3) instead. By running multiple simulations of the trajectory of a single fiber P, one after another in series and after compiling the data, a statistical average orientation may again be found. Fan et al. [22] performed this type of ‘Brownian configuration’ simulation and was able to determine the components of the A2 tensor— without the requirement of a closure approximation for A4 . 2.2. Direct simulation of fiber suspensions Both methods in Section 2.1 have stochastic components and so the fiber orientation calculated in that way is an orientation probability. Information of the motion of individual fibers in the suspension is not available in those methods. Direct simulation on the other hand is a very different method where the state of every fiber in suspension is known. In this method, the position of all suspension fibers are explicitly calculated numerically. Basic rules governing how fibers will translate, and interact with the suspending fluid and other fibers are defined. The fibers are then placed in a flow field and the position and orientation of every fiber is tracked with time. As fibers come in contact with obstacles or interact with another fiber, the rules governing the interaction behavior determine how the fiber will react in that situation. As the simulation is run, many thousands or millions of interactions may take place and eventually structural evolution on a macroscopic level may be observed. The direct numerical simulation used in this work was previously used by Joung et al. in [1]. This numerical simulation was based on the work of Fan et al. [17], however, the method was substantially modified to allow for flexible fibers instead of the standard cylindrical rigid fiber. A brief summary of

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

5

this new simulation follows, however, it is advised that interested readers consult [1] for a complete description. In this simulation method, we proceed initially as indicated by Fan et al. Multiple fibers are placed in a simulation ‘cell’ containing a shear flow field. Each fiber is described by a unit vector P. Fiber center of mass translation occurs in an affine manner. The rotation of P is described by Eq. (1), (this is another form of Jeffreys equation):   ar2 1 T ω∞  P × 2 (6) D·P− 2 D ·P . ar + 1 ar + 1 Fiber interaction forces are calculated given their current relative position and velocity, and local flow conditions. These interactions include short range lubrication forces, and viscous drag with long range hydrodynamic effects (refer to [17] for a full description of the method). From these interaction forces, Fan et al. calculate an additional deviational rotation component which is added to that found using Eq. (6). All fibers are then translated and rotated for that timestep and the process is repeated for subsequent timesteps, thus the suspension evolves in simulation over time. Fan et al. use Ewald summation in their algorithm as a means of modeling for an infinite suspension. The modified simulation of Joung et al. is very similar to Fan et al., however, the fiber model is changed and thus associated theories involving interaction forces are changed accordingly. Also an additional calculation is performed for fiber shape changes. Instead of Ewald summation to simulate an infinite suspension, duplicate cells are repeated on all sides of the simulation cell in a lattice structure (see Fig. 2). The rigid cylindrical rod of Fan et al. is replaced with an inextensible chain of spheres linked with joints centered within each sphere (see Fig. 3). The joints have an elastic rotational stiffness which attempts to straighten the joint if deflected. Internal bending moment Yνb torsion moment Yνt , (see Fig. 4) and the bend and twist angles (θ b , φ t ) may be calculated for each joint using Eq. (7), given the Young’s and shear moduli (θ eq and φ eq are the equilibrium

Fig. 2. The reference cell and periodic repetition of the cell.

6

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Fig. 3. The chain-of-beads fiber model used by Joung et al. The deformation has been exaggerated in this figure.

bend and twist angle at the joint. For straight fibers θ eq = φ eq = 0 radians): |Yνb | = −k b (θ b − θ eq ) |Yνt | = −k t (φ t − φ eq )

EI 2a GJ t when k = . 2a when k b =

(7)

Joint deflections are caused by fiber interaction forces. Fan et al. solve a slender body approximation for long range forces and a lubrication approximation for close range interaction. In this simulation we do not use slender body approximation and instead apply long range interaction forces appropriate to spheres in viscous flow, as described in depth by Karilla and Kim [23]. For a moving sphere m, viscous

Fig. 4. Joining linking connector segments.

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

7

drag force is described by Stokes Eq. (8), Fviscm = −6πηs a[Vm − V∞ ],

(8)

where a is the sphere radius, Vm is sphere velocity and V∞ is local fluid freestream velocity. Undisturbed freestream velocity V∞ = D · rm + Vm is perturbed by the long range hydrodynamic effects of other suspension spheres via the Vm term (rm is the sphere position vector). This perturbation on a sphere m is the cumulative perturbation velocity from all other spheres n in the suspension, as described by Eq. (9),  Ωmn · Fn , Vm = (9) n

where Fn is the drag force acting on all other suspension spheres n, and Ωmn is the Oseen tensor. The Oseen tensor is,  1 rmn rmn Ωmn = , (10) δ+ 8πηs |rmn | |rmn |2 where rmn is the relative position vector between spheres m and n in suspension. Short range lubrication forces are calculated using Eq. (11) Flubmn = −3πηs

a2 Vsqz , h − 2a

(11)

where Vsqz is the relative ‘squeezing’ velocity between spheres m and n given by Vsqz = [rmn /|rmn |] (rmn /|rmn | · [Vn − Vm ]) and h = |rmn |. The simulation of Joung et al. proceeds as follows. The end-to-end vector of the bead-chain is taken to be the vector P. Once the fiber rotation and translation rates are known, their new position and orientation are calculated for the current timestep as per the description above for Fan et al.’s work. Following this stage, the vector P is overlaid by the chain-of-beads. Since the interaction forces are known through Eqs. (8) and (11) and, an additional calculation is performed to determine the deflection of the bead-chain along the fiber length. Hence, at each timestep, fibers may also deform in response to flow conditions and interactions. Fiber deformation effectively results in another small rotation in the end-to-end orientation vector P. Any bulk properties observed to change with stiffness moduli may thus be attributed to flexibility effects. Finally a stress rule is required to calculate the bulk properties given the suspension microstructure. The total stress tensor is σbulk = σsolvent + τfiber .

(12)

The solvent stress component is σsolvent = −p δ + ηs · γ˙ ,

(13)

where p is hydrodynamic pressure. The remaining particle contributed stress τfiber , is appropriately expressed for fiber suspensions by the transversely isotropic fluid equation (TIF) developed by Ericksen and Hand [24,25]. An improved version by Phan-Thien and Graham [26] (TIF-PTG) is τ = 2ηs {D + f (φ, ar )D : PPPP},

(14)

8

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Fig. 5. Curved chain-of-bead fibers. Curve angle from 2◦ < θfiber < 40◦ .

where f (φ, ar ) =

ar2 φ(2 − (φ/G)) , 4(ln(2ar ) − 1.5)(1 − (φ/G))2

G = 0.53 − 0.013ar ,

5 < ar < 30.

From the stress tensor we find suspension relative viscosity ηbulk and the first and second normal stress difference coefficients, Ψ1 , Ψ2 from the relations τyx = ηbulk γ˙yx ,

τxx − τyy = Ψ1 (γ˙yx )2 ,

τyy − τzz = Ψ2 (γ˙yx )2 .

(15)

2.3. Rigid curved fiber simulation While in [1] the joints along the bead-chain were allowed to bend and twist in time, for this study the fiber is rigid and so the algorithms involved with joint bending are disabled and bend angle fixed. Bend angle at each joint along the

fiber was set equally such that the accumulated angle from one end to the other of beads−1 ◦ ranged from 2 < [θfiber = no. θjoint ] < 40◦ . Fibers of varying curvature are shown in Fig. 5. joint=2 2.3.1. Appropriate use of the TIF stress equation on curved fibers In justifying the appropriateness of using the TIF stress equation (14) on these curved fibers, we appeal to the work of Hinch [27] for flexible, inextensible threads. Hinch presents Eq. (16) for thread tension:
 s   s+L L 2 2 1 T = P · E · P 2 (L − s ) − 2εP · E · y− y , (16) 2L −L −L where using the notation of Hinch, T is thread tension, L the half of the fiber length, s the length variable along the fiber centreline, ε the magnitude of the deflection and y(s) is the shape function of the thread as a function of s. Let us describe the arc of the curved fiber in this simulation with the fundamental mode shape y(s) = ε cos((π/2L)s). If the fundamental mode shape is substituted into Eq. (16) and integrated, we get for the

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

9

second term, (2L/π) sin ((π/2L)s) − (2s/L). Since the stress contribution from the fiber is proportional to the integral of tension along the entire fiber length, we may substitute s = L into the resulting integral. The integral goes to zero, and hence the entire second term of Eq. (16) goes to zero. The first order O(ε) effects on tension (and hence stress) for the curved fiber is zero. Therefore, the errors associated with using TIF equation on these curved fibers are sufficiently small as to be considered negligible (second order in deflection magnitude, O(ε2 ) or smaller). While errors associated with the use of TIF may be acceptably small it is, however, acknowledged that many other theories used for this work were also intended only for straight rods. Therefore, as in [1] our confidence in the results decreases as fibers increasingly deviate from dead straight. In the following results of Section 3 one should note that as curvature increases our estimate for viscosity will become inaccurate. This should be kept in mind particularly while reviewing Figs. 8–10. However, by the same token the inaccuracies reduce to zero with curvature and so the pronounced effects shown in Figs. 8–10 for low curvature are not invalidated by a ‘nearly straight’ fiber assumption. 2.4. Conditions The main conditions imposed on this simulation are, • • • •

the solvent is a Newtonian fluid; velocity gradients are homogeneous; beads and fibers are of dimensions such that Brownian effects are negligible; the suspension is assumed to meet conditions set by Batchelor [28] such that an ensemble average of the fiber population is a valid procedure to produce the required bulk properties; • the effects of fiber orientation and concentration on the flow field are not considered, i.e. flow field and fiber orientation are de-coupled on the macro-scale; • beads and fibers are force and torque free at all times. Flow rates are sufficiently slow that inertial effects are negligible; • all fibers are identical in connectivity-configuration and shape. To elaborate on the fourth condition, Batchelor ([28], Section 3) assumes the state of some quantity can be adequately estimated from a single sufficiently extensive realization. Hence, an ensemble average over some region in space may be taken to be equivalent to the integral average of the same quantity over any spatial coordinate with respect to which the quantity is statistically stationary. Therefore, we will assume that the dimensions of our reference simulation cell are sufficiently large, and the number of fibers used in simulation are sufficiently great enough to allow valid averages for stress, viscosity or any other quantity derived from averaging. Furthermore, we will assume that any quantity derived by averaging within our reference cell is representative over the entire infinite suspension since the quantity is assumed statistically stationary in space. On the second and fifth conditions, it has been shown that the presence of particles in a flowing fluid will alter the development of the flow field [29]. In this simulation we assume that fluid velocity field is unaffected by the motion and orientation of fibers. The errors associated with this assumption would presumably decrease as fiber concentration is reduced. In light of Lipscomb et al.’s [29] work this assumption, however, is not entirely correct. We have not considered these effects at this time and thus will also assume that the main focus of this paper, the viscosity–curvature phenomenon of Section 3 is independent of the flow-fiber de-coupling assumption.

10

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

2.5. Simulation initial conditions As shown in Fig. 3, one end of the fiber is arbitrarily chosen as the base and beads and connectors are numbered from that base. Parameters used throughout this work are chosen such that all results are non-dimensional, all lengths are relative to the bead radius a = 1. The basic suspension information is established before commencing. These include terms such as the number of fibers, number of beads, unit bead radius a = 1, fiber length, and type and magnitude of flow, e.g. shear with shear rate γ˙ , or extension with extensional rate "˙ . The initial fiber suspension is then placed in the reference cell. Fibers are oriented randomly and initial fiber velocity and rotation rate is set to be equal to that of the freestream fluid in the locality of the fiber. Care is taken to ensure that no fibers are overlapping before commencing the simulation. For shear flow the simulation is typically run for, 0 ≤ γ ≤ 1000. The increment is set at dγ = 0.01, with a unit shear rate of γ˙ = 1. At each time step information such as the fiber center of mass, end-to-end vector P, global bead positions rm , bead velocity Vm , and external force Fm are determined. From this information the suspension orientation state is known and hence the bulk properties can be calculated. When a fiber translates beyond the reference cell it is repositioned such that it continues its motion from the other side of the cell. The fiber re-enters the cell immediately after it has exited, hence maintaining the suspension volume fraction. The simulation was written in Fortran77 and was compiled and executed on a Digital Alphastation 500 in the SyDCom 1 computing facility. Execution time for 100,000 timesteps ranged from several hours to several days or weeks depending on the volume fraction and periodic depth.

3. Results In the first section we verify that this simulation predicts the expected theoretical result for straight rigid fibers. We then present results for fiber suspensions where all fibers in the suspension are increasingly curved. Viscosity data is presented as a function of curvature for varying concentrations. The fibers used are always of aspect ratio ar = 16.9 to allow comparison to Bibbo’s experimental data and Fan et al.’s numerical data. The number of beads in a fiber is always NE = 16, and bead radius a = 1. The suspension concentration is controlled by varying the number of fibers within the reference cell. The suspension is sheared at the rate γ˙ = 1 with solvent viscosity normalized to ηs = 1. End-to-end curvature of the fibers range from 2◦ ≤ θfiber ≤ 40◦ . Although our use of Eq. (6) is increasingly inaccurate at curvatures as large as 40◦ , results at this curvature magnitude are included in the results for completeness. While the viscosity trends as curvature increases to 40◦ may be interesting for speculative purposes, the reader should be cautious about making definite conclusions at these larger angles. 3.1. Verification: straight fiber suspension viscosity Fig. 6 compares the current simulation relative viscosity results for varying volume fraction (䉬) against an equivalent theoretically based numerical simulation by Fan et al. [17] (夹). On the whole, the two theoretically based methods predict quantitatively similar results and therefore are a close match. 1

Sydney University Distributed Computing Facility.

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

11

Fig. 6. Experimental versus numerical simulation comparison. Bibbo (䊊) and Milliken et al. () results are experimental while Fan et al. (夹) and the current results (䉬) are numerical simulation/theoretical results.

Comparing the current results (straight ‘Rigid Sim’ (䉬)), it predicts lower viscosity than the experimental results—particularly in the concentrated region φ > 0.0592. In the semi-concentrated regions 0.0035 < φ < 0.0592, the current simulation using the chain-of-beads model predicts a slightly higher viscosity than Fan et al.’s continuous straight rod simulation. Bibbo’s viscosity results (䊊) were found using the Couette device (shear flow). Milliken et al.’s () viscosity may be slightly elevated due to the non-shear flow kinematics associated with the falling ball method used to obtain these results. The difference between the current simulation and Fan et al.’s results is most pronounced in the range 0.02 < φ < 0.05. There may be several reasons for this and are believed most likely to be due to the different numerical methods used. The current simulation uses discrete interaction points (spheres) along the fiber length while Fan et al. uses a continuous slender rod. There are also variations in the method used to simulate infinite suspensions. As described earlier, Fan et al. use the Ewald summation while for this simulation, cell repetition is used. Ewald summation is acknowledged to be the more accurate representation of infinite suspensions. Any number of methods may be used of varying complexity and realism. Both methods, however, agree well with each other quantitatively, and the qualitative trend if not identical, is also acceptably close. For the conditions used in this paper we consider these results sufficiently similar to Fan et al. and other previous works. In summary this simulation performs well against experimental data, and quantitatively well against the numerical results of Fan et al. We shall now proceed to curved fiber suspensions. 3.2. Curved fiber suspension—ηbulk sensitivity to curvature Fig. 7 compares the same ‘straight’ fiber simulation result against progressively more curved fiber simulation results.

12

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Fig. 7. Simulation results for varying fiber curvature. Curved fiber suspensions (hollow symbols) consistently result in higher viscosity than straight fibers (䊏).

This result (Fig. 7) is in agreement with the earlier premise that curvature in the fiber is a significant contributor to the elevation in viscosity over straight fiber suspensions. Over the entire volume fraction range and for all curvatures tested, there is a large increase in the viscosity. Further, the curved fiber viscosity magnitude is closer to those produced by Bibbo (䊊) in the laboratory (compare with Fig. 6). Fig. 7 does not reveal the relationship between suspension viscosity and fiber curvature. In Fig. 8 we plot relative viscosity against fiber curvature for varying suspension concentrations. The viscosity data presented features a peak for rather small curvature angles in the range 5◦ < θfiber < 10◦ . Beyond this, viscosity reduces in magnitude and appears to converge to a value much larger than for straight fibers (although we can expect results to become inaccurate at large curvatures). Fig. 8 reveals a peculiar phenomenon in which the viscosity rises rapidly and by a great amount for very small curvature angles. The curve gradient of Fig. 8 at zero fiber curvature is non-zero and in fact quite steep. We also note that the proportional increase in viscosity as curvature deviates from ‘dead-straight’ is more pronounced for higher volume fractions. For example the highly concentrated fiber suspension of volume fraction φ = 0.1562 (symbol in Fig. 8) has a relative viscosity of ηr ≈ 3.5 when fiber curvature is θfiber = 0◦ (dead straight) but for a small curvature of θfiber = 4◦ , the relative viscosity jumps to ηr ≈ 7.7. The relative viscosity is of the order two times greater than the ‘dead-straight’ result and for the same curvatures this order increase is approximately consistent over all the tested volume fractions. This magnitude of viscosity increase is reminiscent of Nawab and Mason’s [8] observed experimental viscosity of order three times the theoretical predictions (Section 1). Fiber deformations may very well have played an important part in their result. This result demonstrates the extreme sensitivity of the suspension viscosity to fiber shape imperfections. For a seemingly small shape distortion, the effect on viscosity is dramatic. The suspension viscosity

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

13

Fig. 8. Suspension relative viscosity versus fiber end-to-end curvature, for varying suspension volume fractions. The fibers are rigid with ar = 16.9, Number of beads per fiber = 16. The suspension contains a population of 100% identically curved fibers.

appears to reach a maxima at a very small curvature. Beyond this critical curvature angle the viscosity reduces markedly again. This behavior would pose the following problem to experimental researchers attempting to produce accurate reproducible viscosity results. First, great care must be taken so that fiber samples are well prepared. Tolerance to bent fibers and other deformation is low and consequences of using imperfect fibers are severe. Yet as we can see from a snapshot of a typical suspension by Milliken et al. [30] (Fig. 1) small curvatures are not uncommon despite the care taken by Milliken in preparing fiber samples. Second, results cannot be much improved by simply removing the most extremely bent fibers from the sample. As we have seen, the greatest effects occur for curvatures of only 5◦ < θfiber < 10◦ . One would have to undertake the much more arduous task of removing all fibers with curvature beyond say, 1◦ . Last, one can expect more flexible fiber materials to deform more than stiffer ones. Therefore, for consistency, stiffer fibers are preferable to more flexible fibers. Likewise for fibers of greater aspect ratio (whereby the tendency to flex is greater than for shorter stout fibers) stiffer fibers are preferable when consistently reproducible viscosity is desirable. 3.3. Viscosity-curvature curve fitting The data presented in Fig. 8 cannot be fitted easily by a low order least square polynomial fit. Instead it is better described with a combination of hyperbolic and F (x)exp functions (where F (x)exp is of the form F (x)exp = x e−x ). Fig. 9 shows the two individual components and the resultant sum of these two parts. The resulting empirical equation fitted to the simulation data is ηbulk (θfiber , φ) = "θfiber e(−αθfiber −ζ ) + δ

eβθfiber − e−βθfiber + γ, eβθfiber + e−βθfiber

(17)

14

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Fig. 9. The F (x)exp and hyperbolic components, and the combined result.

where the coefficients, α, β, γ , δ, " and ζ are functions of volume fraction φ, 0 ≤ φ ≤ 0.16, 0.0007 + 1.2φ + 0.005, β(φ) = 4.2828φ 2 − 0.1138φ + 0.0022, φ γ (φ) = 8904.6471φ 4 − 1716.2639φ 3 + 118.8461φ 2 + 5.2429φ + 0.9922, δ(φ) = −747.8328φ 3 + 237.5086φ 2 − 8.4300φ + 0.3531, 0.01 "(φ) = − 1.5φ + 2, ζ (φ) = 153.3682φ 2 − 61.5718φ + 5.8786. φ α(φ) =

Fig. 10. Curve fitted equation plotted with simulation results (ar = 16.9).

(18)

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

15

Eq. (17) fits most of the simulation data well as can be seen in Fig. 10. One should note that fiber aspect ratio and bulk properties are related. Therefore, the constants presented in Eq. (18) are applicable for fibers with aspect ratio ar = 16.9 only. The actual form of the viscosity–curvature relationship for all fiber aspect ratios and curvatures is yet to be found. This curve fit is presented simply as a glimpse into the true nature of this relationship. We may expect that at any time in shear flow a certain proportion of fibers are permanently bent. A certain percentage may also acquire bends due to mechanical interaction or thermal gradients (warpage). The steep gradient of the viscosity–curvature curve at low curvature angle reveals the sensitivity of suspension viscosity to even the slightest curvature in fibers. The wide scatter in existing experimental viscosity data may well be caused by minor fiber shape imperfections.

4. Discussion and conclusion This work addresses the importance of fiber shape in sheared fiber suspensions. The initial hypothesis is that small deformities in the fibers of a suspension produces a large increase in suspension viscosity. In this work we directly simulate fibers in Newtonian shear flow for varying volume fraction and varying fiber curvature θfiber , and the bulk viscosity is observed. It is found that there is a rapid and large increase in bulk viscosity as fiber curvature is increased. The results of viscosity and curvature can be fitted to a single Eq. (17) given volume fraction and fiber curvature (for fiber aspect ratio ar = 16.9 only). The curious form of the curvature-viscosity relationship (Fig. 10) leads to the obvious question of why this happens. The answer is unknown, but a speculative discussion is offered. It is acknowledged that fiber alignment is central to the resulting viscosity of a fiber suspension. Generally speaking, fiber misalignment increases interaction activity and encourages further misalignment. The resulting suspension is more resistant to flow and so a viscosity increase is observed. When a single infinite aspect ratio fiber has aligned to shear direction it presents no profile in the shear direction. Theoretically it will remain aligned indefinitely. A curved fiber, however, does present a profile in the shear direction. The more curved a fiber, the larger the profile seen in the shear direction. Therefore, regardless of aspect ratio it does not remained aligned indefinitely. Curved fibers in suspension will not exhibit as strong a tendency to remain aligned to shear flow as does a straight fiber suspension. Curved fibers are less able to pack neatly into dense aligned structures than straight fibers. Therefore, a curved fiber suspension has a propensity to misalignment and hence we see a viscosity increase with curvature. This cannot, however, explain the subsequent viscosity decrease beyond approximately θfiber = 10◦ curvature angle. It is yet to be seen whether the behavior in that region also reflects a physical phenomenon or is merely the result of theoretical or numerical inadequacy. The implications are of particular importance to experimenters in fiber suspensions. It would appear that the viscosity (and possibly other bulk properties) are highly sensitive to the quality control maintained in preparing the fiber samples. The rapid increase in bulk viscosity as seen in Fig. 10 shows that even a small deviation from ‘dead straight’ fibers, or the use of highly flexible fibers, or highly elongated fibers (large aspect ratio fibers are more prone to bending) will result in large deviations in bulk property. Further, one simply cannot remove the most deformed from the fiber samples in the hope of reducing this viscosity increase. It appears from Fig. 10 that the fibers which contribute most to the increase are those with a relatively mild curvature of the order O(5–10◦ ).

16

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

Acknowledgements This research is supported by the Australian Research Council (ARC) through a collaborative research grant with Moldflow Pty Ltd. The numerical results were produced on the Sydney University Distributed Computing facility (SyDCom). We wish to thank Professor Roger I. Tanner and Dr. Howard See for their helpful comments and suggestions. References [1] C.G. Joung, N. Phan-Thien, X.J. Fan, Direct simulation of flexible fibers, J. Non-Newtonian Fluid Mech. 99 (2001) 1. [2] G. Potsch, W. Michaeli, Injection Molding—An Introduction, Hanser Publications, 1995. [3] F. Folgar, C.L. Tucker, Orientation behavior of fibers in concentrated suspensions, J. Reinforced Plastics and Composites 3 (1984) 98–119. [4] R. Zheng, P. Kennedy, N. Phan-Thien, X.-J. Fan, Thermoviscoelastic simulation of thermally and pressure induced stresses in injection moulding for the prediction of shrinkage and warpage for fibre-reinforced thermoplastics, J. Non-Newtonian Fluid Mech. 84 (1999) 159. [5] R. Zheng, N. McCaffrey, K. Winch, H. Yu, P. Kennedy, Predicting warpage of injection molded fibre-reinforced plastics, J. Thermoplastic Composite Mater. 9 (1996) 90. [6] E. Anczurowski, S.G. Mason, The kinetics of flowing dispersions: II. Equilibrium orientations for rods and discs (theoretical), J. Colloid and Interface Sci. 23 (1967) 522. [7] F.P. Bretherton, The motion of rigid particles in a shear flow at low Reynolds number, J. Fluid Mech. 14 (1962) 284. [8] M.A. Nawab, S.G. Mason, Viscosity of dilute suspensions of threadlike particles, J. Phys. Chem. 62 (1958) 1248. [9] O.L. Forgacs, S.G. Mason, Particle motions in sheared suspensions: X. Orbits of flexible threadlike particles, J. Colloid Sci. 14 (1959) 473. [10] T. Kitano, T. Kataoka, T. Shirota, An empirical equation of the relative viscosity of polymer melts filled with various inorganic fillers, Rheol. Acta 20 (1981) 207. [11] S. Goto, H. Nagazono, H. Kato, The flow behavior of fiber suspensions in Newtonian fluids and polymer solutions: 1. Mechanical properties, Rheol. Acta 25 (1986) 119. [12] S. Goto, H. Nagazono, H. Kato, The flow behaviour of fibre suspensions in Newtonian fluids and polymer solutions: II. Capillary flow, Rheol. Acta 25 (1986) 246. [13] W.R. Blakeney, The viscosity of suspensions of straight, rigid rods, J. Colloid and Interface Sci. 22 (1966) 324. [14] G. Jeffrey, The motion of ellipsoid particles immersed in a viscous fluid, Proc. R. Soc. London 102 (1923) 161. [15] R.S. Bay, Fibre orientation in injection mol ded composites: a comparison of theory and experiment, Ph.D. thesis, Department of Mechanical Engineering, University of Illinois, IL, USA, 1991. [16] Y. Yamane, Y. Kaneda, M. Doi, Numerical simulation of semi-dilute suspensions of rod-like particles in shear flow, J. Non-Newtonian Fluid Mech. 54 (1994) 405. [17] X. Fan, N. Phan-Thien, R. Zheng, A direct simulation of fibre suspensions, J. Non-Newtonian Fluid Mech. 74 (1998) 113. [18] S.G. Advani, C.L. Tucker, The use of tensors to describe and predict fibre orientation in short fibre composites, J. Rheol. 8 (1987) 751. [19] N. Phan-Thien, R. Zheng, Macroscopic modelling of the evolution of fibre orientation during flow, in: D. Guell, T. Papathanasiou (Eds.), Flow-Induced Alignment in Composite Materials, Woodhead Publishing Cambridge, UK, 1997, Chapter 3, pp. 77–111. [20] S.G. Advani, C.L. Tucker, Closure approximations for three-dimensional structure tensors, J. Rheol. 34 (1990) 367. [21] J.S. Cintra, C.L. Tucker, Orthotropic closure approximations for flow-induced fiber orientation, J. Rheol. 39 (1995) 1095. [22] X. Fan, N. Phan-Thien, R. Zheng, Simulation of fibre suspension flows by the Brownian configuration field method, J. Non-Newtonian Fluid Mech. 84 (1999) 257. [23] S. Kim, S. Karilla, Microhydrodynamics: Principles and Selected Applications, Butterworth Series of Chemical Engineering, Butterworths, London, 1991. [24] J.L. Ericksen, Anisotropic fluids, Arch. Rational Mech. Anal. 4 (1960) 231. [25] G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. 33 (1961) 33–46.

C.G. Joung et al. / J. Non-Newtonian Fluid Mech. 102 (2002) 1–17

17

[26] N. Phan-Thien, A.L. Graham, A new constitutive equation for fibre suspensions: flow past a sphere, Rheol. Acta 30 (1991) 44. [27] E.J. Hinch, The distortion of a flexible inextensible thread in a shearing flow, J. Fluid Mech. 74 (1976) 317. [28] G.K. Batchelor, The stress system in a suspension of force-free particles, J. Fluid Mech. 41 (1970) 545. [29] G.G. Lipscomb, M.M. Denn, D. Hur, D. Boger, The flow of fibre suspensions in complex geometries, J. Non-Newtonian Fluid Mech. 26 (1988) 297. [30] W.J. Milliken, M. Gottlieb, A.L. Graham, L.A. Mondy, R.L. Powell, The viscosity–volume fraction relation for suspensions of rod-like particles by falling ball rheometry, J. Fluid Mech. 202 (1989) 217. [31] Milliken, et al., Electron micrographs of rayon fibers, J. Fluid Mech. 202 (1989) 217–232.