A constitutive equation for fiber suspensions in viscoelastic media

Understanding the fibers suspended in a polymer melt is an important challenge in rheology with viscoelastic properties strongly depending on a fiber orientation. From experimental observations and theoretical analysis, higher fiber alignment leads to lower viscosity compared to random orientation states.^1,2 During the startup shear flows with small shear rates, the change of fiber orientation induces a viscosity overshoot for fiber-filled fluids, whereas no such finding occurs in unfilled fluids.^3–6 Figure 1 illustrates the normalized transient shear viscosity with respect to shear strains for the unfilled and fiber-filled fluids. This implies that the rheological nonlinearity of a viscosity overshoot is relatively increased, when fibers are filled to a neat polymer.

The pioneering work regarding the constitutive equation in suspension mechanics originated from Hinch and Leal.^7,8 Later, Dinh and Armstrong^9,10 developed the typical constitutive equation of fiber orientation-induced extra stress in a Newtonian viscous matrix, as expressed below

where $τ$ is the stress tensor, $ηm$ is the matrix viscosity in the absence of fibers, D is the rate-of-deformation tensor, $A4$ is the fourth order orientation tensor, $D:A4$ is the term for flow-fiber coupling, and $Np$ is the particle number characterizing the fiber-flow coupling or anisotropic degree. For given $Np=0$⁠, the suspensions recover the standard Newtonian viscous fluid. Some studies^11,12 of Tucker, Phan–Thien, and Graham indicated that the number $Np$ is directly proportional to the square of fiber aspect ratio a_r and the fiber volume fraction $ϕ$⁠; namely, $Np∝ar2ϕ$⁠. Recently, Tseng et al.¹³ proposed $Np$ as a function of $ar2ϕ$ and shear rate $γ̇$ for anisotropic flow simulations in the injection molding of fiber-reinforced thermoplastics. In particular, Pipes and coworkers^14–16 constructed the non-Newtonian constitutive relationships for hyper-concentrated fiber suspensions with an oriented fiber assembly.

Significantly, Sepehr et al.⁴ and Eberle et al.⁶ have used the Dinh–Armstrong constitutive equation to predict the transient shear viscosity; the changes of fiber orientation were determined by the famous Folgar–Tucker model^17–20 attached with a slip (or strain reduction) factor of slowing down the orientation response.^20–23 Note that in using such a slip factor, the fiber orientation predictions do not possess rheological objectivity; namely, different results^4,24 were derived in different coordinate systems.^20–22 However, those results^4,24 indicated that the predicted viscosity overshoot did not match the related experimental data. Recently, Favaloro et al.^25–27 and Tseng and Favaloro²⁸ derived the anisotropic fourth-order viscosity tensor to obtain a so-called informed isotropic (IISO) viscosity scalar depending on the flow field and the fiber orientation. Under a shear flow, the analytical shear stress of the IISO model is similar to that of the Dinh–Armstrong model.²⁸ As expected, the IISO model will not provide satisfactory viscosity overshoot predictions.

Furthermore, Férec et al.^29–32 attempted to extend the Dinh–Armstrong model in nonlinear Newtonian viscous fluids, including the Power-law model,³⁰ the Carreau model,²⁹ and the Bingham model.³² In addition, they³¹ further proposed the polymer matrix based on a second-order fluid with normal stress differences for fiber suspensions. Hence, their predictions on the viscosity overshoot were qualitatively verified with respect to model parameters. Moreover, Ait-Kadi and Grmela³³ particularly developed a fiber suspension rheological model in the FENE-S (finitely extensible nonlinear elastic spring) viscoelastic media. Azaiez³⁴ derived the FENE and Giesekus viscoelastic models incorporating the fiber-induced extra stress. Following Azaiez's work,³⁴ Ramazani et al.^35,36 performed that the model predictions of transient shear stress were only in qualitative agreement with the experimental data. However, the viscosity overshoot predicted by the aforementioned fiber suspension models for viscoelastic media remains undetermined.

Most recently, Tseng³⁷ developed a weighted shear/extensional viscosity for the generalized-Newtonian-fluid (GNF), known as GNF-X (eXtended GNF). Based on GNF-X, Tseng³⁸ continued to extend the classical White–Metzner (WM) constitutive equation of viscoelastic fluids, called WMT-X (WM eXtended by Tseng). Thereby, the WMT-X model prediction of transient shear viscosity with overshoot³⁹ matched the experimental data and the K-BKZ (Kaye-Bernstein–Kearsley–Zapas) integral constitutive model prediction of Ebrahimi et al.⁴⁰ for a high-density polyethylene (HDPE) melt. Significantly, the contraction flow patterns were exhibited, and the formation and growth of the corner vortex were explored, which are related to the extensional viscosity and the Weissenberg number. Therefore, the WMT-X model coupled with the Dinh–Armstrong flow-fiber coupling term for fiber suspensions in a viscoelastic media can be proposed, as expressed below

where $λ$ is the relaxation time, $ηW$ is the weighted viscosity of the GNF-X model, and $λ$ and $ηW$ are functions of strain rates $γ̇$⁠; this is called the suspended WMT-X model. When $λ$ = 0, such a model can return to the Dinh–Armstrong fiber suspension model.

Note that $τ□$ is called the Gordon–Schowalter hybrid time derivative,⁴¹ with a slip factor $ξ$⁠, is a linear combination of the Oldroyd upper-convected time derivative $τ∇$ and the Jaumann corotational time derivative⁴² $τ○$

where W is the vorticity tensor, which is the anti-symmetric matrix of the velocity gradient $∇v$⁠, and $DτDt$ is the well-known material derivative convective term. In particular, two conditions, $∂τ∂t$ = 0 and $∇τ$ = 0, indicate the steady state flow and the homogenous flow, respectively; the slip factor $ξ$ = 0 and 1 signify the upper-convected and corotational forms of the WM model, respectively. Additionally, $ξ$ was further interpreted through the framework of the Oldroyd eight-constant model of Saengow et al.^43,44

The GNF-X weighted viscosity $ηW$ of Tseng³⁷ is defined in the coexistence of shear and extensional flows, as expressed below

where W is the weighting function, also called the extension fraction, $ηS$ and $ηE$ are the GNF shear and extensional viscosities with respect to strain rates, respectively, $γ̇S$ and $γ̇E$ are the principal shear rate and principal extensional rate, respectively, and $γ̇T$ is the total strain rate. For a pure shear flow, the GNF-X weighted viscosity can return to the GNF shear viscosity, namely, $ηW=ηS$⁠. The details of the GNF-X weighted viscosity are available elsewhere.³⁷ Additionally, Park⁴⁵ was interested in the GNF-X model in seeking a critical issue: the weighted viscosity $ηW$ is nonobjective for a partitioning of the velocity gradient in a pure extension flow.

The Weissenberg number $Wi$ is the product of the shear rate $γ̇$ and the longest relaxation time $λ$

At a large Weissenberg number, $Wi≫1$⁠, the fluid will respond with the elastic behavior like that of a solid. For $Wi≪1$⁠, a liquid-like viscous state is expected. In a homogenous shear flow, the suspended WMT-X model is further expanded

where D* and W* are dimensionless tensors of the rate-of-deformation tensor D and the vorticity tensor W, respectively. In short, the suspended WMT-X model essentially contains five parameters: $λ$⁠, $Wi$⁠, $ξ$⁠, $ηS$⁠, and $Np$⁠.

In the orientation kinetics of fiber suspensions, the second-order orientation tensor A is the symmetric matrix defined by Advani and Tucker⁴⁶ 

where its trace is A₁₁ + A₂₂ + A₃₃ = 1. Physically, A = I/3 represents the isotropic orientation state, wherein I is the identity matrix. Three diagonal components: A₁₁, A₂₂, and A₃₃, correspond to the directions of flow, crossflow, and thickness, respectively. The off-diagonal components, A₁₂, A₁₃, and A₂₃, of a fiber orientation tensor represent that alignments vary from the coordinate axes. In regard to the Dinh–Armstrong flow-fiber coupling term D : A₄, the fourth order orientation tensor A₄ calculation decoupled in terms of the second-order orientation tensor A is obtained by higher order polynomial closure approximations such as the eigenvalue-based optimal fitting (EBOF) closure^47,48 and the invariant-based optimal fitting (IBOF) closure.⁴⁹ 

A time evolution equation of the second-order orientation tensor A is fixed on the material derivative, denoted as $Ȧ$⁠. At present, the famous fiber orientation models include the Folgar–Tucker IRD (isotropic rotary diffusion) model,¹⁸ the Phelps–Tucker ARD (anisotropic rotary diffusion) model,⁵⁰ the Wang–Tucker RSC (reduced strain closure) model,²⁰ and the ARD-RSC model.⁵⁰ Following the fiber orientation models mentioned above, Tseng et al.^22,23 developed the objective fiber orientation model, known as iARD-RPR (improved ARD model and retarding principal rate model) with just three physically related parameters. It is significant that Favaloro and Tucker²⁷ have reviewed and recommended existing fiber orientation models.

In the present work, the iARD-RPR model is therefore adopted to predict the time evolution of a fiber orientation tensor, which contains three $Ȧ$ terms

First, $ȦHD$ is the Jeffery hydrodynamic (HD)

where $ξf=(ar2−1)/(ar2+1)$ is a shape factor of a particle, a fiber's aspect ratio $ar$ is the ratio of its length $Lf$ to its diameter $Df$⁠, and $ar=Lf/Df$⁠. In general, $ξf$ is close to one for a fiber, $ξf≈1$⁠.

Second, the $ȦiARD$ term describes the anisotropic rotary diffusion of fiber suspensions in a flow field

where $Dr$ is the anisotropic rotary diffusion tensor, C_I is the Folgar–Tucker isotropic parameter due to the fiber–fiber interaction for a short fiber, C_M is the anisotropic parameter due to the fiber–matrix interaction for a long fiber, and the scalar $‖D2‖=12D2:D2$ is the norm of the tensor $D2$⁠. Note that C_M = 0 is suggested for a short fiber.

The later contributor of the RPR model to slow down the rate of an orientation tensor is defined as

where $Λ̇IOK$ is the material derivative of a particular diagonal tensor and its superscript indicates the intrinsic orientation kinetics (IOK) assumption,^22,$R$ is the rotation matrix, $RT$ is the transpose of $R$⁠, the superscript T is the transpose operator of a matrix throughout this paper, $λi$ is the eigenvalues of $A$⁠, $λ1≥λ2≥λ3$⁠, and $R=[e1,e2,e3]$ is defined by eigenvector columns of $A$⁠. The parameter $α$ is in dimensionless units, and the subscripts i, j, and k are indices of permutation. In short, the iARD-RPR model contains three parameters: $CI$⁠, $CM$⁠, and $α$⁠; the details are available elsewhere.²² When $CM$ = 0 and $α$ = 0, the iARD-RPR can return to the standard Folgar–Tucker (FT) model. In addition, Tseng et al.⁵¹ presented valuable reports for dramatic changes in the skin–shell–core structure of the fiber orientation distribution with respect to the model parameters. Those fibers found in the shell region are mainly aligned in the flow direction due to high shear rates near the mold surface, but others at the core are more transversely aligned to the flow due to low shear rates and extensional flows. The fibers in the skin immediately adjacent to the wall of the cavity show a random in-plane orientation as a result of a fountain flow.⁵² 

For a homogenous simple shear flow with an x-axis flow direction, a y-axis in a gradient direction, and a z-axis in a neutral direction, the flow strength of velocity-gradient tensor $∇v$ is given by shear rate $γ̇$ to determine the rate-of-deformation tensor D and the vorticity tensor W, which are the symmetric matrix and the anti-symmetric matrix of the velocity gradient tensor, respectively,

These tensors are input to the WMT-X of Eq. (12) and the iARD-RPR model of Eq. (16). The ordinary differential equations of the WMT-X and iARD-RPR models were solved numerically by the fourth order Runge–Kutta method. The initial condition of an isotropic orientation state, A = I/3, was given. One can first solve the iARD-RPR orientation equation to obtain the second-order orientation tensor A. The fourth order orientation tensor A₄ was approximated by the IBOF closure⁴⁹ based on the predicted second-order orientation tensor A. Thus, the Dinh–Armstrong flow-fiber coupling term D : A₄ was estimated in the suspended WMT-X calculation. For simplification in the present work, the parameters were constrained in the upper-convected time derivative and the short fiber orientation; namely, $ξ$ = 0 and C_M = 0, respectively.

One can examine two fixed fiber orientation states: ideal random and perfect alignment, to analyze the suspended WMT-X model predictions of the transient shear viscosity with respect to shear strains at a shear rate of $γ̇$ = 1.0 s⁻¹. The model parameters are given: $ηS$ = 1000 Pa·s, $λ$ = 1.0 s, $Wi$ = 1.0, N_p = 3.0, C_I = 0.01, and $α$ = 0.8. As shown in Fig. 2, the convergent viscosity of the completely aligned suspension is lower than that of the randomly oriented suspension. This result is in good agreement with experimental observation and model predictions of Powell and coworkers.^1,2 In addition, the time evolution of an orientation tensor (or varied alignment) was performed from the initial random condition to converge into a partial alignment. The “viscosity overshoot” is obviously found. The convergent viscosity of the partially oriented suspension falls between the randomly and completely aligned suspensions. Thus, the changes in the diagonal (A_xx, A_yy, A_zz) and off diagonal orientation tensor components (A_xy) with respect to shear strains, as shown in Fig. 3, can be further investigated. Consequently, while such an overshoot depends on the off-diagonal orientation tensor component in the flow-gradient plan, it is not related to the diagonal orientation components. For the Dinh–Armstrong fiber suspension model, the analytical shear stress^{4,20,24,26,28} is

Thus, one can verify that the change of fourth order orientation tensor component A_xyxy surely dominates the transient shear stress incorporating the overshoot.

Referring to the previous study of Eberle et al.⁶ regarding fiber suspension rheology, the experimental data and the model prediction for the fiber orientation were carried out under the transient shear flow for a 30 wt. % SGF/PBT (polybutylene terephthalate resin containing 30% by weight of short glass fibers) composite at $γ̇$ = 1.0 s⁻¹ and 260 °C. Figure 4 shows the time evolution of the diagonal orientation tensor components with respect to shear strains. In the present work, the iARD-RPR model prediction can better fit the experimental data than the previous work of Eberle et al.,⁶ which used the Folgar–Tucker model attached with a slip factor. The optimal iARD-RPR parameters are addressed in Table I. Note that the parameter C_M was set prior to zero for reducing the complexity; thus, the iARD-RPR model can return to the FT-RPR model. The accuracy of the fiber orientation prediction is positively verified.

Eventually, Fig. 5 shows the transient shear viscosity at different shear rates of $γ̇$ = 1.0 and 6.0 s⁻¹. The model parameters are given in Table II. The suspended WMT-X model predictions are in good agreement with experimental data. Especially for the viscosity overshoot at $γ̇$ = 6.0 s⁻¹, this result enhanced the numerical prediction of Eberle et al.⁶ using the Dinh–Armstrong model in a viscous media. When the relaxation time $λ$ = 0 is given, the suspended WMT-X model returns to the Dinh–Armstrong model. Hence, the Dinh–Armstrong prediction's overshoot is relatively small. In summary, the ultimate goal of this work was accomplished for rheological models of fiber suspensions, namely, proposing the suspended WMT-X model to fundamentally explore the effect of the fiber orientation on viscoelastic behaviors. It is well known that fourth order tensor-based viscosity models have seen much difficulty in numerical applications. Based on the fourth order viscosity tensor, the IISO viscosity model is more useful in practical injection and compression molding simulations of fiber-reinforced thermoplastics.^25,28 Similar to the IISO model applications, peculiar complex flows induced by the anisotropic effect of the fiber orientation will be investigated via the suspended WMT-X viscoelastic fluid model in future work.