An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics

J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

Contents lists available at ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm

An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics Jay H. Phelps, Charles L. Tucker III ∗ Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, 1206 West Green Street, Urbana, IL 61801, United States

a r t i c l e

i n f o

Article history: Received 9 May 2008 Received in revised form 14 August 2008 Accepted 14 August 2008 Keywords: Fiber orientation Rotary diffusion Long-fiber thermoplastics Injection molding

a b s t r a c t The Folgar–Tucker model, which is widely-used to predict fiber orientation in injection-molded composites, accounts for fiber–fiber interactions using isotropic rotary diffusion. However, this model does not match all aspects of experimental fiber orientation data, especially for composites with long discontinuous fibers. This paper develops a fiber orientation model that incorporates anisotropic rotary diffusion. From kinetic theory we derive the evolution equation for the second-order orientation tensor, correcting some errors in earlier treatments. The diffusivity is assumed to depend on a second-order space tensor, which is taken to be a function of the orientation state and the rate of deformation. Model parameters are selected by matching the experimental steady-state orientation in simple shear flow, and by requiring stable steady states and physically realizable solutions. Also, concentrated fiber suspensions align more slowly with respect to strain than models based on Jeffery’s equation, and we incorporate this behavior in an objective way. The final model is suitable for use in mold filling and other flow simulations, and it gives improved predictions of fiber orientation for injection molded long-fiber composites. © 2008 Elsevier B.V. All rights reserved.

1. Introduction There is considerable value in being able to predict the fiber orientation state in a flowing fiber suspension. Example uses are modeling the processing of fiber-reinforced plastics, and models of paper making. An accurate treatment of fiber orientation is also an essential element in any rheological model of a fiber suspension. For non-Brownian fibers, the model proposed by Folgar and Tucker [1] has been widely used. This phenomenological model uses Jeffery’s single-fiber equation to represent flow effects [2], and adds an isotropic rotary diffusion term, with a diffusivity proportional to the scalar rate of deformation, to represent fiber–fiber interactions. Adding the rotary diffusion term gives steady orientation states that are not perfectly aligned, which is consistent with experimental results. Rotary diffusion is proportional to the interaction coefficient CI , a scalar parameter that is adjusted to fit experimental data. The Folgar–Tucker model was originally formulated as a Fokker–Planck type equation for the probability density function (PDF) of fiber orientation. However, numerically calculating the

∗ Corresponding author. Tel.: +1 217 244 3822. E-mail addresses: [email protected] (J.H. Phelps), [email protected] (C.L. Tucker III). 0377-0257/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2008.08.002

PDF for tens of thousands of nodes in an injection molding simulation is prohibitively expensive. Recasting the model in terms of the second-order moment tensor of the orientation distribution [3] makes the calculation affordable for engineering purposes. In the years since its introduction, there have been few major modifications to the Folgar–Tucker model. Experimental evidence has shown that the kinetics of orientation change in concentrated suspensions is much slower than Jeffery-based models predict [4,5], and an objective model to account for this has recently been proposed [6]. That model primarily modifies the Jeffery terms in the model. Another type of modification, which is the subject of this paper, concerns the rotary diffusion term. Rotary diffusion models are interesting for both practical and theoretical reasons. On the practical side, the Folgar–Tucker model has proved useful for modeling fiber orientation in short-fiber/thermoplastic composites (SFTs) that are processed by injection molding, e.g. [7,8]. However, for long-fiber thermoplastics (LFTs) the Folgar–Tucker model is less accurate. LFT pellets are usually prepared by a pultrusion process, and the pellets typically have 10–13 mm long fibers aligned along their length. This is a contrast to SFTs, where the fiber length is typically 0.2–0.4 mm. When LFTs are injection molded they develop fiber orientation patterns that are qualitatively similar to SFT orientation patterns, but quantitatively they achieve lower flow-direction alignment than SFTs. Accurately predicting the fiber orientation in LFT composites requires a better model for fiber ori-

166

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

entation, and it appears that modifying the rotary diffusion term could provide the requisite model behavior. On the theoretical side, one modification to the rotary diffusion term was proposed by Ranganathan and Advani [9]. They suggested that diffusivity should be dependent on the orientation state, being higher when the fibers are randomly oriented and lower when the fibers are aligned. They proposed a phenomenological model in which the interaction coefficient is inversely proportional to the average interfiber spacing. Because their rotary diffusion is isotropic, their model gives the same steady-state distributions of fiber orientation as the Folgar–Tucker model; only the transient approach to steady state is affected. Also, their model reaches an upper limit in volume fraction and aspect ratio when the average interfiber spacing becomes negative. LFT materials fall well beyond this limit, so Ranganathan and Advani’s model cannot be applied to LFT materials. A different modification of rotary diffusion was proposed by Fan et al. [10] and Phan-Thien et al. [11] who performed direct numerical simulations of multiple interacting fibers in a concentrated suspension undergoing simple shear flow. They could not fit their calculated steady-state orientation distribution using the Folgar–Tucker model, so they developed a rotary diffusion model in which the scalar interaction coefficient CI was replaced by a secondorder tensor C. This makes the rotary diffusion anisotropic. Fan and Phan-Thien developed a moment-tensor equation for their model, and found values of the C tensor to fit their steady orientation states. Significantly, they did not explore the dynamic properties of their moment-tensor equation. We build on their approach in this work. Koch [12] developed a detailed, mechanistic model of long-range hydrodynamic interactions in semi-concentrated fiber suspensions, and from this produced a model for the rotary diffusion. Like the Folgar–Tucker model, the diffusivity scales with the rate of deformation. However, in Koch’s model the diffusivity varies with the orientation state, and is anisotropic. No moment-tensor form of this model was developed, so the properties of this model have been largely unexplored. In this paper we explore anisotropic rotary diffusion models for fiber orientation. The treatment is phenomenological. In Section 2, after reviewing some background material, we compare predictions of the Folgar–Tucker model to fiber orientation measurements in molded LFT samples, showing what types of model improvements are needed for LFTs. We then examine the Fan and Phan-Thien model, showing that it fails a simple test for correct diffusive behavior, and the Koch model, showing that it has limited usefulness for modeling LFTs. With these items as motivation, in Section 3 we derive the correct moment-tensor form for an anisotropic rotary diffusion (ARD) model that is characterized by a second-order tensor. The dynamic behavior of this model proves to be very sensitive to the choice of model parameters, so we present a systematic procedure for choosing robust parameter sets. Finally, in Section 4 we show that the new ARD model gives improved fiber orientation predictions for LFT materials. The paper closes with a short discussion.

2. Background and motivation The orientation of a single, straight fiber can be characterized by a unit vector p directed along the fiber axis. For collections of fibers, the complete description of orientation is a probability density function (p) where the probability of a fiber being oriented between p and p + dp equals (p) dp. Experimental measurements of orientation are samples drawn from the distribution function . A useful way to report the data is to compute the average values of certain functions of p. The second-order orientation tensor A [3]

is a general and concise way to describe orientation. This tensor is given by



A=

pp dp,

(1)



where dp denotes an integral over all possible fiber orientations, and pp is the tensor product of the fiber orientation vector p with itself. The primary impetus for using A instead of to describe orientation comes from the numerical simulation of fiber orientation in complex flows. It is computationally more efficient to track the components of A over multiple time steps and many spatial nodes, than it is to track the distribution function over the same discretization. 2.1. Folgar–Tucker fiber orientation model Turning attention to modeling the evolution of orientation, we first review the concept of orientation space. Since p is a unit vector, the space of all possible orientations (the orientation space) is the surface of a unit sphere. The probability density function is defined on this space. The probability density is conserved (as the total probability is constant), and the probability density function must satisfy a continuity equation at any time t. This continuity equation is given by d ˙ = −∇ s · ( p), dt

(2)

where s represents the gradient operator on the surface of the unit sphere, or the gradient operator in orientation space. The probability flux term p˙ may be decomposed into a hydrodynamic term h p˙ and a diffusive flux vector q, giving d h = −∇ s · ( p˙ + q). dt

(3)

Typically, the hydrodynamic contribution is modeled by Jeffery’s equation for particle motion in a dilute suspension [2], or h p˙ = W · p + (D · p − D : ppp),

(4)

where W = (1/2)(L − LT ) is the vorticity tensor, D = (1/2)(L − LT ) is the rate-of-deformation tensor, and  is the particle shape parameter (for long, slender fibers  → 1). L represents the velocity gradient tensor of the dilute suspension, with components Lij = ∂vi /∂xj where vi represents the component of the velocity in the xi direction. In order to model the diffusion term appearing in Eq. (3), Folgar and Tucker [1] suggest the phenomenological relationship: q = −CI ˙ ∇ s ,

(5)

where CI is the fiber–fiber interaction coefficient (a fitting parameter), and ˙ = (2D : D)1/2 is the scalar magnitude of D. In this model, Folgar and Tucker scale q with , ˙ assuming that the rate of fiber–fiber collisions is proportional to , ˙ and the orientation perturbation per collision is independent of . ˙ Advani and Tucker [3] utilize Eqs. (1) and (3)–(5) to develop a time evolution equation for A. This is the standard Folgar–Tucker model, and is given by h d A˙ = A˙ + A˙ ,

(6)

where the hydrodynamic contribution is h A˙ = (W · A − A · W) + (D · A + A · D − 2A : D),

(7)

and the diffusive contribution is d A˙ = 2CI (I ˙ − 3A).

(8)

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

In this model, A˙ represents the material derivative of A, I is the identity tensor, and A is the fourth-order orientation tensor, defined as



A=

pppp dp.

model, and it is given by d h = −∇ s · [ p˙ − CI ˙ ∇ s dt

(9)

+

3 

wi (ei ei · p − ei ei : ppp)],

(15)

i=1

h

In practice, A is approximated as a closure function of components of A. The present work utilizes the ORE orthotropic closure. The eigenvalue-based orthotropic closure family was introduced by Cintra and Tucker [13], and the coefficients for the ORE version are given by VerWeyst [14]. d The A˙ term approximately accounts for the fiber–fiber interactions in a concentrated suspension. This interaction term is a rotary diffusion term. It acts isotropically, providing a flux that pulls fibers toward a random state. 2.2. The RSC fiber orientation model Recent experience has shown that the rate of orientation development in SFT materials is much slower than this model predicts [4,5]. Wang et al. [6] introduce the reduced-strain closure (RSC) model, which slows the orientation kinetics in an objective fashion, in order to achieve better agreement between experiment and prediction. To develop the RSC model, Wang et al. [6] write the symmetric orientation tensor as a function of its eigenvalues i and 3  e e . This expression may then be subeigenvectors ei , A = i=1 i i i stituted into the Folgar–Tucker model, and the material derivatives of both the eigenvalues ˙ i and the eigenvectors e˙ i may be isolated. The authors then make the phenomenological assumption that the eigenvalue kinetics for the RSC model are slowed by a factor of , compared to the eigenvalue kinetics derived from the Folgar–Tucker model, while the expressions for the eigenvector kinetics remain unchanged. This is given by ˙ RSC = ˙ FT , i i

(10)

and RSC FT e˙ i = e˙ i .

(11) RSC

Using Eqs. (10) and (11), the tensor material derivative A˙ re-constructed. The result is the RSC model, given by

may be

RSC = (W · A − A · W) + {D · A + A · D A˙

− 2[A + (1 − )(L − M : A)] : D} + 2CI (I ˙ − 3A).

(12)

The fourth-order tensors L and M are analytical functions of the eigenvalues and eigenvectors of A, and are given by L=

167

3 

i (ei ei ei ei ),

(13)

i=1

and M=

3 

ei ei ei ei .

(14)

I=1

The  parameter, which slows the orientation kinetics, must be chosen to fit experimental data. Typical values of  for SFT materials range from 0.05 to 0.2 [6]. Although the RSC model was developed at the level of a moment-tensor equation, there is an equivalent model at the level of kinetic theory. Wang et al. [6] also derive such a kinetic theory

where p˙ is Jeffery’s equation for particle motion. The term with the summation represents an orientation-dependent probability flux. The scalar coefficients wi are calculated as wi = −

5(1 − ) ˙ − 3i )], [(i D : ei ei − ei ei : A : D) + CI (1 4

(16)

where i and ei are once again the eigenvalues and eigenvectors of the orientation tensor A. This kinetic theory RSC model corresponds to Eq. (12) in the same manner that Eq. (3) corresponds to the Folgar–Tucker moment-tensor model of Eqs. (6)–(8). While the RSC-equivalent kinetic theory model does contain the moment tensors A and A, these tensors may be evaluated at each time step from their definitions, so no closure approximation is needed to solve Eq. (15) [6]. For the remainder of this work, we return to the moment-tensor form of the RSC model. 2.3. Application to LFT injection molding The Folgar–Tucker model family has been implemented in almost all commercial software programs for injection molding. For research purposes we have made a similar implementation in a research software package known as ORIENT. This software and its implementation have a structure similar to commercial injection molding codes. A detailed description of the theory and numerical methods behind this program is found in Bay and Tucker [7]. A brief summary follows here. ORIENT is a finite difference program that calculates the velocity field and second-order orientation tensor A during mold filling for either an end-gated strip or a center-gated disk. ORIENT also takes into account non-isothermal conditions and solves for the temperature field within the mold. ORIENT uses the Hele–Shaw approximation [15] for solving for the velocity field in mold-filling operations, and accounts for the dependence of viscosity on both strain rate and temperature. In solving Eq. (6) or (12), ORIENT needs a boundary condition for orientation at the mold inlet. Typically, experimental data taken near the inlet is used as a boundary condition for each orientation tensor component. ORIENT accepts values of , CI , and  as input. ˙ The outputs of ORIENT include vx (flow direction velocity), T, , and A as functions of x (flow direction coordinate) and z (thickness direction coordinate), and pressure as a function of x, all at each time step during filling. Typically, we only use the data for the last time step, when the cavity is full. There is no post-filling analysis. To obtain experimental LFT data against which to compare ORIENT calculations, a series of glass-reinforced (40 wt%) specimens with polypropylene matrices were molded using a range of injection speeds. Two geometries were considered: a 90 mm long × 80 mm wide × 3 mm thick end-gated strip, also called an ISO plaque [16], and a 180 mm diameter, 3 mm thick center-gated disk. After molding, the number-average fiber length (L¯ n ) was 2.866 mm, and the weight-average fiber length (L¯ w ) was 5.167 mm, for a slow-filled ISO plaque specimen [17]. Experimental orientation measurements were performed on polished cross-sections using the image analysis methods developed by Hine et al. [18]. Further details of these moldings and orientation measurements can be found in Nguyen et al. [17]. We summarize some key findings here. Let the A11 orientation direction correspond to the flow direction during mold filling, while the A22 and A33 components correspond

168

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

steady-state orientation in simple shear flow. CI is usually chosen to match A11 in this region. Typical CI parameter values for SFT composites are 0.006–0.01; the best-fit value for a given data set will depend on which closure approximation is used in the model. For the LFT composites studied here,  and CI were set to 1/30 and 0.03, respectively, in order to fit the experimental values of A11 . Compared to SFTs, this is an unusually large value of CI . In Fig. 1(a), note that the A11 component is predicted with a reasonable level of accuracy, while in Fig. 1(b) the A22 and A33 components are underand over-predicted, respectively. The same trends were observed for other filling speeds, for both ISO plaque and disk geometries. More importantly, Nguyen et al. [17] showed that micromechanical predictions of material properties based on the predicted fiber orientation distributions do not adequately capture the mechanical anisotropy of the solid parts, while predictions using the measured fiber orientation distributions do. With a single, scalar CI parameter available to model the fiber–fiber interactions, the Folgar–Tucker model does a good job of predicting only one of the orientation components in our LFT moldings. We hypothesize that a model which can anisotropically model fiber–fiber interactions will be better able to predict all of the orientation components Aij . This should be especially true for LFT composites. In order to anisotropically model these interactions one must introduce a global tensorial representation of the interaction (rotary diffusion) term. Previous attempts at such a model will now be reviewed. 2.4. Previous ARD models

Fig. 1. Second-order orientation tensor components A11 (a), A22 , and A33 (b) for a 3 mm thick, slow-filled, long-glass-reinforced ISO plaque 90 mm long × 80 mm wide. The ORIENT predictions are performed with both the standard Folgar–Tucker model (CI = 0.03) and the RSC variant ( = 1/30 and CI = 0.03).

2.4.1. The Koch model Koch [12] considered the influence of long-range hydrodynamic fiber–fiber interactions on orientation in semi-dilute fiber suspensions. Koch developed a mechanistic model for orientational diffusivity, and derived an expression for a spatial tensor C such that C˙ represents rotary diffusion. Using the notation of our work, this tensor is given by C=

to the cross-flow and thickness directions, respectively. Fig. 1(a) and (b) show both the predicted and experimental values of A11 , A22 , and A33 versus the non-dimensional thickness coordinate z/h for a 3 mm thick, slow-filled, glass-reinforced ISO plaque at a position approximately halfway down the length of the channel. In these plots, we observe a region near the center of the molding where the A11 component has a low value. This is referred to as the core. The regions to either side of the core, with A11 values ranging from approximately 0.65 to 0.7, are referred to as the shell. Some moldings also show thin skin layers near the surface, where the fibers are less aligned in the flow direction than in the shell. In general, LFT composites exhibit wider cores and lower flow-direction (A11 ) alignment in the shell than their SFT counterparts. Fig. 1(a) and (b) also highlight the need for the RSC model. The standard Folgar–Tucker model predicts a very narrow core. However, by slowing the orientation kinetics with the RSC model, the core remains wide, consistent with the experimental data. While the need for this model is illustrated here with LFT composites, the same treatment is also needed with SFTs [4]. The value of CI is typically chosen to fit experimental data by considering the steady-state simple shear flow solution to the Folgar–Tucker or RSC model. Away from the mold inlet, the shellregion flow closely approximates a simple shear flow, and fibers in the shell have experienced very large shear strains during mold filling. Thus, we expect the orientation in the shell to match the

nL3 [ˇ1 (D : A : D)I + ˇ2 D : A : D], ˙ 2 ln2 ˛

(17)

where n is the number of fibers per unit volume, L is the fiber length, and ˛ = L/d is the fiber aspect ratio. While C is a spatial tensor, rotary diffusion only operates in orientation space, and Koch points out that “only the components of the tensor tangent to the unit sphere have physical significance.” The model assumes a semi-dilute suspension, and so is expected to be valid for nL2 d less than about 3. In terms of fiber volume fraction vf , this requirement is vf < 3/4˛. Koch obtains values of ˇ1 = 3.16 × 10−3 and ˇ2 = 1.13 × 10−1 by fitting the model to theoretical orientation distributions in extensional flows. In Eq. (17), A is the sixth-order orientation tensor, defined as



A=

pppppp dp.

(18)

˙ the sixth-order tenIn our work, where we seek an equation for A, sor must be approximated from the second-order tensor using a closure approximation. Our work utilizes the invariant-based INV6 closure approximation of Jack and Smith [19]. In order to explore the effectiveness of this model, the orientad tion tensor C may be incorporated into a model for A˙ , the diffusive d

˙ We develop a general expression for A˙ in Section contribution of A. 3.1 of this paper. For Koch’s model, this is d ˙ − 2(tr C)A − 5(C · A + A · C) + 10A : C]. A˙ = [2C

(19)

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

169

2.4.2. The Fan and Phan-Thien model d Phan-Thien et al. [11] propose replacing the A˙ term in Eq. (8) with d ˙ − 2 tr(C)A − 3(C · A + A · C) + 6A : C], A˙ = [2C

(20)

where C is a spatial tensor describing fiber–fiber interactions. The hydrodynamic term from Eq. (7) is unaffected. This model does sat˙ and Phan-Thien et al. [11] isfy the symmetry requirement on A, corrected an earlier problem from Fan et al. [10] so that the model maintains tr A = 1. However, our own preliminary calculations with Eq. (20) produced unrealistic dynamic behavior, and we discovered that this model suffers from a major flaw. When fibers are isotropically distributed, the second-order orientation tensor A takes on the exact value of Fig. 2. Dynamic solution for three cases of the Koch model in simple shear flow. L = 1000 ␮m and d = 17 ␮m.

Combining Eqs. (6), (7), (17) and (19) gives a complete description of orientation kinetics utilizing Koch’s anisotropic rotary diffusion model. Fig. 2 shows the orientation evolution of the A11 and A33 tensor components in simple shear flow (v1 = x ˙ 3 , where we use ˙ = 1) for the Koch model, for three examples of fiber volume fraction in simple shear flow. In each case, L = 1000 ␮m and d = 17 ␮m, giving ˛ = 58.8, which is representative of LFT fiber dimensions after injection molding [17]. vf = 0.19 corresponds to the volume fraction typically seen in LFT composites. One can compare the simple shear flow results at steady state to the experimental data shown in Fig. 1. In Fig. 1, the shell region (which is dominated by simple shear flow) of the LFTs exhibits experimental A11 and A33 components of approximately 0.65 and 0.01, respectively. In Fig. 2, however, the Koch model for vf = 0.19 gives A11 and A33 equal to approximately 0.43 and 0.23, respectively, at steady state. This result is overly diffusive, as the A11 and A33 components do not change significantly from their initial condition (A = (1/3)I). However, corresponds to nL2 d = 14.23, which is well outside of the semi-dilute regime for which Koch’s model was developed. If we reduce vf to 0.05, the Koch model does match the target of A11 = 0.65 at steady state, but it gives A33 = 0.10, which is almost identical to the Folgar–Tucker model. This suggests that the ˇ1 term in Eq. (17), which is an isotropic diffusivity whose magnitude varies with orientation state, is dominating the steady-state behavior compared to the anisotropic ˇ2 term. Fig. 2 also gives the Koch model orientation evolution for a fiber volume fraction of vf = 0.01 (nL2 d = 0.749). In this case, the A11 component rises above unity, while the A33 component falls below zero, results that are non-physical. The C tensor given by the Koch model will always have positive eigenvalues, which should guarantee a physically plausible solution for y(p,t). However, we have imple˙ which requires closure mented Koch’s model in an equation for A, approximations for both the fourth-order and sixth-order orientation tensors. The non-physical results at low vf are almost certainly caused by the closure approximations for A or A. These errors might be alleviated by using other closures, though the ORE closure for A and the INV6 closure for A are among the best available closures to date. Because of the need to calculate the sixth-order tensor, the Koch model is computationally much more expensive than the Folgar–Tucker model, yet in the best case it provides no better match to experimental data. These observations indicate that the Koch model provides no advantage for predicting fiber orientation in LFT composites.

A=

1 I, 3

(21)

and the components of the fourth-order orientation tensor A are given exactly by Aijkl =

1 (ı ı + ıik ıjl + ıil ıjk ), 15 ij kl

(22)

where ıij is the Kronecker delta. Introducing these tensors into Eq. (20) gives d 4 4 A˙ = C − (tr C)I. 5 15

(23)

d Thus, for an isotropic distribution of fibers, A˙ = / 0 for a general C. In other words, when the fibers are randomly distributed, the diffusion term of this model would pull the orientation away from isotropy. Diffusion provides a flux that is proportional to a gradient—in this case a gradient of . For an isotropic orientation state the gradient ∂ /∂p is zero, so diffusive flux must be zero, and d A˙ must be zero. The model of Phan-Thien et al. [11] does not have

this property, so it must be incorrect. Further examining the derivation of Phan-Thien et al. [11], we conclude that their projection of C onto the orientation space was incorrect. They correctly note that the rotary diffusion tensor C should be projected onto a surface normal to the orientation vector p. In projecting the tensor C, they left-multiply it by the tensor (I − pp). However, the correct projection involves both a left- and a right-multiplication by (I − pp), and the right-multiplication is omitted from their derivation. Rotary diffusion can neither translate a fiber nor change a fiber’s length. Thus, diffusion is defined only on the unit spherical surface traced by all possible fiber orientations, and not in the direction that the fiber is pointing. Rotary diffusion for fiber orientation is properly defined two-dimensionally in the orientation space. The derivation of Phan-Thien et al. [11] does not draw a clear distinction between spatial quantities (such as p and C) and surface quantities in the orientation space (such as s , q, or ). We correct this problem below by developing a model that incorporates a proper two-dimensional understanding of rotary diffusion. 3. Theory 3.1. Model development In order to properly introduce anisotropic rotary diffusion, we consider a two-dimensional rotary diffusivity surface tensor, denoted Dr and defined in surface spherical coordinates. The diffusive flux vector q for the ARD model is then q = −D ˙ r · ∇s ,

(24)

170

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

where s is the surface gradient operator defined in the orientation space. Here, we retain the assumption of Folgar and Tucker that q ∝ . ˙ Although diffusion is properly recognized as two-dimensional on the surface of a unit sphere, the orientation tensor whose evolution we are attempting to model is a space tensor, expressed in global Cartesian spatial coordinates. Some method is needed to pass diffusion information between these two coordinate systems. To do this, we define Dr as the projection of a three-dimensional space tensor B onto the local surface coordinates of the unit sphere. This projection of B is given by Dr = tT Bt,

(25)

where, using spherical coordinates for the orientation space, t is the 3 × 2 matrix, e.g. [20]:



t=

cos  cos cos  sin −sin 

−sin cos 0



.

(26)

The tensor B could presumably change with orientation p, so we write it as a function of a symmetric global tensor C and a tensorial representation of orientation, pp. Hand [21] gives the representation theorem for the ways that B can depend on C and pp. This is B = [a1 + a2 tr C + a3 tr C2 + p · (a4 C + a5 C2 ) · p]I + [a6 + a7 tr C + a8 p · C · p]C + a9 C2 .

(27)

Each ai could, most generally, be a scalar function of the joint invariants of C and pp. However, to keep the analysis tractable, we impose the constraint that B be no more than quadratic in C, which constrains each ai parameter to be a constant. Using Eqs. (3), (24), (25) and (27), we can write the rotary diffusion portion of the distribution function’s continuity equation as



d dt

d

= ˙ ∇ s · [(tT Bt) · ∇ s ].

(28)

d dt

d



pp dp

d = A˙ = ˙

 (29)

d

and . Thus, A˙ may be examined component-by-component, using the scalar equation:



pi pj ∇ s · [(tT Bt) · ∇ s ] dp.

(30)

Expanding the surface divergence via the product rule gives



A˙ dij = ˙

∇ s · [pi pj (tT Bt) · ∇ s ] dp

 − ˙

= −˙ = −˙ = ˙







(∇ s ) · [(tT Bt) · (∇ s pi pj )] dp

∇ s · [ (tT Bt) · ∇ s pi pj ]dp+˙



∇ s · [(tT Bt)·∇ s pi pj ] dp

∇ s · [(tT Bt) · ∇ s pi pj ] dp. (32)

The periodicity argument once again applies to the first integral in the second line, reducing it to zero. The surface gradient operator can be written in surface spherical coordinates as

∇s =

∂ 1 ∂ e + e , ∂  sin  ∂

(33)

where e and e are unit vectors in the local  and directions respectively. With t, B (via its dependence on p), and the components of p all trigonometric function of  and , Eq. (33) may be substituted into Eq. (32). The components of C which enter into the resulting equation may be factored outside the integral, leaving only trigonometric functions of  and inside the integral. The remaining trigonometric functions inside the integral can be rearranged so that, upon integrating over all p, only components of the A tensor remain. These multiply components of the C tensor that had previously been factored out of the integral. From that point, an extensive algebraic exercise yields each term of the rotary diffusion portion of the fiber orientation tensor evolution model. Reconstruction of the tensor representation gives d A˙ = (a ˙ 1 + a2 tr C + a3 tr C2 )(2I − 6A) + a4 [2(A ˙ : C)I

˙ : C2 )I + 2(C2 · A + A · C2 ) + 2(C · A + A · C) − 10A : C] + a5 [2(A − 10A : C2 ] + (a ˙ 6 + a7 tr C)[2C − 2(tr C)A − 5(C · A + A · C) ˙ tr (A · C) + 2(C2 · A + A · C2 ) + 10A : C] + a8 {2C

2 + 14C : A : C} + a9 [2C ˙ − 2(tr C2 )A − 5(C2 · A + A · C2 )

+ 10A : C2 ].

(34) h

{∇ s · [(tT Bt) · ∇ s ]}pp dp.

Here, note that the quantity inside the braces on the right-hand side is a scalar, the surface divergence of a surface flux vector q. This scalar multiplies the space tensor pp. Although pp is spatial, its components pi pj are scalar functions of the spherical coordinates 

A˙ dij = ˙

A˙ dij

− 2(tr C)(A : C) − 4A : C2 − 7[C · (A : C) + (A : C) · C]

Multiplying each side of this equation by the tensor pp, integrating over all orientations p, and factoring the time derivative out of the left-hand side yields:



The second integral may then be re-arranged and expanded such that

(∇ s pi pj ) · [(tT Bt) · ∇ s ] dp.

(31)

Note that (and its gradient), the components of p, and the surface diffusion tensor Dr = tT Bt are all periodic over the orientation space, so the first integral evaluates to zero. With B (and by extension Dr ) symmetric, we can re-order the dot product in the second integral.

For the entire orientation tensor evolution equation, A˙ = A˙ + d h A˙ , the hydrodynamic term A˙ remains unchanged from the Folgar–Tucker model. Note that setting a1 = CI and a2–9 = 0 recovers the Folgar–Tucker model, and that setting a6 = 1 and a1–5,7–9 = 0 recovers a corrected form of the Phan-Thien model, Eq. (20). Moreover, this model satisfies the symmetry requirements on A˙ (with d the assumption that C is symmetric) and also gives A˙ = 0 for a random distribution of fibers (A = (1/3)I), for all choices of a1–9 and all choices of C. Equation (34) is our most general form of the ARD model. However, owing to our earlier decision to constrain each ai parameter to be a scalar constant in order to preserve tractability in the derivation of the model, this is not the most general model. By treating the ai parameters as joint invariants of the C and pp tensors, a host of models may be developed using this derivation template. However, we refrain from developing any more terms for this model, and instead focus our attention on simplifying the model at hand to make it useful for applications. Preliminary numerical experiments indicated that not all of the ai terms from Eq. (34) were useful in predicting fiber orientation, and we chose to retain only the a6 term. While this is an arbitrary choice, we can argue for it heuristically as follows. First, note that the a1–3 terms retain the same form as the isotropic Folgar–Tucker

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

model. Although the C tensor can affect these terms as a contribution to their leading scalar factor, these terms remain expressions for isotropic rotary diffusion. From an application standpoint, these terms are equivalent to a Folgar–Tucker model with an adjustable CI parameter and can only reach the same steady orientation states as Folgar–Tucker. Thus, we do not further consider the a1–3 terms. Next, the remaining a5 and a7–9 terms are quadratic with respect to C, as opposed to the a4 and a6 terms, which are linear with respect to C. Solving for the components of C which best fit experimental orientation data is decidedly more difficult with these quadratic terms than with the linear terms. Thus, we set the a5 and a7–9 terms aside as a matter of convenience. Finally, by adjusting the remaining, linear, a4 and a6 parameters relative to one another, a host of models that exhibit anisotropic rotary diffusion and that are linear in the C tensor may be developed. However, the kernel of the a4 term, given by Eq. (27), is (p·C·p)I. Thus, the a4 term allows the local rate of diffusion to depend on orientation. However, at any single orientation p, rotary diffusion remains locally isotropic. The a6 term, though, exhibits true local anisotropy. By allowing C itself to be a function of orientation, we may then also introduce orientation dependent changes in the degree of isotropy. This choice is equivalent to setting B = a6 C in Eq. (27). We thus discard the a4 term, and focus attention on the a6 term. Thus far, we have imposed no restrictions on the values that the components of C may take. For the overall model to be objective, C must be a function of other objective tensors that appear in Eq. (34). We could create a quite general objective model by requiring C = C(D, A, A, A). The Koch model, Eq. (17), is an example of this. However, again for simplicity, we further assume that C = C(D, A). The representation theorem once again comes from Hand [21], where C = b1 I + b2 A + b3 A2 +

b4 b6 b5 D + 2 D2 + (D · A + A · D) ˙ ˙ ˙

+

b8 b7 (D · A2 + A2 · D) + 2 (D2 · A + A · D2 ) ˙ ˙

+

b9 2 2 (D · A + A2 · D2 ). ˙ 2

b5 b4 D + 2 D2 . ˙ ˙

Together, Eqs. (36) and (37) form what is referred to in the remainder of this work as the ARD fiber orientation model. This model’s equivalent kinetic theory model is d h = −∇ s · [ p˙ − (t ˙ T Ct) · ∇ s ], dt

(38)

where t is defined in Eq. (26). Just as the RSC model of Eq. (12) gives an objective means to slow the orientation kinetics of the Folgar–Tucker model of Eq. (6), one can also develop an ARD-RSC model which slows the kinetics of Eq. (37). Applying the RSC treatment embodied in Eqs. (10) and (11) to Eq. (37) gives

+ (1 − )(L − M : A)] : D} + {2[C ˙ − (1 − )M : C] − 2(trC)A − 5(C · A + A · C) (35)

(36)

As with the parameters of the general ARD model, our scalar constant constraint on the parameters of the definition of the C tensor simplifies the model for our applications. The second choice to set b6–9 equal to zero simplifies the process of fitting the bi parameters to experimental data. The choice to eliminate the terms with dependencies on products of A and D was pragmatic, and no theoretical consideration was given to including the other terms. The selection of the remaining b1–5 parameters will be covered in Secttion 3.2. Returning to our ARD fiber orientation model in which we have eliminated all but the a6 term, one may now set a6 equal to unity without loss of generality, since C remains adjustable via the b1–5 parameters. This simplifies our general ARD model to A˙ ARD = (W · A − A · W) + (D · A − A · D − 2A : D) + [2C ˙ − 2(tr C)A − 5(C · A + A · C) + 10A : C].

Fig. 3. Transient orientation kinetics in simple shear flow using the ARD model with bi = (2.0257, −26.9707, 43.0677, 1.1676 × 10−5 , −7.1671).

ARD,RSC A˙ = (W · A − A · W) + {D · A + A · D − 2[A

d This model ensures the objectivity of C, and also makes A˙ proportional to , ˙ consistent with the Folgar–Tucker model. Once again, each scalar bi could most generally be a function of the nine joint invariants of A and D. However, to keep the model tractable, we constrain the bi coefficients to be scalar constants and set b6–9 = 0, using:

C = b1 I + b2 A + b3 A2 +

171

(37)

+ 10[A + (1 − )(L − M : A)] : C}.

(39)

Note that setting  equal to 1 reduces the ARD-RSC model to the ARD model. The ARD-RSC model is the model whose properties we will explore. An effective technique for choosing the b1–5 parameters will now be introduced. 3.2. Model parameter selection Properly assigning values to the bi fitting parameters is a critical step in implementing the ARD model. The consequences of poor parameter selection are shown in Fig. 3. Here, the orientation approaches a reasonable-looking steady state and then falls into an oscillatory pattern, indicative of a dynamic instability in the model. Such behavior is clearly undesirable. The process of parameter selection first considers the ability of the model to produce the desired final orientation state in a molded part. In the injection molding process, the flow very near the upper and lower mold walls closely approximates a simple shear flow (SSF). Moreover, the strains accumulated there by the end of the filling are very large. Thus, the fiber orientation state close to the walls should be the steady state in simple shear flow. Phan-Thien et al. [11] choose the components of C in order to fit a desired steady-state orientation in simple shear flow. However, our C tensor in Eq. (36) is a function of A and D in order to retain ˙ Thus, we adopt a similar technique objectivity in our model for A.

172

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

Fig. 4. (a) Parameter search for ARD-RSC model in simple shear flow. (♦) A solution that both is stable and has positive steady-state eigenvalues for C, () a stable solution with negative eigenvalues in C, and () an unstable steady state. Here, the target SSF steady-state orientation matrix is given in Eq. (41). (b-d) Stable solutions with positive steady-state eigenvalues of C are successively analyzed for PEF, UEF, and BEF flows. Physically realizable solutions that give positive eigenvalues of C solutions are marked with a plus (+) symbol, and advance to the next flow under consideration. Physically realizable solutions that do not meet the eigenvalue of C requirement are marked with a mid-dot (·) symbol, and non-physical solutions are denoted with a times (×) symbol. (d) The valid solution located farthest from the boundary of invalid solutions is additionally marked with a oplus (⊕) symbol.

of fitting the bi parameters to a desired steady-state orientation tensor. Once the model has the desired steady-state behavior, we consider its dynamic properties. These include stability of the steady state, positive eigenvalues in the C tensor, and physically plausible transient solutions. Each of these criteria will now be discussed in detail. 3.2.1. Matching steady-state orientation data In a 1–3 simple shear flow, the velocity gradient tensor L is given by



LSSF =

0 0 0

0 1 0 0 0 0



.

(40)

For the LFT samples examined in our experiments, the orientation data near the walls is well represented by



ASSF =

0.65 0 0.03

0 0.03 0.34 0 0 0.01



.

(41)

With target values for the components of A and a given flow (i.e., a given L tensor), the only unknown quantities in the ARD or ARDRSC models are the bi parameters (and  in the case of the ARD-RSC ˙ If the targeted orientation tensor model) and the components of A. A is a steady-state solution, the value of the components of A˙ are all then zero, leaving a set of algebraic equations that is linear in bi .

For simple shear flow, the target orientation and velocity gradient tensors give A˙ 23 = A˙ 12 = 0 at all times. Also, with the requirement that tr A = 1 (and, by extension, tr A˙ = 0), one of the equations for the three diagonal components of A˙ is redundant. This leaves three independent equations, and five unknown model parameters. Thus, one can choose any two of the bi parameters independently, and then solve directly for the remaining three. We choose b3 and b5 as the independent parameters, and solve for b1 , b2 , and b4 to match the target steady-state orientation in simple shear flow. The necessary linear equations are derived in Appendix A. This method is guaranteed to yield a set of parameters for either the ARD or ARD-RSC model in which the target orientation tensor is a steady state of the chosen model in simple shear flow. However, depending on the values b3 and b5 , the model may or may not behave in other desirable ways, for either the simple shear flow or for other flows. 3.2.2. Stability Fig. 3 showed the oscillatory pattern into which an ARD model solution may settle when the fitting parameters are not carefully chosen. Note that, at around time t = 8, the model appears to be settling into a nearly steady solution, before beginning to oscillate. This is consistent with a saddle point instability. In order to eliminate instabilities such as those observed in Fig. 3, we perform a linear stability analysis of the model at the desired steady state. In a 1–3 simple shear flow, the independent components of the fiber orientation tensor are A11 , A33 , and A31 . We

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

173

numerically calculate the derivatives of the independent components of A˙ with respect to the independent components of A at the steady state, giving

⎡ ˙ ∂A11 ⎢ ∂A11 ⎢ ⎢ ∂A˙ 33 ∂A˙ =⎢ ∂A ⎢ ∂A11 ⎢ ⎣ ˙ ∂A31 ∂A11

∂A˙ 11 ∂A33

∂A˙ 11 ∂A31

∂A˙ 33 ∂A33

∂A˙ 33 ∂A31

∂A˙ 31 ∂A33

∂A˙ 31 ∂A31

⎤

⎥ ⎥ ⎥ ⎥. ⎥ ⎥ ⎦

(42)

For a stable steady state, the eigenvalues of this tensor must all be negative. 3.2.3. Eigenvalues of diffusivity tensor Next, one should consider that negative diffusivity in the 2D orientation space has no physical meaning. Thus, the two eigenvalues of Dr must be positive. Given the definition of Dr as a projection of B, ensuring that B (and by extension C) has positive eigenvalues will ensure that Dr will also have positive eigenvalues. This check on C is made at every time step during a dynamic orientation calculation for each flow of interest. 3.2.4. Validity in other flows Finally, the model’s behavior under other flow conditions must be assessed. Although there is usually no experimental data from other flows against which to compare orientation results, we want the solutions generated to be physically valid. The criteria for validity are: 0 ≤ A11 , A22 , A33 ≤ 1 and −0.5 ≤ A23 , A31 , A12 ≤ 0.5. We consider transient orientation behavior in planar elongational flow (PEF), uniaxial elongational flow (UEF), and biaxial elongational flow (BEF). Their respective velocity gradient tensors are



LPEF =

 L

UEF

=

and BEF

L

 =

1 0 0 0 −1 0 0 0 0 1 0 0

0 −0.5 0

0.5 0 0 0.5 0 0



,

(43)

0 0 −0.5

0 0 −1

 ,

(44)

 .

(45)

3.3. Searching for satisfactory model parameters Summarizing all of the criteria, a valid set of parameters bi must give a solution to either the ARD or ARD-RSC model in which • the steady, simple shear flow solution is equal to the target orientation tensor, • the steady, simple shear flow solution is linearly stable, • the C tensor has positive eigenvalues at all times in all flows, and • the transient solution is physically valid in simple shear, as well as planar, uniaxial, and biaxial elongation flows. Stability is not explicitly verified in the extensional flows. However, we have not observed any instabilities in any of the extensional flows for any parameter set. If the choice of the two independent bi parameters is done systematically, one can utilize the stability criteria and the positive eigenvalue of C criteria to develop a family of solutions which satisfies the first three described solution criteria for simple shear flow. Then, using these same parameters, a new velocity gradient tensor

Fig. 5. Dynamic solution to the ARD-RSC model in simple shear flow with  = 1/30, and the bi parameter set given in Eq. (46).

can be inserted into the model, and the validity of the model for the other flows may be assessed. As previously “good” solutions fail to meet these increasing criteria, they are discarded from further consideration. Eventually we find a family of parameter sets that satisfies all of the criteria above. We choose to systematically select b3 and b5 , and linearly solve for b1 , b2 , and b4 for each (b3 , b5 ) pair. The ability of each solution to meet the stability and positive eigenvalue of C requirements in simple shear is then recorded. It is convenient to select a grid of (b3 , b5 ) values, and display the ability of each solution to meet the SSF criteria in a plot in (b3 , b5 ) space. The ability of the parameter sets to meet the criteria of physically plausible solutions and positive eigenvalues of C in the extensional flows can be displayed in subsequent plots. Fig. 4(a) gives a criteria plot for the SSF solution given the target A of Eq. (41). Here we have utilized the ARD-RSC model with  = 1/30, as our initial mold-filling simulations with the RSC model indicated a significant slowing of the orientation kinetics in LFT materials. There are many solutions that meet the stability and C eigenvalue criteria for simple shear. The different elongational flows (Fig. 4(b) and (d)) gradually eliminate some of these b3 –b5 combinations, but do leave multiple solutions that meet all of the previously described criteria for all of the flows. To select a single parameter set from the remaining valid solutions (Fig. 4(d)), we implement a boundary erosion technique common in image processing to isolate the point furthest from the boundary of the region of valid solutions. This point is also indicated in Fig. 4(d). Fig. 5 gives the orientation kinetics in simple shear flow for the chosen parameter set, which is bi =(1.924 × 10−4 , 5.839 × 10−3 , 4.0 × 10−2 , 1.168 × 10−5 , 0). (46) The initial condition in Fig. 5 is A = (1/3)I, a random arrangement of fibers. The evolution of the A23 and A12 components are not reported, as they remain at their initial zero values. The orientation kinetics of all of the valid bi solution sets are quite similar, so Fig. 5 is representative of this entire family of solutions for the ARD-RSC model with  = 1/30 and the target steady state of Eq. (41). Suitable parameter sets for other target orientation states can be found using the same procedure. This ARD-RSC model gives the type of behavior we desire for incorporation into a mold-filling simulation for LFT materials.

174

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

5. Conclusion

Fig. 6. Second-order orientation tensor components A11 (a), A22 , and A33 (b) for a 3 mm thick, slow-filled, glass-reinforced ISO plaque 90 mm long × 80 mm wide. The ORIENT predictions use the RSC variant of the Folgar–Tucker model with  = 1/30 and CI = 0.03, while the ORIENT-ARD predictions use the ARD-RSC model with  = 1/30 and the bi parameter set given in Eq. (46).

4. Application to LFT injection molding The ARD model family was implemented in ORIENT to predict the development of orientation during mold filling in LFT composites. Using the ARD-RSC model with  = 1/30 and the bi parameter set given in Eq. (46), we compared the predicted orientation tensor components to the experimental components of the same specimen and the same location that was examined in Fig. 1. This comparison is shown in Fig. 6. In this comparison, the A11 data is reasonably well predicted by both the RSC and ARD-RSC models, as seen in Fig. 6(a). However, using the new ARD-RSC fiber orientation model (Fig. 6(b)), the A22 and A33 components are much better predicted than before (Fig. 1(b)). There is still a poor fit near z/h = 0 in both the A11 and A22 comparisons, which seems due to asymmetry in the experimental data (ORIENT assumes symmetry about the midplane z = 0). This asymmetry could come from mold wall temperature differences, resulting in asymmetric cooling of the polymer melt, or from the fact the gate in this cavity is not located symmetrically about that cavity midplane. Finally, using the orientation distribution predicted by the ARDRSC model to predict elastic moduli in the flow and crossflow directions gives results that compare favorably with the measured moduli [17,22]. This is an improvement compared to the agreement observed when the RSC model was used to predict the fiber orientation distribution.

Adding anisotropic rotary diffusion (ARD) to a Jeffery-type fiber orientation model offers greatly improved ability to match experimental data, compared to the isotropic Folgar–Tucker model. In this paper we have derived the moment-tensor equation for an anisotropic rotary diffusion that is described by a secondorder space tensor, correcting an error in the equation given by Phan-Thien et al. [11]. Combining this equation with a closure approximation for the sixth-order orientation tensor allowed us to explore the model of Koch [12]. That model was found to offer no improvement over the Folgar–Tucker form, at least for modeling long-fiber thermoplastic composites. Moreover, at low fiber volume fractions our implementation of Koch’s model gave non-physical results, probably induced by the closure approximation. Proceeding with a phenomenological model, we allowed the rotary diffusion tensor to depend on fiber orientation state and on the rate-of-deformation tensor. The resulting model has five scalar parameters: three are chosen to match a target steady-state orientation in simple shear flow, and the other two are chosen to ensure a stable steady state, positive eigenvalues in the diffusion tensor, and physically valid results in various elongational flows. We show that the new model provides improved predictions of fiber orientation in injection molded long-fiber thermoplastic composites. While long-fiber thermoplastics provided the practical motivation for this work, our anisotropic rotary diffusion model may have other applications as well. Fan et al. [10] and Phan-Thien et al. [11] found that their direct numerical simulation results for semi-concentrated fiber suspensions could not be matched by the Folgar–Tucker model. Our ARD model can easily be fitted to this type of result using the parameter-selection procedure shown here. For short-fiber reinforced thermoplastics, the Folgar–Tucker model can match experimental rheology data for the shear viscosity, but in some cases cannot simultaneously match the normal stress data [6]. Because the ARD model offers the ability to tailor all components of the orientation tensor at steady state, it should provide a better match to rheological data for fiber-reinforced materials. Our model for the rotary diffusion tensor, Eq. (36), is purely phenomenological. It would be very helpful to have a mechanistic model that would give the rotary diffusion as a function of local orientation, fiber orientation state, and type of deformation. When such models are developed, the framework of this paper will allow those models to be implemented in moment-tensor form, so that their properties can be explored and their usefulness determined. Acknowledgements The work presented here was supported by the Department of Energy through the Pacific Northwest National Laboratory. Comments and suggestions from Drs. Mark T. Smith, Ba Nghiep Nguyen, and James D. Holbery have been most helpful. The experimental orientation measurements made by Vlastimil Kunc of the Oak Ridge National Laboratory were integral to this work, and we are grateful for his input. Appendix A. Derivation of linear equations for bi parameters Consider the ARD-RSC model, Eq. (39), and the definition of the global, spatial rotary diffusion tensor C in Eq. (36). Combining these two equations gives ARD,RSC A˙ = (W · A − A · W) + (D · A + A · D − 2[A + (1 − )(L

˙ − (1 − )M : I] − 2(tr I)A − M : A)] : D) + b1 (2[I

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

− 5(I · A + A · I) + 10[A + (1 − )(L − M : A)] : I) ˙ − (1 − )M : A] − 2(tr A)A − 5(A · A · A · A) + b2 (2[A

(51) by choosing the rows corresponding to A˙ 11 , A˙ 33 , and A˙ 31 . We write this as

2 + 10[A + (1 − )(L − M : A)] : A) + b3 (2[A ˙

˜ [B]

− (1 − )M : A2 ] − 2(tr A2 )A − 5(A2 · A + A · A2 ) + 10[A + (1 − )(L − M : A)] : A2 )

b5 + 10[A + (1 − )(L − M : A)] : D) + (2[D2 ˙

= (W · A − A · W) + (D · A + A · D − 2[A

2 ˙ : A − A2 ) + 2b3 (A ˙ − (tr A2 )A + 10b2 (A



ˆ = {H},

(53)

ˆ is a 3 × 3 matrix, and {H} ˆ is the 3 × 1 column vector with where [B] components: ˆi = H ˜ i − B˜ i3 b3 − B˜ i5 b5 , H

(54)

(55)

+ 5A : A − 5A ) + b4 (2[D − (1 − )M : D]

˙ − 3A33 ), Bˆ 21 = 2(1

(56)

− 5(D · A + A · D) + 10[A + (1 − )(L − M : A)] : D)

˙ Bˆ 31 = −6A 31 ,

(57)

˙ Bˆ 12 = 10(A 11ij Aij − A1i Ai1 ),

(58)

˙ Bˆ 22 = 10(A 33ij Aij − A3i Ai3 ),

(59)

˙ Bˆ 32 = 10(A 31ij Aij − A3i Ai1 ),

(60)

+

3

b5 (2[D2 − (1 − )M : D2 ] − 2(tr D2 )A − 5(D2 · A ˙

+ A · D2 ) + 10[A + (1 − )(L − M : A)] : D2 ).

(48)

We now utilize a notation in which the components of a symmetric ARD,RSC tensor are written as a column vector. For A˙ (now denoted A˙ for brevity), this is

⎧ ⎫ A˙ 11 ⎪ ⎪ ⎪ ⎪ A˙ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ˙ 22 ⎪ ⎬

A33 . A˙ 23 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ˙ ⎪ ⎪ ⎪ ⎩ A˙ 31 ⎪ ⎭ A12

(49)

⎪ ⎪ ⎪ ⎩ b4 ⎪ ⎭

Bˆ 23 = 2[D33 − (1 − )M33ij Dij ] − 5(D3i Ai3 + A3i Di3 ) + 10[A33ij Dij + (1 − )(L33ij Dij − M33ij Aijkl Dkl )],

+ 10[A31ij Dij + (1 − )(L31ij Dij − M31ij Aijkl Dkl )], (50)

we may combine Eqs. (48) and (50), adopt the notation of Eq. (49), and write the ARD-RSC model as

⎧ ⎫ b1 ⎪ ⎪ ⎪ ⎨ b2 ⎪ ⎬

(62)

(51)

b5

Here [B] is a 6 × 5 matrix whose components Bij are the coefficients of the parameters bi in Eq. (48), and {H} is a 6 × 1 column vector. ˙ A, D, W, and  (since A, L, The components Bij and Hi depend on A, and M are functions of A). For a steady-state orientation, A˙ = 0. Then, given a simple shear flow and a target steady-state value of A, all of the Bij and Hi components are determined. However, due to the property that tr A = 1 and to the symmetries of the simple shear flow (SSF), there are only three independent scalar equations in Eq. (51). Thus, we modify Eq.

(63)

ˆ 1 = A˙ 11 − (W1i Ai1 − A1i Wi1 ) − {D1i Ai1 + A1i Di1 − 2[A11ij Dij H + (1 − )(L11ij Dij − M11ij Aijkl Dkl )]} − 2b3 (A ˙ 1i Ai1 − A11 Aij Aji + 5A11ij Aik Akj − 5A1i Aij Aj1 ) −

= {H}.

(61)

Bˆ 33 = 2[D31 − (1 − )M31ij Dij ] − 5(D3i Ai1 + A3i Di1 )

H = A˙ − (W · A − A · W) − (D · A + A · D − 2[A + (1 − )(L − M : A)] : D),

Bˆ 13 = 2[D11 − (1 − )M11ij Dij ] − 5(D1i Ai1 + A1i Di1 ) + 10[A11ij Dij + (1 − )(L11ij Dij − M11ij Aijkl Dkl )],

Defining H as

b3

b1 b2 b4

˙ − 3A11 ), Bˆ 11 = 2(1

2

[B]

(52)

where i ranges from 1 to 3. ˆ Following through with the algebra, the components of the [B] ˆ vector are matrix and the {H}

+ (1 − )(L − M : A)] : D) + 2b1 (I ˙ − 3A)

˙ = {A}

˜ = {H},

˜ is a 3 × 5 matrix, and {H} ˜ is a 3 × 1 column vector. where [B] Now, five bi parameters must satisfy three equations. We choose b3 and b5 as independent variables, and modify Eq. (52) such that ˆ [B]

(47)

which is linear in bi . Noting that M : I = I, L : I = A, A : I = A, M : A = A, L : A = A2 , tr A = 1, tr I = 3, and tr D = 0, the model can be simplified, leaving A˙

b3

⎪ ⎪ ⎪ ⎩ b4 ⎪ ⎭



− (1 − )M : D2 ] − 2(tr D2 )A − 5(D2 · A + A · D2 )

ARD,RSC

⎧ ⎫ b1 ⎪ ⎪ ⎪ ⎨ b2 ⎪ ⎬ b5

+ b4 (2[D − (1 − )M : D] − 2(tr D)A − 5(D · A + A · D)

+ 10[A + (1 − )(L − M : A)] : D2 ),

175

b5 {2[D1i Di1 − (1 − )M11ij Dik Dkj ] ˙

− 2A11 Dij Dji − 5(D1i Dij Aj1 + A1i Dij Dj1 ) + 10[A11ij Dik Dkj + (1 − )(L11ij Dik Dkj − M11ij Aijkl Dkm Dlm )]},

(64)

ˆ 2 = A˙ 33 − (W3i Ai3 − A3i Wi3 ) − {D3i Ai3 + A3i Di3 − 2[A33ij Dij H + (1 − )(L33ij Dij − M33ij Aijkl Dkl )]} − 2b3 (A ˙ 3i Ai3 − A33 Aij Aji + 5A33ij Aik Akj − 5A3i Aij Aj3 ) −

b5 {2[D3i Di3 − (1 − )M33ij Dik Dkj ] ˙

− 2A33 Dij Dji − 5(D3i Dij Aj3 + A3i Dij Dj3 ) + 10[A33ij Dik Dkj + (1 − )(L33ij Dik Dkj − M33ij Aijkl Dkm Dlm )]},

(65)

176

J.H. Phelps, C.L. Tucker III / J. Non-Newtonian Fluid Mech. 156 (2009) 165–176

ˆ 3 = A˙ 31 − (W3i Ai1 − A3i Wi1 ) − {D3i Ai1 + A3i Di1 − 2[A31ij Dij H + (1 − )(L31ij Dij − M31ij Aijkl Dkl )]} − 2b3 (A ˙ 3i Ai1 − A31 Aij Aji + 5A31ij Aik Akj − 5A3i Aij Aj1 ) −

b5 {2[D3i Di1 − (1 − )M31ij Dik Dkj ] ˙

− 2A31 Dij Dji − 5(D3i Dij Aj1 + A3i Dij Dj1 ) + 10[A31ij Dik Dkj + (1 − )L31ij Dik Dkj − M31ij Aijkl Dkm Dlm )]},

(66)

Eq. (53) may now be solved to find b1 , b2 , and b4 for any choice of b3 and b5 . References [1] F. Folgar, C.L. Tucker III, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plast. Compos. 3 (1984) 98–119. [2] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. A 102 (1922) 161–179. [3] S.G. Advani, C.L. Tucker III, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 31 (1987) 751–784. [4] H.M. Huynh, Improved fiber orientation predictions for injection-molded composites, Master’s Thesis, University of Illinois at Urbana-Champaign, 2001. [5] M. Sepehr, G. Ausias, P.J. Carreau, Rheological properties of short fiber filled polypropylene in transient shear flow, J. Non-Newtonian Fluid Mech. 123 (2004) 19–32. [6] J. Wang, J.F. O’Gara, C.L. Tucker III, An objective model for slowing orientation kinetics in concentrated fiber suspensions: theory and rheological evidence, J. Rheol. 52 (5) (2008) 1179–1200. See also US Patent 7,266, 469 (2007). [7] R.S. Bay, C.L. Tucker III, Fiber orientation in simple injection moldings. Part I. Theory and numerical methods, Polym. Compos. 13 (4) (1992) 317–331.

[8] R.S. Bay, C.L. Tucker III, Fiber orientation in simple injection moldings. Part II. Experimental results, Polym. Compos. 13 (4) (1992) 332–342. [9] S. Ranganathan, S.G. Advani, Fiber–fiber interaction in homogeneous flows of nondilute suspensions, J. Rheol. 35 (1991) 1499–1522. [10] X. Fan, N. Phan-Thien, R. Zheng, A direct simulation of fibre suspensions, J. Non-Newtonian Fluid Mech. 74 (1998) 113–135. [11] N. Phan-Thien, X. Fan, R.I. Tanner, R. Zheng, Folgar–Tucker constant for a fibre suspension in a Newtonian fluid, J. Non-Newtonian Fluid Mech. 103 (2002) 251–260. [12] D.L. Koch, A model for orientational diffusion in fiber suspensions, Phys. Fluids 7 (1995) 2086–2088. [13] J.S. Cintra Jr., C.L. Tucker III, Orthotropic closure approximations for flowinduced fiber orientation, J. Rheol. 39 (1995) 1095–1122. [14] B.E. VerWeyst, Numerical predictions of flow-induced fiber orientation in 3-D geometries, Ph.D. Thesis, University of Illinois at Urbana-Champaign, Urbana, IL, 1998. [15] C.A. Hieber, S.F. Shen, A finite-element/finite-difference simulation of the injection-molding filling process, J. Non-Newtonian Fluid Mech. 7 (1980) 1–32. [16] International Organization for Standardization, Plastics—injection moulding of test specimens for thermoplastic materials. Part 5. Preparation of standard specimens for investigation anisotropy, ISO 294-5:2001 (2001). [17] B.N. Nguyen, Pacific Northwest National Laboratory, Personal communication (2008). [18] P.J. Hine, N. Davidson, R.A. Duckett, A.R. Clarke, I.M. Ward, Hydrostatically extruded glass–fiber reinforced polyoxymethylene. I. The development of fiber and matrix orientation, Polym. Compos. 17 (1996) 720–729. [19] D.A. Jack, D.E. Smith, An invariant based fitted closure of the sixth-order orientation tensor for modeling short-fiber suspensions, J. Rheol. 49 (2005) 1091–1115. [20] R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, PrenticeHall, Englewood Cliffs, NJ, 1962. [21] G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. 13 (1962) 33–64. [22] B.N. Nguyen, V. Kunc, B. Frame, J.H. Phelps, C.L. Tucker III, S.K. Bapanapalli, J.D. Holbery, M.T. Smith, Fiber length and orientation in long-fiber injection-molded thermoplastics. Part I. Modeling of microstructure and elactic properties, J. Compos. Mater. 42 (2008) 1003–1029.