Abstract

Translate

Fiber orientation is an important descriptor of the microstructure for short fiber polymer composite materials where accurate and efficient prediction of the orientation state is crucial when evaluating the bulk thermo-mechanical response of the material. Macroscopic fiber orientation models employ the moment-tensor form in representing the fiber orientation state, and they all require a closure approximation for the higher-order orientation tensors. In addition, various models have more recently been developed to account for rotary diffusion due to fiber-fiber and fiber-matrix interactions which can now more accurately simulate the experimentally observed slow fiber kinematics in polymer composite processing. It is common to use explicit numerical initial value problem-ordinary differential equation (IVP-ODE) solvers such as the 4th- and 5th-order Dormand Prince Runge–Kutta (RK45) method to predict the transient and steady-state fiber orientation response. Here, we propose a computationally efficient method based on the Newton-Raphson (NR) iterative technique for determining steady state orientation tensor values by evaluating exact derivatives of the moment-tensor evolution equation with respect to the independent components of the orientation tensor. We consider various existing macroscopic-fiber orientation models and several closure approximations to ensure the robustness and reliability of the method. The performance and stability of the approach for obtaining physical solutions in various homogeneous flow fields is demonstrated through several examples. Validation of our orientation tensor exact derivatives is performed by benchmarking with results of finite difference techniques. Overall, our results show that the proposed NR method accurately predicts the steady state orientation for all tensor models, closure approximations and flow types considered in this paper and was relatively faster compared to the RK45 method. The NR convergence and stability behavior was seen to be sensitive to the initial orientation tensor guess value, the fiber orientation tensor model type and complexity, the flow type and extension to shear rate ratio.

Full text

Turn on search term navigation

Translate

1. Introduction

Characterization and evaluation of fiber suspensions has received considerable attention over the past four decades, particularly in the area of short fiber polymer composites produced in flow processes such as extrusion and molding. Understanding the motion and orientation of suspended fibers is critical for predicting the thermal, mechanical, and electrical performance of products made of these materials. Fiber orientation within a fiber suspension is often simulated using fiber orientation tensors, which provide a metric for assessing the random nature of sets of fibers and the degree of alignment. The goal of this paper is to compute the steady state orientation tensor of a fiber suspension in a homogeneous flow field using a Newton-Raphson iteration method with exact Jacobians.

Common tensor-based models used to compute the fiber orientation of a fiber suspension in polymer composite flows originated in the pioneering work by Jeffery [1], which described the evolution of a single rigid ellipsoid in a purely viscous flow of a Newtonian fluid. Jeffery’s hydrodynamic (HD) single ellipsoid model is limited to dilute suspension and ignores the effect of momentum diffusion due to fiber-fiber interactions. Moreover, Jeffery’s assumption ignores the flexural and fracture behavior of the fiber and assumes no slip contact between the fiber and surrounding fluid. Studies have shown that fibers suspended in industrial polymer melt flows tend to orient relative to the flow field [2], which is not captured in Jeffery’s work. As a result, enhancements to Jeffery’s single-fiber model have been made to better capture the bulk behavior of fibers in semi-dilute and concentrated suspension which include fiber-fiber interaction.

Although theoretically feasible, it is computationally expensive and nearly impractical to simulate the behavior of every individual particle in an industrially relevant short fiber polymer suspension flow. Alternatively, Folgar and Tucker [3,4] introduced isotropic rotary diffusion (IRD) and an orientation probability distribution function (ODF) to model interacting fiber suspensions based on Jeffery’s single-fiber model. The time evolution of the ODF is defined by the Fokker-Plank equation for probability distribution function (PDF) of fiber orientation. Conventionally, a numerical method such as finite volume (cf. Bay [5]) and, more recently, a more computationally efficient exact spherical harmonics method (cf. Montgomery-Smith et al. [6]) has been used to solve the Folgar-Tuckers (FT) equation of change for fiber orientation, but these methods have yet to see significant application in molding or extrusion processes. Advani and Tucker [7] proposed a tensorial representation of fiber orientation evolution based on the series expansion of even order moments of the ODF. Their approach simplifies computations, which has led to the widespread use of orientation tensors as the preferred method of evaluating fiber orientation for short fiber polymer composites. The Advani-Tucker (AT) orientation evolution model is commonly used to compute the 2nd order orientation tensor $a_{i j}$ , which requires a closure approximation of the 4th order orientation tensor $a_{i j k l}$ . Due to experimentally observed disparity in the fiber orientation kinematics compared to those computed from the Advani and Tucker orientation model [7], more recent modifications have been proposed in an effort to reduce the rate of alignment in polymer melt flow. Huynh [8] applied a strain reduction factor (SRF) to slow the transient response of the orientation tensor, which unfortunately lacked material objectivity. Alternatively, Wang et al. [9] introduced a phenomenological reduced strain closure (RSC) model that applies the reduction factor solely to the spectral decomposed principal rates of the orientation tensor without modifying effecting the rotation tensor, thus maintaining material objectivity. Similarly, Tseng et al. [10,11] proposed a retarding principal rate (RPR) model that involved a coaxial correction to the FT model by assuming the intrinsic orientation kinetics (IOK) describing the behavior of the fiber suspension involved a nonlinear modification to the principal directions of the material derivative.

Prediction of fiber orientation with IRD-based models has been validated experimentally for short-fiber/thermoplastic composites (SFT) with a fiber length typically in the range of 0.2 mm to 0.4 mm [12]. However, for long-fiber/thermoplastic composites (LFT) (e.g., with fiber lengths above 10 mm), IRD models tend to be less accurate with decreased alignment. Modifications to the rotary diffusion term have been proposed for improving accuracy when predicting the components of the orientation tensor. Ranganathan et al. [13] assumed an isotropic rotary diffusivity that inversely varies with the degree of alignment of the orientation tensor which was implemented using a phenomenological interaction parameter that depends on the reciprocal of the inter fiber spacing. The applicability of their model is limited to the transient orientation state while being well suited for long-range fiber-fiber interaction. Unfortunately, their model is known to be unsuitable for LFT steady state orientation prediction as with other IRD models, likely due to its isotropic diffusivity.

Fan et al. [14] and Phan-Thien [15] were the first to propose an anisotropic rotary diffusion (ARD) moment-tensor model by replacing the FT scalar phenomenological interaction parameter with a second order anisotropic rotary diffusion tensor. In a similar manner, Koch [16] developed an ARD model suited for semi-dilute suspensions with an anisotropic spatial diffusion tensor that depends on the orientation state and the rate of deformation tensor. However, their model was based on the more complicated PDF form for the orientation tensor representation rather than the moment-tensor form and proved ineffective in LFT modelling. Phelps and Tucker [12] extended the work of Fan [14] and Phan-Thien et al. [15] by developing a more general moment-tensor anisotropic diffusion model that depends on the spatial diffusion and orientation tensors called the Phelps-Tucker (PT) model. The derivation of the spatial diffusion tensor was written as a function of the orientation state and rate of deformation tensor in a similar manner to that proposed by Hand [17]. Phelps’s model showed remarkable improvements in predicting orientation states of LFTs. Tseng et al. [18] proposed the improved anisotropic rotary diffusion model (iARD) which defines a two-parameter spatial diffusion tensor that couples the effect of fiber-matrix interaction and fiber-fiber interaction. Unfortunately, the iARD model lacked material objectivity, limiting its applicability. More recently, Tseng et al. [19] proposed the principal anisotropic rotary diffusion (pARD) model assuming a principal spatial diffusion tensor that corotates with the orientation tensor. Bakharev [20] proposed the moldflow anisotropic rotational diffusion (mARD) model based on a reduction of the terms from the generic ARD model by Phelps to linear terms only with a spatial diffusion tensor like Tseng’s model. Latz et al. [21] developed a fully coupled flow-orientation model for concentrated suspensions by replacing the diffusion term in the FT model with an effective collision tensor that incorporates both and isotropic diffusion interaction term and a topological exclusion interaction term based on a nematic ‘Onsager’ potential of non-Brownian Maier-Saupe form (NEM-MS). They found the influence of the topological interaction on the fiber orientation state to be flow dependent, having a significant effect on channel and contraction flows, and found a relatively lesser influence for flow around cylinders. Kugler et al. [22], Favaloro et al. [23], Agboola et al. [24] and Park and Park [25] present detailed review and comparison of existing fiber orientation models. The foregoing ARD models are applied in polymer composite industry and have been incorporated in mold-filling flow computations in injection molding process simulations [26,27,28,29,30,31,32].

All fiber orientation tensor evolution equations include a tensor of higher order that must be evaluated or approximated. The common approach is to introduce a closure approximation that defines a higher-order fiber orientation tensor as a function of lower order information. Closure approximations are either based on an analytical expression or formed through a curve fitting process. The non-fitted analytical closure approximations include the those constructed from combinations of lower-order orientation and tensors and include the linear (LIN), quadratic (QDR), and the hybrid (HYB) closures. The general class of Hinch and Leal’s [33] composite closure approximations precontracted with the deformation rate tensor also falls under this category. Fitted closures are derived in terms of independent constants which are obtained using test data from solutions to the PDF for multiple standard flow conditions. The fitted closures include Eigenvalue-Based Fitted (EBF) and Invariant-Based Fitted (IBF) closures [34,35,36,37]. Higher-order optimal fitted closures such as the Eigenvalue-Based Optimal Fitted (EBOF) and Invariant-Based Optimal Fitted (IBOF) closures [34,35,36,37], which are extensions to the EBF and IBF closures, respectively, have also been developed for improved model accuracy. More details on different closure approximations and their associated advantages and disadvantages is discussed in later sections and also in [22,38,39,40]. Other fitted closure approximations include the neural network (NN)-based closures by Jack et al. [41,42] and the 6th-order Invariant-based orthotropic fitted closure by Jack and Smith [37,43].

The steady state values of orientation tensors in homogeneous flows have traditionally been computed with time evolution numerical IVP-ODE techniques such as the 4th-order Runge–Kutta method or predictor-corrector methods (see, e.g., [7,40]). However, when the objective is to obtain a steady-state orientation state, the transient solution process can require excessive computations for orientation tensor solution that are not needed. The main objective of this article is to present a computationally efficient approach for determining the steady state orientation tensors of fiber suspension in Newtonian homogenous flow field based on Newton-Raphson iteration with exact Jacobians. This approach zeros the time derivative of the 2nd order orientation tensor using explicit derivatives of the 2nd order orientation tensor rate of change equation with respect to its 2nd order orientation tensor components. Here, we consider various fiber orientation models and closure approximations and compare their performance in homogeneous flow fields. We benchmark the results of the explicit derivatives with those obtained using finite differences to ensure accuracy, and our Newton-Raphson iteration results are compared with time evolution numerical solutions in terms of both accuracy and efficiency. With the growing application of short carbon fiber-reinforced polymer composites in the additive manufacturing industry, the development of reliable and computationally efficient modelling tools for modelling fiber suspension in order to better understand the development of the inherent bead microstructure is crucial. Our NR technique can be incorporated into existing commercially available short fiber orientation modelling software as a viable alternative to conventional numerical IVP-ODE techniques for determining the fiber orientation steady state which can improve the overall computational performance.

2. Materials and Methods

2.1. Representation of Fiber Orientation State

Figure 1 shows an orientation state of a prolate ellipsoid including the local coordinate systems $x_{1}^{'}$ , $x_{2}^{'}$ , $x_{3}^{'}$ and Euler angles definition where the orientation vector $\underline{p}$ lies along the major axis of the ellipsoid.

The orientation state of a fiber, as shown in Figure 1, can be described by the unit vector $\underline{p}$ as follows [7]:

(1) $\underline{p} = {[\begin{matrix} \sin θ \cos ϕ & \sin θ \sin ϕ & \cos θ \end{matrix}]}^{T}$

The orientation of the unit vector

\underline{p}

is commonly described in terms of the probability distribution function (PDF)

ψ (\underline{p})

over all possible directions of

\underline{p}

. It can be shown that

ψ (\underline{p}) = ψ (- \underline{p})

satisfies the normalization condition.

(2) $\oint ψ (\underline{p}) d \underline{p} = \int_{θ = 0}^{π} \int_{ϕ = 0}^{2 π} ψ (θ, ϕ) \sin θ d θ d ϕ = 1$

The PDF

ψ (\underline{p})

also satisfies the continuity condition [7].

(3) $\frac{D ψ}{D t} = - \frac{\partial}{\partial \underline{p}} (ψ \dot{\underline{p}})$

A set of higher-order, even order tensors, can also be used to describe the orientation state of a single fiber or fiber suspension such as the 2nd and 4th-order orientation tensors, given, respectively, as

(4) $a_{i j} = \oint p_{i} p_{j} ψ (\underline{p}) d \underline{p} a n d a_{i j k l} = \oint p_{i} p_{j} p_{k} p_{l} ψ (\underline{p}) d \underline{p}$

Orientation tensors in Equation (4) are completely symmetric, i.e.,

(5) $a_{i j} = a_{j i} a_{i j k l} = a_{j i k l} = a_{k i j l} = a_{l i j k} = a_{i k j l} = a_{i l j k} = \dots, 24 p e r m u t a t i o n s$

and a normalization condition requires that

(6) $a_{i i} = 1, a_{i j k k} = a_{i j}$

Consequently, there are only 5 independent components in

a_{i j}

and 9 independent components in

a_{i j k l}

2.2. Determination of Steady State Orientations

The numerical approach developed here for determining steady state values of the orientation vector $p_{i}$ and orientation tensor components $a_{i j}$ computes a zero rate of change of each corresponding evolution equation using the Newton-Raphson method. The rate of change equation (often referred to as the residual, $R_{i}$ ) for both $p_{i}$ and $a_{i j}$ are set to zero, respectively:

(7) $R_{i} = \frac{D p_{i}}{D t} = 0 a n d R_{i j} = \frac{D a_{i j}}{D t} = 0$

where

i, j \in 1,2, 3

, and summation convention is implied unless otherwise stated. Newton-Raphson iteration computes the root (in our case the orientation state) as [44]

(8) $p_{i}^{+} = p_{i}^{-} - [J_{i j} (p_{i}^{-})]^{- 1} R_{j} (p_{i}^{-})$

for the orientation vector and

(9) ${a_{i j}}^{+} = {a_{i j}}^{-} - [J_{m n i j} ({a_{i j}}^{-})]^{- 1} R_{m n} ({a_{i j}}^{-})$

for the 2nd order orientation tensor, where the

(+)

and (

-)

superscripts denote the current and previous iteration, respectively. The implication of Equations (8) and (9) is the need to compute the Jacobians

J_{i j}

and

J_{m n i j}

, which are the derivatives of the orientation vector and tensor, respectively, with respect to its components that define them, i.e.,

(10) $J_{i j} = \frac{\partial R_{i}}{\partial p_{j}} a n d J_{m n i j} = \frac{\partial R_{m n}}{\partial a_{i j}}$

The above computations, including the 2nd and 4th order tensor inverse calculations, are conveniently performed using contracted tensor notation. Since there are only 5 independent components of

a_{i j}

(see discussion below), the residual

R_{m n}

may be represented in contracted form as the vector

R_{r}

(and similarly, the tensor

a_{i j}

may be represented as the vector

a_{r}

), and the Jacobian

J_{m n i j}

as a matrix

J_{r s}

. The contracted notation is related to the tensor components through index mapping.

(11) $r (m, n) = n - \frac{1}{2} (m - 1) (m - 6), f o r n = m \dots 3, f o r m = 1, 2$

It follows that Equation (9) in contracted notation becomes

(12) ${a_{r}}^{+} = {a_{r}}^{-} - {[J_{r s} ({a_{r}}^{-})]}^{- 1} R_{r} ({a_{r}}^{-})$

In the following sections, we present existing models for rate of change equations of the orientation vector

p_{i}

and tensor

a_{i j}

based on those considered in the review by Kugler et al. [22]. Equating each rate of change equation to zero yields the Newton-Raphson residual in Equation (1), which is then differentiated to obtain expressions for the associated Jacobian.

2.3. Fiber Orientation Modelling in the Dilute Regime

Jeffery’s hydrodynamic model for the motion of a single rigid ellipsoidal particle in an incompressible Newtonian viscous fluid flow field forms the basis for evaluating fiber orientation for dilute suspensions. (Dilute suspension is a heterogenous suspension where the average interparticle spacing is relatively large such that there is no restriction on the fibers motion due to hydrodynamic forces or mechanical contact. More details on fiber suspension concentration regimes can be found in [45,46].) Jeffery assumed that an ellipsoidal particle is convected with the bulk motion of the undisturbed surrounding fluid where the rate of change of $\underline{p}$ written in component form is given as [1,22,23]

(13) ${\dot{p}}_{i} = ω_{i j} p_{j} + ξ ({\dot{γ}}_{i j} p_{j} - {\dot{γ}}_{k l} p_{k} p_{l} p_{i})$

where

ω_{i j}

and

{\dot{γ}}_{i j}

are the anti-symmetric and symmetric decomposition of the rate of deformation tensor

L_{i j}

, given as

(14) $ω_{i j} = \frac{1}{2} (\frac{\partial v_{i}}{\partial x_{j}} - \frac{\partial v_{j}}{\partial x_{i}}) a n d {\dot{γ}}_{i j} = \frac{1}{2} (\frac{\partial v_{i}}{\partial x_{j}} + \frac{\partial v_{j}}{\partial x_{i}})$

such that

L_{i j} = \partial v_{i} / \partial x_{j} = {\dot{γ}}_{i j} + ω_{i j}

, and

ξ

is a particle shape parameter given as

ξ = (r_{e}^{2} - 1) / (r_{e}^{2} + 1)

for a particle with equivalent ellipsoidal aspect ratio

r_{e}

. It is worth noting that

r_{e}

can be related to the geometric aspect ratio of various axisymmetric particle shapes such as cylindrical particles [47]. We define components

R_{i}

of the Newton-Raphson residual in Equation (8) for computing the steady-state orientation of the vector

\underline{p}

in Equation (9) as

(15) $R_{i}^{J F} = {\dot{p}}_{i}$

where the superscript

J F

denotes the ‘Jeffery’ type residual. The Jacobian

J_{m n}^{J F} = \frac{\partial R_{m}^{J F}}{\partial p_{n}}

is obtained by taking derivatives of the components

R_{m}^{J F}

with respect to

p_{n}

, i.e.,

(16) $J_{m n}^{J F} = \frac{\partial {\dot{p}}_{m}}{\partial p_{n}} = ω_{m j} δ_{j n} + ξ [{\dot{γ}}_{m j} δ_{j n} - {\dot{γ}}_{k l} (δ_{k n} p_{l} p_{m} + p_{k} δ_{l n} p_{m} + p_{k} p_{l} δ_{m n})]$

where we note that derivatives of

\underline{p}

, with respect to itself, form the identity tensor., i.e.,

\partial p_{i} / \partial p_{j} = δ_{i j}

. Equation (16) contains only 2 independent components of the orientation vector

R_{i}^{J F}

, making it possible to represent the components of

J_{m n}^{J F}

as a

2 \times 2

matrix. Although Jeffery’s model for dilute suspensions has limited direct application since fiber-fiber interaction is ignored, the NR method can be used to find the steady state solution for flows with a high extension to shear rate ratio, such as those considered in [48].

2.4. Fiber Orientation Modelling in Semi-Dilute and Concentrated Regime

In the concentrated regime, the average interparticle spacing is very small (orders of magnitude less than the smallest particle dimension) such that the fiber motion is affected by hydrodynamic forces and possible direct mechanical contact with other particles. Definitions for the semi-dilute and concentrated fiber suspension regimes can be found in [45,46]. The effect of the interparticle interaction on the single particle’s motion is typically modeled by including a momentum diffusion term in Equation (13). Here, we consider various diffusion models which are commonly used in commercial SFRP simulations having fiber suspensions in the concentrated regime. In the following sections, the NR method for resolving the steady state fiber orientation is developed for multiple 2nd order fiber orientation tensor and rotary diffusion models and include several 4th order closure approximations.

2.4.1. The Folgar-Tucker Diffusion Model

A model that accounts for the effect of momentum diffusion due to short- and long-range fiber-fiber interaction in non-dilute fiber suspensions was first proposed by Folgar-Tucker model [3,4] by incorporating a rotary diffusion term into Jeffery’s single fiber model as

(17) ${\dot{p}}_{i}^{F T} = {\dot{p}}_{i}^{J F} - D_{r} \frac{1}{ψ} \frac{\partial ψ}{\partial p_{i}}$

where

D_{r}

is the rotary diffusivity constant that introduces a Brownian diffusion effect among contacting particles. For long slender particles (i.e.,

ξ \approx 1

), Folgar and Tucker set

D_{r} = C_{I} \dot{γ}

, where

C_{I}

is a phenomenological interaction coefficient and

\dot{γ} = \sqrt{2 {\dot{γ}}_{i j} {\dot{γ}}_{j i}}

is the scalar magnitude of the strain rate tensor

{\dot{γ}}_{i j}

. The orientation PDF

ψ

in Equation (17) defines the probability that a given fiber has a particular orientation (see [3,4,7]). Advani and Tucker [7] derived the orientation tensor form of the Folger-Tucker model defined in terms of the 2nd and 4th order orientation tensors. Advani and Tucker developed an equation for the rate of change of

a_{i j}

, identified as the orientation material derivative tensor

{\dot{a}}_{i j}

in terms of the 2nd and 4th order tensor:

(18) $\frac{D a_{i j}^{A T}}{D t} = {\dot{a}}_{i j}^{A T} = {\dot{a}}_{i j}^{H D} + {\dot{a}}_{i j}^{I R D}$

where

{\dot{a}}_{i j}^{H D}

is the hydrodynamic component derived from Jeffery’s equation (cf. Equation (11)) given as

(19) ${\dot{a}}_{i j}^{H D} = - (ω_{i k} a_{k j} - a_{i k} ω_{k j}) + ξ ({\dot{γ}}_{i k} a_{k j} - a_{i k} {\dot{γ}}_{k j} - 2 {\dot{γ}}_{k l} a_{i j k l})$

and

{\dot{a}}_{m n}^{I R D}

is the added isotropic rotary diffusion:

(20) ${\dot{a}}_{i j}^{I R D} = 2 D_{r} (δ_{i j} - α a_{i j})$

In the above,

α

is a dimension factor (i.e.,

α = 3

for 3D orientation and

α = 2

for 2D planar orientation). It follows that the residual in Equation (8) for the Advani-Tucker model is

(21) $R_{m n}^{A T} = {\dot{a}}_{i j}^{A T}$

The associated Jacobian

J_{m n i j}^{F T}

for the Advani-Tucker model is obtained by differentiating

R_{m n}^{A T}

with respect to the components of

a_{i j}

, yielding

(22) $J_{m n i j}^{A T} = \frac{\partial R_{m n}^{A T}}{\partial a_{i j}} = \frac{\partial {\dot{a}}_{m n}^{H D}}{\partial a_{i j}} + \frac{\partial {\dot{a}}_{m n}^{I R D}}{\partial a_{i j}}$

Given that the derivative of

a_{i j}

with respect to itself is

(23) $\frac{\partial a_{m n}}{\partial a_{i j}} = δ_{m i} δ_{n j}$

the derivative terms on the right-hand side of Equation (22) may be written as

(24) $\frac{\partial {\dot{a}}_{m n}^{H D}}{\partial a_{i j}} = (- ω_{m k} + ξ {\dot{γ}}_{m k}) δ_{k i} δ_{n j} + (ω_{k n} + ξ {\dot{γ}}_{k n}) δ_{m i} δ_{k j} - 2 ξ {\dot{γ}}_{k l} \frac{\partial a_{m n k l}}{\partial a_{i j}}$

(25) $\frac{\partial {\dot{a}}_{m n}^{I R D}}{\partial a_{i j}} = - 2 D_{r} α δ_{m i} δ_{n j}$

Expressions for derivatives of

a_{m n k l}

(i.e.,

\frac{\partial a_{m n k l}}{\partial a_{i j}}

) in Equation (24) and others using well known closure approximations are presented in the following sections.

2.4.2. Strain Reduction Factor (SRF) Model

The SRF model was developed by Huynh [8] to reduce the rate of alignment as compared to the AT model by introducing a strain reduction factor $κ$ ( $0 < κ < 1$ ) into Equation (16) (i.e., multiplying the right-hand-side of Equation (16) by $κ$ ) to slow down the orientation kinetics as observed experimentally. He based his premise on a reduced bulk strain of fiber clusters in a concentrated suspension flow. The residual and Jacobian in this case is just a multiplication of $κ$ with that previously obtained for the FT model.

(26) $R_{m n}^{S R F} = κ R_{m n}^{A T}, J_{m n i j}^{F T} = κ J_{m n i j}^{A T}, 0 < κ < 1$

Steady state orientations predicted by the SRF model have been experimentally validated [9]; however, it tends to overshoot measured alignments in the initial transient. The SRF model also does not satisfy the rheological test of material objectivity and, therefore, cannot be applied to complex flows.

2.4.3. Reduced Strain Closure (RSC) Model

Wang et al. [9] addressed material objectivity with their reduced strain closure (RSC) model which applies a reduction factor only to the evolution rate of the spectral decomposed principal directions $\underline{\underline{Φ}}$ of the orientation tensor $\underline{\dot{λ}}$ , without modifying the rate of rotation $\dot{\underline{\underline{Φ}}}$ , given as

(27) ${\dot{λ}}_{i}^{R S C} = κ {\dot{λ}}_{i}^{A T}, {\dot{Φ}}_{i j}^{R S C} = {\dot{Φ}}_{i j}^{A T}, a_{m n} | a_{m n} = λ_{i} Φ_{m i} Φ_{n i}$

where

|

denotes the mathematical abbreviation for ‘such that’, and repeated indices imply summation. It follows that the material derivative for the RSC model may be written as [9]

(28) ${\dot{a}}_{m n}^{R S C} = {\dot{a}}_{m n}^{A T} - (1 - κ) {\dot{a}}_{m n}^{Δ A T}$

where the kernel of

{\dot{a}}_{m n}^{Δ A T}

defines the part of the standard AT model affected by the modification factor

κ

to yield the RSC model, given as

(29) ${\dot{a}}_{m n}^{Δ A T} = {\dot{λ}}_{i}^{A T} Φ_{m i} Φ_{n i} = 2 ξ {\dot{γ}}_{k l} (L_{m n k l} - M_{m n r s} a_{r s k l}) + {\dot{a}}_{m n}^{I R D}$

where

(30) $L_{m n k l} = {\dot{λ}}_{i} Φ_{m i} Φ_{n i} Φ_{k i} Φ_{l i}, M_{m n k l} = Φ_{m i} Φ_{n i} Φ_{k i} Φ_{l i}$

The Newton-Raphson residual for the RSC model is $R_{m n}^{R S C} = {\dot{a}}_{m n}^{R S C}$ , and the Jacobian is obtained by taking partial derivatives of Equation (31) as follows:

(31) $J_{m n i j}^{R S C} = J_{m n i j}^{F T} - (1 - κ) \frac{\partial {\dot{a}}_{m n}^{Δ A T}}{\partial a_{i j}}$

where

(32) $\frac{\partial {\dot{a}}_{m n}^{Δ A T}}{\partial a_{i j}} = 2 ξ {\dot{γ}}_{k l} \frac{\partial}{\partial a_{i j}} \{L_{m n k l} - M_{m n r s} a_{r s k l}\} + \frac{\partial {\dot{a}}_{m n}^{I R D}}{\partial a_{i j}}$

which can be re-written as

(33) $\frac{\partial {\dot{a}}_{m n}^{Δ A T}}{\partial a_{i j}} = 2 ξ {\dot{γ}}_{k l} [\frac{\partial L_{m n k l}}{\partial a_{i j}} - a_{r s k l} \frac{\partial M_{m n r s}}{\partial a_{i j}} - M_{m n r s} \frac{\partial a_{r s k l}}{\partial a_{i j}}] + \frac{\partial {\dot{a}}_{m n}^{I R D}}{\partial a_{i j}}$

In the above, the derivatives of 4th order tensors

M_{m n k l}

and

L_{m n k l}

are obtained by differentiating Equation (33), respectively, as

(34) $\begin{matrix} \frac{\partial M_{m n k l}}{\partial a_{r s}} = \frac{\partial}{\partial a_{r s}} & \{Φ_{m i} Φ_{n i} Φ_{k i} Φ_{l i}\} \\ = Φ_{n i} Φ_{k i} Φ_{l i} \frac{\partial Φ_{m i}}{\partial a_{r s}} + Φ_{m i} Φ_{k i} Φ_{l i} \frac{\partial Φ_{n i}}{\partial a_{r s}} + Φ_{m i} Φ_{n i} Φ_{l i} \frac{\partial Φ_{k i}}{\partial a_{r s}} \\ + Φ_{m i} Φ_{n i} Φ_{k i} \frac{\partial Φ_{l i}}{\partial a_{r s}} \end{matrix}$

and

(35) $\frac{\partial L_{m n k l}}{\partial a_{r s}} = Φ_{m i} Φ_{n i} Φ_{k i} Φ_{l i} \frac{\partial {\dot{λ}}_{i}}{\partial a_{r s}} + {\dot{λ}}_{i} \frac{\partial}{\partial a_{r s}} \{Φ_{m i} Φ_{n i} Φ_{k i} Φ_{l i}\}$

where the derivatives of the eigenvalues

λ_{i}

and eigenvectors

Φ_{j k}

appearing in Appendix A have been used to simplify the final expression.

2.4.4. Retarding Principal Rate (RPR) Model

Tseng et al. [10,11] developed a Retarding Principal Rate (RPR) model, which, like the RSC model, reduces fiber orientation kinetics. Their approach is based on a coaxial modification to the principal directions of the orientation tensor evolution rate via a nonlinear correlation. The material derivative tensor ${\dot{a}}_{m n}^{X}$ of any standard model X (where X can define the AT, RSC, or other orientation tensor evolution model) is linearly combined with the RPR correction to slow down the orientation response rate, i.e.,

(36) ${\dot{a}}_{m n}^{X - R P R} = {\dot{a}}_{m n}^{X} + {\dot{a}}_{m n}^{R P R}$

where the RPR correction

{\dot{a}}_{m n}^{R P R}

is written in terms of its eigenvalue derivative matrix

{\dot{Λ}}_{k l}^{I O K}

(37) ${\dot{a}}_{m n}^{R P R} = - Φ_{m k} {\dot{Λ}}_{k l}^{I O K} Φ_{n l}, {\dot{Λ}}_{k l}^{I O K} = {\dot{Λ}}_{k l}^{I O K} ({\underline{\dot{λ}}}^{X})$

where

{\dot{Λ}}_{k l}^{I O K}

is the correction to the principal values of the orientation tensor material derivative based on the IOK assumption (see, e.g., Equation (41) below). The superscript on

{\dot{Λ}}_{k l}^{I O K}

indicates the intrinsic orientation kinetics (IOK) assumption. Given that the correction is coaxial, the rotation tensor growth rate is unaffected and is obtained from the spectral decomposition of

a_{m n}

, where

\underline{\underline{Φ}}

satisfies

(38) $Λ_{m n} = Φ_{k m} a_{k l} Φ_{l n}$

The columns of the eigenmatrix obtained in this way are reordered in descending order with respect to the magnitude of the eigenvalues, i.e.,

(39) $Φ_{i j} = \{Φ_{i j} | λ_{j} = Λ_{j j}, λ_{1} \geq λ_{2} \geq λ_{3}\}$

where

|

denotes the mathematical abbreviation for ‘such that’ and the repeated indices in

Λ_{j j}

do not imply summation but rather represent the

j t h

diagonal component of

Λ_{i j}

. The growth rate of the principal eigenvalues of the standard model X, i.e.,

{\dot{Λ}}_{k l}^{X}

, is obtained from

(40) ${\dot{Λ}}_{k l}^{X} = Φ_{k m} {\dot{a}}_{k l}^{X} Φ_{l n}, {\dot{λ}}_{k}^{X} = {\dot{Λ}}_{k k}^{X}$

The diagonal tensor

{\dot{Λ}}_{k k}^{I O K}

is defined by a 2-parameter non-linear correlation to the principal values of the standard model

{\dot{Λ}}_{k l}^{X}

such that

(41) ${\dot{Λ}}_{k k}^{I O K} = {\dot{λ}}_{k}^{I O K} = α [{\dot{λ}}_{k}^{X} - β ({\{{\dot{λ}}_{k}^{X}\}}^{2} + 2 {\dot{λ}}_{l}^{X} {\dot{λ}}_{m}^{X})], {{\dot{Λ}}_{k l}^{I O K}|}_{k \neq l} = 0$

For an RPR corrected model, the NR residual

R_{m n}^{X - R P R}

is simply the material derivative

(42) $R_{m n}^{X - R P R} = {\dot{a}}_{m n}^{X - R P R}$

and the Jacobian

J_{m n i j}^{X - R P R}

is given as

(43) $J_{m n i j}^{X - R P R} = \frac{\partial {\dot{a}}_{m n}^{X}}{\partial a_{i j}} + \frac{\partial {\dot{a}}_{m n}^{R P R}}{\partial a_{i j}}$

The partial derivative of the RPR correction term

{\dot{a}}_{m n}^{R P R}

(44) $\frac{\partial {\dot{a}}_{m n}^{R P R}}{\partial a_{i j}} = - \{\frac{\partial Φ_{m k}}{\partial a_{i j}} {\dot{Λ}}_{k l}^{I O K} Φ_{n l} + Φ_{m k} \frac{\partial {\dot{Λ}}_{k l}^{I O K}}{\partial a_{i j}} Φ_{n l} + Φ_{m k} {\dot{Λ}}_{k l}^{I O K} \frac{\partial Φ_{n l}}{\partial a_{i j}}\}$

and the partial derivative of the modified growth rate of eigenvalues tensor

{\dot{Λ}}_{k l}^{I O K}

is obtained by taking the derivative of Equation (41) as

(45) $\frac{\partial {\dot{Λ}}_{k k}^{I O K}}{\partial a_{i j}} = \frac{\partial {\dot{λ}}_{k}^{I O K}}{\partial a_{i j}} = α [\frac{\partial {\dot{λ}}_{k}^{X}}{\partial a_{i j}} - 2 β ({\dot{λ}}_{k}^{X} \frac{\partial {\dot{λ}}_{k}^{X}}{\partial a_{i j}} + \frac{\partial {\dot{λ}}_{l}^{X}}{\partial a_{i j}} {\dot{λ}}_{m}^{X} + {\dot{λ}}_{l}^{X} \frac{\partial {\dot{λ}}_{m}^{X}}{\partial a_{i j}})], {\frac{\partial {\dot{Λ}}_{k l}^{I O K}}{\partial a_{i j}}|}_{k \neq l} = 0$

where

(46) $\frac{\partial {\dot{Λ}}_{k l}^{X}}{\partial a_{i j}} = \{\frac{\partial Φ_{k m}}{\partial a_{i j}} {\dot{a}}_{k l}^{X} Φ_{l n} + Φ_{k m} \frac{\partial {\dot{a}}_{k l}^{X}}{\partial a_{i j}} Φ_{l n} + Φ_{k m} {\dot{a}}_{k l}^{X} \frac{\partial Φ_{l n}}{\partial a_{i j}}\}, \frac{\partial {\dot{λ}}_{k}^{X}}{\partial a_{i j}} = \frac{\partial {\dot{Λ}}_{k k}^{X}}{\partial a_{i j}}$

In the above,

\frac{\partial {\dot{a}}_{m n}^{X}}{\partial a_{i j}}

is the partial derivative of material derivative

{\dot{a}}_{m n}^{X}

with respect to the components of

a_{i j}

obtained a priori, and the partial derivatives of the eigenvector matrix with respect to the same

(i . e ., \frac{\partial Φ_{m n}}{\partial a_{i j}})

can be obtained as shown in Appendix A (cf. [49]).

2.4.5. Anisotropic Rotary Diffusion (ARD) Models

While the IRD models were experimentally observed to be accurate in predicting the orientation state of short-fiber/thermoplastic composites (SFT), they were ineffective in accurately predicting the complete set of orientation tensor components for the long-fiber/thermoplastic composites (LFT), which motivated ARD model development. Various ARD models and modified forms have been developed based on the definition of the spatial diffusion tensor. Most models utilize the moment-tensor form for the ARD representation developed by Phelps and Tucker [12]. The general expression for the 2nd order orientation tensor evolution rate is a linear combination of the Jeffery’s-based model and the rotary diffusion term, given as

(47) ${\dot{a}}_{m n}^{P T} = {\dot{a}}_{m n}^{H D} + {\dot{a}}_{m n}^{A R D}$

where the rotary diffusion term

{\dot{a}}_{m n}^{A R D}

is defined in terms of the spatial diffusion coefficient and the orientation state, given as

(48) ${\dot{a}}_{m n}^{A R D} = \dot{γ} [2 C_{m n} - 2 C_{r s} δ_{r s} a_{m n} - 5 (C_{m k} a_{k n} + a_{m k} C_{k n}) + 10 a_{m n k l} C_{k l}]$

and

\underline{\underline{C}}

is the spatial diffusion tensor. Based on this model, the NR residual and Jacobian are, respectively, given as

(49) $R_{m n}^{P T} = {\dot{a}}_{m n}^{P T}$

(50) $J_{m n i j}^{P T} = \frac{\partial {\dot{a}}_{m n}^{P T}}{\partial a_{i j}} = \frac{\partial {\dot{a}}_{m n}^{H D}}{\partial a_{i j}} + \frac{\partial {\dot{a}}_{m n}^{A R D}}{\partial a_{i j}}$

where the derivative of the rotary diffusion (

A R D)

term is obtained using product rule as follows:

(51) $\begin{matrix} \frac{\partial {\dot{a}}_{m n}^{A R D}}{\partial a_{i j}} = \dot{γ} [2 \frac{\partial C_{m n}}{\partial a_{i j}} & - 2 (\frac{\partial C_{r s}}{\partial a_{i j}} δ_{r s} a_{m n} + C_{r s} δ_{r s} δ_{m i} δ_{n j}) + \dots \\ - 5 (\frac{\partial C_{m k}}{\partial a_{i j}} a_{k n} + C_{m k} δ_{k i} δ_{n j} + δ_{m i} δ_{k j} C_{k n} + a_{m k} \frac{\partial C_{k n}}{\partial a_{i j}}) + \dots \\ + 10 (\frac{\partial a_{m n k l}}{\partial a_{i j}} C_{k l} + a_{m n k l} \frac{\partial C_{k l}}{\partial a_{i j}})] \end{matrix}$

In the Mold-Flow ARD (mARD) model developed by Bakharev [20], the Phelps and Tucker’s rotary diffusion ARD expression is truncated to include just the linear terms, i.e.,

(52) ${\dot{a}}_{m n}^{m A R D} = \dot{γ} [2 C_{m n} - 2 C_{k l} δ_{k l} a_{m n}]$

(53) $\frac{\partial {\dot{a}}_{m n}^{m A R D}}{\partial a_{i j}} = \dot{γ} [2 \frac{\partial C_{m n}}{\partial a_{i j}} - 2 (\frac{\partial C_{k l}}{\partial a_{i j}} δ_{k l} a_{m n} + C_{k l} δ_{k l} δ_{m i} δ_{n j})]$

The corresponding evolution rate equation for the 2nd order orientation tensor based on mARD model is given as

(54) ${\dot{a}}_{m n}^{m P T} = {\dot{a}}_{m n}^{H D} + {\dot{a}}_{m n}^{m A R D}$

Various forms of the spatial diffusion coefficient

C_{m n}

used in the ARD model (indicated by a superscript in the following) have been developed. The Phelps and Tucker [12] representation

C_{m n}^{P T}

based on a modification to Hand’s anisotropic tensor [17] is given as a function of the rate of deformation tensor

{\dot{γ}}_{i j}

and orientation state

a_{i j}

(55) $C_{m n}^{P T} = b_{1} δ_{m n} + b_{2} a_{m n} + b_{3} a_{m k} a_{n k} + \frac{b_{4}}{\dot{γ}} {\dot{γ}}_{m n} + \frac{b_{5}}{{\dot{γ}}^{2}} {\dot{γ}}_{m k} {\dot{γ}}_{n k}$

where

b_{i}

are dimensionless constants and are determined experimentally. The derivative of the

C_{m n}^{P T}

with respect to

a_{i j}

that is substituted into Equation (56) is

(56) $\frac{\partial C_{m n}^{P T}}{\partial a_{i j}} = b_{2} δ_{m i} δ_{n j} + b_{3} (δ_{m i} δ_{k j} a_{n k} + a_{m k} δ_{n i} δ_{k j})$

The sensitivity of the numerical stability of the model response on the PT model parameters

b_{i}

coupled with the complicated process involved in their determination are major limitations to the application of this model. Tseng et al. [18] developed an improved anisotropic rotary diffusion model (iARD) with a two-parameter spatial diffusion tensor representation written in terms of

{\dot{γ}}_{i j}

and its scalar magnitude

\dot{γ}

that couples the effect of fiber-matrix interaction and fiber-fiber interaction, given as

(57) $C_{m n}^{i A R D} = C_{I} (δ_{m n} - 4 C_{M} \frac{{\dot{γ}}_{m k} {\dot{γ}}_{n k}}{{\dot{γ}}^{2}})$

where

C_{I}

and

C_{M}

are the fiber-fiber and fiber-matrix interaction parameters, respectively. An alternate form of Equation (60) is given as

(58) $C_{m n}^{i A R D} = C_{I} (δ_{m n} - C_{M} {\tilde{L}}_{m n}), {\tilde{L}}_{m n} = (L_{m k} L_{n k}) / (L_{r s} L_{r s})$

where

L_{m n}

is the rate of deformation tensor. The derivative

C_{m n}^{i A R D}

with respect to

a_{i j}

is simply zero, i.e.,

(59) $\frac{\partial C_{m n}^{i A R D}}{\partial a_{i j}} = 0$

Due to the lack of material objectivity of the rate of deformation tensor

L_{m n}

, in Equation (61), Tseng et al. [19] developed an improved objective principal spatial tensor ARD model (pARD) that is coaxial with the orientation tensor, given as

(60) $C_{m n}^{p A R D} = \{C_{I} Φ_{m k} D_{k l} Φ_{n l}, \underline{\underline{Φ}} | a_{m n} = Φ_{m k} {\bar{a}}_{k l} Φ_{n l}\}$

where the tensor

D_{k l}

contains only diagonal terms, and its trace is unity, i.e.,

(61) $D_{k k} = 1, {D_{k l}|}_{k \neq l} = 0$

The derivative of

C_{m n}^{p A R D}

with respect to the 2nd order orientation tensor for substitution into Equation (56) is

(62) $\frac{\partial C_{m n}^{p A R D}}{\partial a_{i j}} = C_{I} \{\frac{\partial Φ_{m k}}{\partial a_{i j}} D_{k l} Φ_{n l} + Φ_{m k} D_{k l} \frac{\partial Φ_{n l}}{\partial a_{i j}}\}$

Another ARD model form suggested by Wang [50] called the Wang-Phelps-Tucker (WPT) model involved truncating the terms of the PT model to obtain

(63) $C_{m n}^{W P T} = b_{1} δ_{m n} + b_{3} a_{m k} a_{n k}$

Falvoro et al. [23] provided an alternative form of

C_{m n}^{W P T}

, where he replaced the coefficients in Equation (66) with a weighted superposition of the interaction coefficient

C_{I}

, written as

(64) $C_{m n}^{W P T} = C_{I} ((1 - w) δ_{m n} + w a_{m k} a_{n k})$

where

w

is the weighting factor. The derivative of

C_{m n}^{p A R D}

with respect to

a_{i j}

for substitution into Equation (56) is

(65) $\frac{\partial C_{m n}^{W P T}}{\partial a_{i j}} = w C_{I} (δ_{m k} a_{n k} + a_{m k} δ_{n k})$

Lastly, we consider the

D_{z}

ARD model (cf. Falvoro et al. [23]) developed by Moldflow for simulating

2.5

D flow processes. Their model is defined in terms of the interaction coefficient

C_{I}

, a moment of interaction thickness parameter

D_{z}

, and the unit normal to the mold surface

\hat{n}

. The expression for

C_{m n}

for the

D_{z}

ARD model is

(66) $C_{m n}^{D z} = C_{I}^{D z} (δ_{m n} - (1 - D_{z}) {\hat{n}}_{m} {\hat{n}}_{n})$

and

(67) $\frac{\partial C_{m n}^{D z}}{\partial a_{i j}} = 0$

2.4.6. Nematic Potential (NEM) Model

Latz et al. [21] proposed a 2-parameter nematic potential ARD (NEM) model that couples the phenomenological effect of the momentum diffusion due to fiber-fiber interaction and a topological interaction effect of diffusion due to an exclusion volume mechanism. The material derivative of the 2nd order orientation tensor based on the nematic diffusion model is

(68) ${\dot{a}}_{m n}^{n e m} = {\dot{a}}_{m n}^{H D} + {\dot{a}}_{m n}^{I R D - M S}$

where

(69) ${\dot{a}}_{m n}^{I R D - M S} = \dot{γ} [C_{I} (δ_{m n} - α a_{m n}) + U_{0} (a_{m k} a_{k n} - a_{k l} a_{m n k l})]$

In the above,

U_{0}

is the ‘Onsager’ nematic topological interaction coefficient of the Maier-Saupe (MS) potential. Typically, for stability,

U_{0} \leq 4 C_{I}

for 2D analysis and

U_{0} > 8 C_{I}

for 3D analysis. The NR residual and Jacobian are, respectively, given as

(70) $R_{m n}^{n e m} = {\dot{a}}_{m n}^{n e m}$

(71) $J_{m n i j}^{n e m} = \frac{\partial {\dot{a}}_{m n}^{H D}}{\partial a_{i j}} + \frac{\partial {\dot{a}}_{m n}^{I R D - M S}}{\partial a_{i j}}$

where the derivative of the nematic diffusion term is obtained from Equation (72) as

(72) $\frac{\partial {\dot{a}}_{m n}^{I R D - M S}}{\partial a_{i j}} = \dot{γ} [- C_{I} α δ_{m i} δ_{n j} + U_{0} (δ_{m i} δ_{k j} a_{k n} + a_{m k} δ_{k i} δ_{n j} - δ_{k i} δ_{l j} a_{m n k l} - a_{k l} \frac{\partial a_{m n k l}}{\partial a_{i j}})]$

It is common to combine multiple diffusion models when predicting the orientation state with the goal of improved accuracy. One such combination is the ARD-RSC models, whose material time derivative is expressed as [19]

(73) ${\dot{a}}_{m n}^{p A R D - R S C} = {\dot{a}}_{m n}^{R S C} - κ {\dot{a}}_{m n}^{I R D} + {\dot{a}}_{m n}^{A R D} + {\dot{a}}_{m n}^{Δ R S C}$

where

(74) ${\dot{a}}_{m n}^{Δ R S C} = - 2 \dot{γ} (1 - κ) [M_{m n k l} - δ_{k l} a_{m n} - 5 (L_{m n k l} - M_{m n r s} a_{r s k l})] C_{k l}$

and

{\dot{a}}_{m n}^{I R D}

{\dot{a}}_{m n}^{R S C}

, and

{\dot{a}}_{m n}^{A R D}

are given in Equations (20), (28) and (48), respectively. In this case the NR residual

R_{m n}^{p A R D - R S C}

is the material derivative

{\dot{a}}_{m n}^{p A R D - R S C}

, i.e.,

(75) $R_{m n}^{p A R D - R S C} = {\dot{a}}_{m n}^{p A R D - R S C}$

and the Jacobian is obtained by taking the partial derivative of Equation (76) with respect to

a_{i j}

(76) $J_{m n i j}^{p A R D - R S C} = \frac{\partial {\dot{a}}_{m n}^{R S C}}{\partial a_{i j}} - κ \frac{\partial {\dot{a}}_{m n}^{I R D}}{\partial a_{i j}} + \frac{\partial {\dot{a}}_{m n}^{A R D}}{\partial a_{i j}} + \frac{\partial {\dot{a}}_{m n}^{Δ R S C}}{\partial a_{i j}}$

where

(77) $\frac{\partial {\dot{a}}_{m n}^{Δ R S C}}{\partial a_{i j}} = - 2 \dot{γ} (1 - κ) \{\begin{matrix} [\frac{\partial M_{m n k l}}{\partial a_{i j}} - δ_{k l} δ_{m i} δ_{n j} - 5 \frac{\partial}{\partial a_{i j}} \{L_{m n k l} - M_{m n r s} a_{r s k l}\}] C_{k l} \\ + [M_{m n k l} - δ_{k l} a_{m n} - 5 (L_{m n k l} - M_{m n r s} a_{r s k l})] \frac{\partial C_{k l}}{\partial a_{i j}} \end{matrix}\}$

The partial derivative terms on the right-hand side (R.H.S) of Equation (90) appear in the previous sections.

2.5. Closure Approximations and Their Explicit Derivatives

Derivatives of the orientation tensor closure approximation are used in the Newton-Raphson iteration method to compute the steady-state fiber orientation tensor state. The derivatives used in this study are described below.

2.5.1. Quadratic Closure Approximation

The quadratic (QDR) closure ${\tilde{a}}_{i j k l}$ was introduced by Doi [51] and Lipscomb [52] and is defined as dyadic product of the 2nd order orientation tensor $a_{i j}$ :

(78) ${\tilde{a}}_{i j k l} = a_{i j} a_{k l}$

The derivative of

{\tilde{a}}_{i j k l}

with respect to

a_{m n}

(79) $\frac{\partial {\tilde{a}}_{i j k l}}{\partial a_{m n}} = \frac{\partial a_{i j}}{\partial a_{m n}} a_{k l} + a_{i j} \frac{\partial a_{k l}}{\partial a_{m n}} = δ_{i m} δ_{j n} a_{k l} + a_{i j} δ_{k m} δ_{l n}$

The QDR closure inherently lacks symmetry but preserves the required symmetry of the computed

a_{m n}

when employed in the Advani-Tucker equation (cf. Equation (18)).

2.5.2. Linear Closure Approximation

The linear (LIN) closure approximation of ${\hat{a}}_{i j k l}$ first proposed by Hand [17] uses all possible products of the 2nd order orientation tensor $a_{i j}$ and the identity tensor $δ_{i j}$ as follows:

(80) $\begin{matrix} {\hat{a}}_{i j k l} = - h_{1} (δ_{i j} δ_{k l} & + δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) \\ + h_{2} (a_{i j} δ_{k l} + a_{i k} δ_{j l} + a_{i l} δ_{j k} + δ_{i j} a_{k l} + δ_{i k} a_{j l} + δ_{i l} a_{j k}) \end{matrix}$

where

h_{1}

and

h_{2}

are numerical factors which vary based on spatial dimensionality.

h_{1} = 1 / 35

for three-dimensional orientation and

h_{1} = 1 / 24

for planar orientation, while

h_{2} = 1 / 7

for three-dimensional orientation and

h_{2} = 1 / 6

for planar orientation [7]. It follows that the derivative of

{\hat{a}}_{i j k l}

with respect to

a_{m n}

(81) $\frac{\partial {\hat{a}}_{i j k l}}{\partial a_{m n}} = h_{2} (δ_{i m} δ_{j n} δ_{k l} + δ_{i m} δ_{k n} δ_{j l} + δ_{i m} δ_{l n} δ_{j k} + δ_{i m} δ_{j n} δ_{k l} + δ_{i m} δ_{k n} δ_{j l} + δ_{i m} δ_{l n} δ_{j k})$

The LIN closure is exact for random orientation distribution, while the QDR closures are exact for a uniaxially aligned fiber orientation.

2.5.3. Hybrid Closure Approximation

The hybrid closure approximation $a_{i j k l}$ is a weighted combination of the ${\hat{a}}_{i j k l}$ and ${\tilde{a}}_{i j k l}$ defined above, written in terms of Herman’s Orientation factor (HOF), $f$ , as follows [7]:

(82) $a_{i j k l} = f {\tilde{a}}_{i j k l} + (1 - f) {\hat{a}}_{i j k l}$

Advani & Tucker [7] proposed an expression for

f

as an invariant of the orientation state given as

f = a_{f} a_{i j} a_{j i} - b_{f}

, where

a_{f}

and

b_{f}

are constants that depend on the spatial dimension.

a_{f} = 3 / 2

for three-dimensional orientation and

a_{f} = 2

for planar orientation, while

b_{f} = 1 / 2

for three-dimensional orientation and

b_{f} = 1

for planar orientation.

The derivative of the hybrid closure approximation $a_{i j k l}$ with respect to components of the 2nd order tensor $a_{m n}$ is

(83) $\frac{\partial a_{i j k l}}{\partial a_{m n}} = f \frac{\partial {\tilde{a}}_{i j k l}}{\partial a_{m n}} + (1 - f) \frac{\partial {\hat{a}}_{i j k l}}{\partial a_{m n}} + \frac{\partial f}{\partial a_{m n}} ({\tilde{a}}_{i j k l} - {\hat{a}}_{i j k l})$

where

(84) $\frac{\partial f}{\partial a_{m n}} = a_{f} (δ_{i m} δ_{j n} a_{j i} + a_{i j} δ_{j m} δ_{i n})$

We denote the hybrid closure based on this HOF form as HYB₁. An alternative form of the HOF,

f

, proposed by Advani & Tucker [7], is

(85) $f = 1 - α^{α} e_{i j k} a_{i 1} a_{j 2} a_{k 3}$

where it follows that

$\frac{\partial f}{\partial a_{m n}} = - α^{α} e_{i j k} \{δ_{i m} δ_{1 n} a_{j 2} a_{k 3} + a_{i 1} δ_{j m} δ_{2 n} a_{k 3} + a_{i 1} a_{j 2} δ_{k m} δ_{3 n}\}$

We denote hybrid closure based on this alternative form of the HOF as HYB₂. The hybrid model is observed to perform better during the transient orientation prediction; however, the hybrid closure was shown to over-predict the steady state fiber alignment compared with simulation results based on the more accurate orientation distribution function, which are computationally expensive since they require a finite difference grid in space and time.

2.5.4. Hinch and Leal Closure Approximation

Hinch and Leal [33] developed numerous composite closure approximations for the 4th order orientation tensor in precontracted forms with the deformation rate tensor. The accuracy of their predictions was shown to depend on the flow field and flow magnitude considered. Hinch and Leal (HL) closure approximations are not explicit expressions of the 4th order orientation tensor $a_{i j k l}$ but are instead written in contracted form with the deformation rate tensor, i.e., $γ_{k l} a_{i j k l}$ . Advani and Tucker provide a general explicit expression of $a_{i j k l}$ (cf. Equation (86)), summarizing all Hinch and Leal (HL) closures as

(86) $\begin{matrix} a_{i j k l} = β_{1} (δ_{i j} δ_{k l}) & + β_{2} (δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) + β_{3} (δ_{i j} a_{k l} + a_{i j} δ_{k l}) \\ + β_{4} (a_{i k} δ_{j l} + a_{j l} δ_{i k} + a_{i l} δ_{j k} + a_{j k} δ_{i l}) + β_{5} (a_{i j} a_{k l}) + β_{6} (a_{i k} a_{j l} + a_{i l} a_{j k}) \\ + \dots + β_{7} (δ_{i j} a_{k m} a_{m l} + a_{i m} a_{m j} δ_{k l}) + β_{8} (a_{i m} a_{m j} a_{k n} a_{n l}) \end{matrix}$

where

β_{i}

are defined in Table 1. The partial derivative of the above expression with respect to

a_{r s}

(87) $\begin{matrix} \frac{\partial a_{i j k l}}{\partial a_{r s}} = [\frac{\partial β_{1}}{\partial a_{r s}} (δ_{i j} δ_{k l}) & + \frac{\partial β_{2}}{\partial a_{r s}} (δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) + \frac{\partial β_{3}}{\partial a_{r s}} (δ_{i j} a_{k l} + a_{i j} δ_{k l}) + \frac{\partial β_{4}}{\partial a_{r s}} (a_{i k} δ_{j l} + a_{j l} δ_{i k} + a_{i l} δ_{j k} + a_{j k} δ_{i l}) \\ + \frac{\partial β_{5}}{\partial a_{r s}} (a_{i j} a_{k l}) + \frac{\partial β_{6}}{\partial a_{r s}} (a_{i k} a_{j l} + a_{i l} a_{j k}) + \frac{\partial β_{7}}{\partial a_{r s}} (δ_{i j} a_{k m} a_{m l} + a_{i m} a_{m j} δ_{k l}) \\ + \frac{\partial β_{8}}{\partial a_{r s}} (a_{i m} a_{m j} a_{k n} a_{n l})] \\ + [β_{3} (δ_{i j} δ_{k r} δ_{l s} + δ_{i r} δ_{j s} δ_{k l}) + β_{4} (δ_{i r} δ_{k s} δ_{j l} + δ_{j r} δ_{l s} δ_{i k} + δ_{i r} δ_{l s} δ_{j k} + δ_{j r} δ_{k s} δ_{i l}) \\ + β_{5} (δ_{i r} δ_{j s} a_{k l} + a_{i j} δ_{k r} δ_{l s}) + β_{6} (δ_{i r} δ_{k s} a_{j l} + a_{i k} δ_{j r} δ_{l s} + δ_{i r} δ_{l s} a_{j k} + a_{i l} δ_{j r} δ_{k s}) \\ + β_{7} (δ_{i j} δ_{k r} δ_{m s} a_{m l} + δ_{i j} a_{k m} δ_{m r} δ_{l s} + δ_{i r} δ_{m s} a_{m j} δ_{k l} + a_{i m} δ_{m r} δ_{j s} δ_{k l}) \\ + β_{8} (δ_{i r} δ_{m s} a_{m j} a_{k n} a_{n l} + a_{i m} δ_{m r} δ_{j s} a_{k n} a_{n l} + a_{i m} a_{m j} δ_{k r} δ_{n s} a_{n l} + a_{i m} a_{m j} a_{k n} δ_{n r} δ_{l s})] \end{matrix}$

Mullens [40] provided a summary table (cf. Table 1) for the

β_{i}

factors of the Hinch and Leal (HL) closures subdivided into weak flow (WF) (Isotropic (ISO), Linear (LIN) and Quadratic (QDR)) and strong flow (SF), and Hinch and Leal first composite (HL1) and second composite (HL2) closure forms.

Table 1

Summary of the Hinch and Leal closure $β_{i}$ factors for the different flow classifications.

		$β_{1}$	$β_{2}$	$β_{3}$	$β_{4}$	$β_{5}$	$β_{6}$	$β_{7}$	$β_{8}$
WF	ISO	$\frac{1}{15}$	$\frac{1}{15}$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
	LIN	$- \frac{1}{35}$	$- \frac{1}{35}$	$\frac{1}{7}$	$\frac{1}{7}$	$\dots$	$\dots$	$\dots$	$\dots$
	QDR	$\dots$	$\dots$	$\dots$	$\dots$	$1$	$\dots$	$\dots$	$\dots$
SF	SF2	$\dots$	$\dots$	$\dots$	$\dots$	$1$	$1$	$\dots$	$- \frac{2}{〈a^{2}〉}$
HL	HL1	$\dots$	$\dots$	$\frac{2}{5}$	$\dots$	$- \frac{1}{5}$	$\frac{3}{5}$	$- \frac{2}{5}$	$\dots$
HL	HL2	$\frac{26}{315} α$	$\frac{26}{315} α$	$\frac{16}{63} α$	$- \frac{4}{21} α$	$1$	$1$	$\dots$	$- \frac{2}{〈a^{2}〉}$

where the parameters $〈a^{2}〉$ and $α$ are, respectively,

(88) $〈a^{2}〉 = a_{i j} a_{j i}, α = \exp [2 \frac{1 - 3 〈a^{2}〉}{1 - 〈a^{2}〉}]$

and the partial derivatives are, respectively, given as

(89) $\frac{\partial 〈a^{2}〉}{\partial a_{r s}} = δ_{i r} δ_{j s} a_{j i} + a_{i j} δ_{j r} δ_{i s} \frac{\partial α}{\partial a_{r s}} = - \frac{4 α}{{(1 - 〈a^{2}〉)}^{2}} \frac{\partial 〈a^{2}〉}{\partial a_{r s}}, \frac{\partial}{\partial a_{r s}} \{\frac{k}{〈a^{2}〉}\} = - \frac{k}{{〈a^{2}〉}^{2}} \frac{\partial 〈a^{2}〉}{\partial a_{r s}}$

2.5.5. Eigenvalue-Based Fitted (EBF) Closure Approximations

More recently, fitted closure approximations have emerged, including the Eigenvalue-Based Fitted (EBF) closures that require a principal axis transformation which ensures material objectivity. Orthotropic closure approximations for the 4th order tensor implicitly impose objectivity such that the approximation is independent of the coordinate frame selection. Indeed, the Orthotropic Linear (ORT), Orthotropic Smooth (ORS), Orthotropic Fitted (ORF) closures, ORF closure for low fiber-fiber interaction coefficient (ORL), ORF closure for wide range of fiber-fiber interaction coefficient values with 2nd order orientation tensor eigenvalue polynomial expansions (ORW), ORW closure with 3^rd order orientation tensor eigenvalue polynomial expansions (ORW3), orthotropic natural closure—exact midpoint fit (NAT₁) and orthotropic natural closure—and extended quadratic fit (NAT₂) are included in the EBF family of closures originally developed by Cintra and Tucker [39] based on the assumption of coincident principal axis of the 2nd and 4th order tensor. The ORF has been shown to have better performance compared to non-fitted closure approximations; however, the ORF behaved poorly for flows with very low interaction coefficients and sometimes gave non-physical oscillations similar to that seen with the Hinch and Leal closure (HL2) under the same conditions. The ORL behaves better with a low interaction coefficient in simple shear flow but overpredicts flow direction alignment and is unstable for radial diverging flows. Due to the inherent symmetries of $a_{i j k l}$ , its components can be represented in (6 × 6) contracted notation as

(90) $A_{r s} = a_{i j k l}$

where the index of the contracted notation is related to the index notation according to

(91) $r = \{\begin{matrix} i = j & {f o r, δ}_{i j} = 1 \\ (9 - i - j) & {f o r, δ}_{i j} = 0 \end{matrix} and s = \{\begin{matrix} k = l & {f o r, δ}_{k l} = 1 \\ (9 - k - l) & f o r, δ_{k l} = 0 \end{matrix}$

It follows that the derivative of the 4th order tensor with respect to

a_{m n}

in contracted notation is

(92) $\frac{\partial A_{r s}}{\partial a_{m n}} = \frac{\partial a_{i j k l}}{\partial a_{m n}}$

Symmetry of the 4th order tensor requires

a_{i j k l} = a_{k l i j}

, which implies that

A_{r s} = A_{s r}

. The contracted tensor

A_{r s}

transformed to the principal axes has the orthotropic form

{\bar{A}}_{r s}

given as

(93) $\bar{\underline{\underline{A}}} = [\begin{matrix} {\bar{A}}_{11} & {\bar{A}}_{12} & {\bar{A}}_{13} \\ {\bar{A}}_{21} & {\bar{A}}_{22} & {\bar{A}}_{23} \\ {\bar{A}}_{31} & {\bar{A}}_{32} & {\bar{A}}_{33} \\ {\bar{A}}_{44} \\ {\bar{A}}_{55} \\ {\bar{A}}_{66} \end{matrix}]$

The contracted tensor transforms from its principal reference frame to the original coordinate axes according to

(94) $A_{r s} = M_{r i} M_{s j} {\bar{A}}_{i j}$

The 6 × 6 transformation matrix

M_{i j}

is given as

M_{i j} = F_{i m} Q_{m n} F_{n j}^{- 1}

, where

F_{i j} = k δ_{i j}, k = \{\begin{matrix} \begin{matrix} 1, & i \leq 3 \end{matrix} \\ \begin{matrix} 2, & i > 3 \end{matrix} \end{matrix}

, and

Q_{m n}

is the contracted form of

Q_{i j k l}

, where

Q_{i j k l} = Φ_{i k} Φ_{j l} + (1 - δ_{k l}) Φ_{j k} Φ_{i l}

. The modal matrix

Φ_{i j}

whose kth column is the corresponding eigenvector

{\underline{x}}^{k}

of eigenvalues

λ_{k} = Λ_{k k}

(no summation implied) is obtained from the spectral decomposition of

a_{i j}

(95) $\underline{\underline{Φ}} | a_{m n} = Φ_{m k} Λ_{k l} Φ_{n l}$

A more direct way is to transform the 4th order orientation tensor given in its principal reference frame as

{\bar{a}}_{m n p q}

to that in the original reference frame

a_{i j k l}

according to

(96) $a_{i j k l} = Φ_{i m} Φ_{j n} Φ_{k p} Φ_{l q} {\bar{a}}_{m n p q}$

where

{\bar{a}}_{m n p q}

is reconstructed from the contracted form

{\bar{A}}_{r s} .

Using the product rule

(97) $\begin{matrix} \frac{\partial a_{i j k l}}{\partial a_{r s}} = Φ_{i m} Φ_{j n} & Φ_{k p} Φ_{l q} \frac{\partial {\bar{a}}_{m n p q}}{\partial a_{r s}} + \dots \\ + (\frac{\partial Φ_{i m}}{\partial a_{r s}} Φ_{j n} Φ_{k p} Φ_{l q} + Φ_{i m} \frac{\partial Φ_{j n}}{\partial a_{r s}} Φ_{k p} Φ_{l q} + Φ_{i m} Φ_{j n} \frac{\partial Φ_{k p}}{\partial a_{r s}} Φ_{l q} \\ + Φ_{i m} Φ_{j n} Φ_{k p} \frac{\partial Φ_{l q}}{\partial a_{r s}}) {\bar{a}}_{m n p q} \end{matrix}$

The derivative of the eigenvector tensor

\underline{\underline{Φ}}

appears in [49] (cf. Appendix A). Symmetry requirements of the transformed orthotropic tensor reduces the total number of independent non-zero components to 9, and additional special symmetry properties of the exact 4th order tensor require that

{\bar{a}}_{i j k l} = {\bar{a}}_{k j i l} = {\bar{a}}_{l j k i} = {\bar{a}}_{i k j l} = {\bar{a}}_{i l k j}

, reducing the non-zero independent components in

A_{i j}

to the 6 diagonal terms, i.e.,

(98) ${\bar{A}}_{i j} = {\bar{A}}_{k k} k = 9 - i - j, i \neq j$

where repeated indices for

{\bar{A}}_{k k}

do not imply summation. The normalization property

a_{i j k k} = a_{i j}

of the exact 4th order tensor further requires that

(99) $[\begin{matrix} {\bar{A}}_{44} \\ {\bar{A}}_{55} \\ {\bar{A}}_{66} \end{matrix}] = {\underline{\underline{ẞ}}}^{- 1} \{[\begin{matrix} λ_{1} \\ λ_{2} \\ λ_{3} \end{matrix}] - [\begin{matrix} {\bar{A}}_{11} \\ {\bar{A}}_{22} \\ {\bar{A}}_{33} \end{matrix}]\}$

where

λ_{i}

are the eigenvalues of

a_{i j}

, with

\sum_{i} λ_{i} = 1

and

ẞ_{i j} = 1 - δ_{i j}

. Based on the foregoing conditions, the only three surviving non-zero independent terms are

{\bar{A}}_{11}

{\bar{A}}_{22}

and

{\bar{A}}_{33}

. The general form for orthotropic closure is to express the three surviving non-zero independent components (

{\bar{A}}_{11}

{\bar{A}}_{22},

{\bar{A}}_{33}

) of the contracted 4th order tensor in the principal reference frame after imposing all symmetric and normalization conditions of the exact 4th order tensor, as a scalar function

F_{k} (λ_{1}, λ_{2})

of the two largest eigenvalues (

λ_{1}, λ_{2}

) of

a_{i j}

. Most fitted closures take the form of an nth-order binomial function in

λ_{1}

and

λ_{2}

(note that

λ_{3} = 1 - λ_{1} - λ_{2})

to represent

F_{k}

(100) ${\bar{A}}_{k k} = F_{k} (λ_{1}, λ_{2}) = f_{k}^{(n)} (λ_{1}, λ_{2}), λ_{1} \geq λ_{2} \geq λ_{3}, k = 1, 2, 3$

where, again, no summation is implied by the repeated indices for

{\bar{A}}_{k k}

. Polynomial order exceeding

n = 4

falls under the class of eigenvalue-based optimal fitting (EBOF) closures. The Verweyst (VST) [53] closure and some members of the family of Least Absolute Residuals (LAR) regression closures developed by Mullens M. [40], such as the 4th order optimal fitting LAR closure (LAR4) and Flow-Field LAR4 closure (FFLAR4), fall under the EBOF category. Generally, we can represent the function

f_{k}^{(n)}

as a tensor product of a constant coefficient matrix

C_{i j}^{(n)}

and an nth order permuted bivariate polynomial vector

⩓_{j}^{(n)} = ⩓_{j}^{(n)} (λ_{1}, λ_{2})

, i.e.,

(101) $f_{k}^{(n)} (λ_{1}, λ_{2}) = C_{k j}^{(n)} ⩓_{j}^{(n)}$

$k = 1,2, 3, j = 1 \dots (n + 1) (n + 2) / 2$

Representations of

C_{i j}^{(n)}

and

⩓_{j}^{(n)}

for various closure approximations appear in Appendix B. The derivative of the components of the orthotropic closure in contracted notation with respect to

a_{r s}

(102) $\frac{\partial {\bar{A}}_{k k}}{\partial a_{r s}} = \frac{\partial {\bar{A}}_{k}}{\partial a_{r s}} = C_{k j}^{(n)} \frac{\partial ⩓_{j}^{(n)}}{\partial a_{r s}} = C_{k j}^{(n)} ⩓_{j r s}^{' (n)} = C_{k j}^{(n)} {\tilde{⩓}}_{j l}^{(n)} λ_{l r s}^{'} k = 1,2, 3, l = 1,2$

The nth order binomial permutation vector

⩓_{k}^{(n)}

and its derivative coefficient matrix

{\tilde{⩓}}_{k l}^{(n)}

for the quadratic closure are given from terms of binomial expansion, respectively, as

(103) $\{\begin{matrix} \begin{matrix} ⩓_{k}^{(n)} (λ_{1}, λ_{2}) = λ_{1}^{i - j} λ_{2}^{j} \\ {\tilde{⩓}}_{k l}^{(n)} = \frac{\partial}{\partial λ_{l}} ⩓_{k}^{(n)} = \{\begin{matrix} (i - j) \cdot λ_{1}^{i - j - 1} λ_{2}^{j} & l = 1 \\ j \cdot λ_{1}^{i - j} λ_{2}^{j - 1} & l = 2 \end{matrix} \end{matrix}, k = 1 + j + i (i + 1) / 2 \end{matrix} j = 0 \dots i, i = 0 \dots n$

For a special case of orthotropic fitted closure called rational elliptical (RE) closure, such as the Wetzel (WTZ) RE closure [54] and the Rational LAR closure with 3^rd order numerator and 2nd order denominator optimal fitting (LAR32) [40], the scalar function for the 3 independent tensor components is

(104) $F_{k} (λ_{1}, λ_{2}) = \frac{f_{k}^{(n)} (λ_{1}, λ_{2})}{f_{k}^{(m)} (λ_{1}, λ_{2})}$

The derivative of the components of the above with respect to

a_{r s}

(105) $\frac{\partial {\bar{A}}_{k k}}{\partial a_{r s}} = \frac{1}{{[f^{(m)}]}^{2}} [f^{(m)} \frac{\partial f^{(n)}}{\partial a_{r s}} - f^{(n)} \frac{\partial f^{(m)}}{\partial a_{r s}}]$

Again, no summation is implied by the repeated indices in

{\bar{A}}_{k k}

. From the normalization condition of the 4th order tensor, we obtain the following for the derivative of

{\bar{A}}_{k k}

(k = 4,5, 6)

(106) $\frac{\partial {\bar{A}}_{k k}}{\partial a_{r s}} = ẞ_{k i}^{- 1} \{\frac{\partial λ_{i}}{\partial a_{r s}} - \frac{\partial {\bar{A}}_{i i}}{\partial a_{r s}}\}, \frac{\partial λ_{i}}{\partial a_{r s}} = E_{i l} λ_{l r s}^{'}, \underline{\underline{E}} = {[\begin{matrix} 1 & 0 & - 1 \\ 0 & 1 & - 1 \end{matrix}]}^{T} i = 1,2, 3, l = 1,2$

For the partial derivative of the eigenvalues with respect to the components of the 2nd order orientation tensor, kindly refer to Appendix A. EBF closures are computationally more expensive because of the principal axis transformation.

2.5.6. Invariant-Based Fitted (IBF) Closure Approximations

Of the class of IBF closures, the natural (NAT) closure approximation of Verleye et al. [55] was built on the work of Lipscomb et al. [52] and formed the basis for other IBF developments. Verleye et al. [46] developed a general expression for the fully symmetric 4th order tensor in terms of invariants of the 2nd order tensor, the identity tensor and fitted coefficients which were determined through a least square fitting process. The NAT closure assumed the absence of fiber-fiber interaction and infinitely long fiber geometry. The closure is based on the foregoing assumptions; however, it has been reported to possess singularities for axisymmetric orientation states. The Invariant-based optimal fitting (IBOF) closure developed by Chung and Kwon [34] was an extension to the NAT closure development; however, the independent coefficients are derived from regression analysis based on actual data obtained from distribution function calculation (DFC) for multiple flow types (similar to the procedure used for EBF closures). In the contracted form, the 4th order tensor $a_{i j k l}$ based on symmetry properties is given as

(107) $\underline{\underline{A}} = [\begin{matrix} A_{11} & A_{12} & A_{13} & A_{14} & A_{15} & A_{16} \\ A_{22} & A_{23} & A_{24} & A_{24} & A_{26} \\ A_{33} & A_{34} & A_{35} & A_{36} \\ A_{44} & A_{45} & A_{46} \\ A_{55} & A_{56} \\ \dots S y m & A_{66} \end{matrix}]$

which conforms to the symmetry requirement

(108) $A_{44} = A_{23}, A_{45} = A_{36}, A_{46} = A_{25} A_{55} = A_{13}, A_{56} = A_{14}, A_{66} = A_{12}$

and normalization condition

(109) $\sum_{n = 1}^{3} A_{n m} = a_{m}, a_{m} = a_{i j}, m = \{\begin{matrix} i = j & i = j \\ 9 - i - j & i \neq j \end{matrix}$

It follows that under these conditions, the relationship between

a_{i j k l}

and

a_{i j}

can be written as

(110) $A_{11} + A_{12} + A_{13} = a_{11}, A_{12} + A_{22} + A_{23} = a_{22} A_{13} + A_{23} + A_{33} = a_{33}, A_{14} + A_{24} + A_{34} = a_{23} A_{15} + A_{25} + A_{35} = a_{13}, A_{16} + A_{26} + A_{36} = a_{12}$

Taking partial derivatives of Equations (108) and (109), we obtain the following:

(111) $\frac{\partial A_{m n}}{\partial a_{r s}} = \frac{\partial A_{i j}}{\partial a_{r s}}, a n d \sum_{n = 1}^{3} \frac{\partial A_{n m}}{\partial a_{r s}} = \frac{\partial a_{m}}{\partial a_{r s}}$

There are thus only 9 independent components for the 4th order tensor. The IBOF is developed in terms of the fully symmetric 4th order expansion of

a_{i j k l}

as a combination of the 2nd order tensor

a_{i j}

and identity tensor

δ_{k l}

, where the Cayley-Hamilton theory is used to obtain

(112) $a_{i j k l} = β_{1} S (δ_{i j} δ_{k l}) + β_{2} S (δ_{i j} a_{k l}) + β_{3} S (a_{i j} a_{k l}) + β_{4} S (δ_{i j} a_{k m} a_{m l}) + β_{5} S (a_{i j} a_{k m} a_{m l}) + β_{6} S (a_{i m} a_{m j} a_{k n} a_{n l})$

In the above, the

S

operator represents the symmetric permutation expansion of its argument, for example,

(113) $\begin{matrix} S (T_{i j k l}) = \frac{1}{24} [T_{i j k l} & + T_{i j l k} + T_{i k j l} + T_{i k l j} + T_{i l j k} + T_{i l k j} + T_{j i k l} + T_{j i l k} + T_{j k i l} + T_{j k l i} + T_{j l i k} + T_{j l k i} + T_{k i j l} \\ + T_{k i l j} + T_{k j i l} + T_{k j l i} + T_{k l i j} + + T_{k l j i} + T_{l i j k} + T_{l i k j} + T_{l j i k} + T_{l j k i} + T_{l k i j} + T_{l k j i}] \end{matrix}$

To incorporate the IBOF closure in our NR solution approach, we differentiate

a_{i j k l}

with respect to

a_{r s}

to obtain

(114) $\begin{matrix} \frac{\partial}{\partial a_{r s}} \{a_{i j k l}\} = & [\frac{\partial β_{1}}{\partial a_{r s}} S (δ_{i j} δ_{k l}) + \frac{\partial β_{2}}{\partial a_{r s}} S (δ_{i j} a_{k l}) + \frac{\partial β_{3}}{\partial a_{r s}} S (a_{i j} a_{k l}) + \frac{\partial β_{4}}{\partial a_{r s}} S (δ_{i j} a_{k m} a_{m l}) \\ + \frac{\partial β_{5}}{\partial a_{r s}} S (a_{i j} a_{k m} a_{m l}) + \frac{\partial β_{6}}{\partial a_{r s}} S (a_{i m} a_{m j} a_{k n} a_{n l})] + \dots \\ + [β_{2} S (δ_{i j} δ_{k r} δ_{l s}) + β_{3} \{S (δ_{i r} δ_{j s} a_{k l}) + S (a_{i j} δ_{k r} δ_{l s})\} \\ + β_{4} \{S (δ_{i j} δ_{k r} δ_{m s} a_{m l}) + S (δ_{i j} a_{k m} δ_{m r} δ_{l s})\} \\ + β_{5} \{S (δ_{i r} δ_{j s} a_{k m} a_{m l}) + S (a_{i j} δ_{k r} δ_{m s} a_{m l}) + S (a_{i j} a_{k m} δ_{m r} δ_{l s})\} \\ + β_{6} {S (δ_{i r} δ_{m s} a_{m j} a_{k n} a_{n l}) + S (a_{i m} δ_{m r} δ_{j s} a_{k n} a_{n l}) \\ + S (a_{i m} a_{m j} δ_{k r} δ_{n s} a_{n l}) + S (a_{i m} a_{m j} a_{k n} δ_{n r} δ_{l s})}] \end{matrix}$

The

β_{i}

coefficients are expressed as functions of the second and third invariants (

I I a n d I I I, r e s p e c t i v e l y

) of

a_{i j}

. Employing Equations (111) and (112), and applying the Cayley-Hamilton theorem, it can be shown that only 3 independent coefficients remain. It follows that the expressions for the IBOF dependent coefficients

(β_{1}, β_{2}, β_{5})

are

(115) $β_{1} = \frac{3}{5} [- \frac{1}{7} + \frac{1}{5} β_{3} (\frac{1}{7} + \frac{4}{7} I I + \frac{8}{3} I I I) - β_{4} (\frac{1}{5} - \frac{8}{15} I I - \frac{14}{15} I I I) - β_{6} (\frac{1}{35} - \frac{4}{35} I I - \frac{24}{105} I I I + \frac{16}{15} I I I I I + \frac{8}{35} {I I}^{2})] β_{2} = \frac{6}{7} [1 - \frac{1}{5} β_{3} (1 + 4 I I) + \frac{7}{5} β_{4} (\frac{1}{6} - I I) - β_{6} (- \frac{1}{5} + \frac{4}{5} I I + \frac{2}{3} I I I - \frac{8}{5} {I I}^{2})] β_{5} = - \frac{4}{5} β_{3} - \frac{7}{5} β_{4} - \frac{6}{5} β_{6} (1 - \frac{4}{3} I I)$

which can be differentiated with respect to

a_{r s}

to yield

(116) $\begin{matrix} \frac{\partial β_{1}}{\partial a_{r s}} = \frac{3}{5} [\frac{1}{5} \frac{\partial β_{3}}{\partial a_{r s}} & (\frac{1}{7} + \frac{4}{7} I I + \frac{8}{3} I I I) - \frac{\partial β_{4}}{\partial a_{r s}} (\frac{1}{5} - \frac{8}{15} I I - \frac{14}{15} I I I) + \dots \\ - \frac{\partial β_{6}}{\partial a_{r s}} (\frac{1}{35} - \frac{4}{35} I I - \frac{24}{105} I I I + \frac{16}{15} I I I I I + \frac{8}{35} {I I}^{2})] + \dots \\ + \frac{3}{5} [[\frac{4}{35} β_{3} - β_{4} + β_{6} (\frac{4}{35} - \frac{16}{35} I I - \frac{16}{15} I I I)] \frac{\partial I I}{\partial a_{r s}} + [\frac{1}{5} \frac{8}{3} β_{3} + \frac{14}{15} β_{4} + β_{6} (\frac{24}{105} - \frac{16}{15} I I)] \frac{\partial I I I}{\partial a_{r s}}] \end{matrix} \begin{matrix} \frac{\partial β_{2}}{\partial a_{r s}} = \frac{6}{7} [- \frac{1}{5} \frac{\partial β_{3}}{\partial a_{r s}} (1 + 4 I I) + \frac{7}{5} \frac{\partial β_{4}}{\partial a_{r s}} (\frac{1}{6} - I I) - \frac{\partial β_{6}}{\partial a_{r s}} (- \frac{1}{5} + \frac{4}{5} I I + \frac{2}{3} I I I - \frac{8}{5} {I I}^{2})] \\ + \frac{6}{7} [- \frac{1}{5} [4 β_{3} + 7 β_{4} + β_{6} (4 - 16 I I)] \frac{\partial I I}{\partial a_{r s}} - \frac{2}{3} β_{6} \frac{\partial I I I}{\partial a_{r s}}] \end{matrix} \frac{\partial β_{5}}{\partial a_{r s}} = - \frac{4}{5} \frac{\partial β_{3}}{\partial a_{r s}} - \frac{7}{5} \frac{\partial β_{4}}{\partial a_{r s}} - \frac{6}{5} \frac{\partial β_{6}}{\partial a_{r s}} (1 - \frac{4}{3} I I) + \frac{8}{5} β_{6} \frac{\partial I I}{\partial a_{r s}}$

The independent coefficients

(β_{3}, β_{4}, β_{6})

in the study by Chung et al. [34] are obtained from a 5th order binomial fitted function in terms of

I I a n d I I I

(117) $β_{m} = \sum_{i = 0}^{5} \sum_{j = 0}^{i} a_{k}^{m} \cdot {I I}^{i - j} {I I I}^{j}, k = j + i (i + 1) / 2, m = 1, 2, 3$

where the coefficients of the binomial terms appear in Table A1. The non-unity invariants of

a_{2}

are given as

(118) $I I = λ_{1} λ_{2} + λ_{2} λ_{3} + λ_{3} λ_{1}, I I I = λ_{1} λ_{2} λ_{3}$

It follows that derivative of the independent coefficients in Equation (120) with respect to the components

a_{r s}

are

(119) $\frac{\partial β_{m}}{\partial a_{r s}} = \sum_{i = 0}^{5} \sum_{j = 0}^{i} a_{k}^{m} \{(i - j) \cdot {I I}^{i - j - 1} {I I I}^{j} \frac{\partial I I}{\partial a_{r s}} + j \cdot {I I}^{i - j - 1} {I I I}^{j - 1} \frac{\partial I I I}{\partial a_{r s}}\}$

where

(120) $\frac{\partial I I}{\partial a_{r s}} = (λ_{2} + λ_{3}) \frac{\partial λ_{1}}{\partial a_{r s}} + (λ_{1} + λ_{3}) \frac{\partial λ_{2}}{\partial a_{r s}} + (λ_{1} + λ_{2}) \frac{\partial λ_{3}}{\partial a_{r s}}$

$\frac{\partial I I I}{\partial a_{r s}} = (λ_{2} λ_{3}) \frac{\partial λ_{1}}{\partial a_{r s}} + (λ_{1} λ_{3}) \frac{\partial λ_{2}}{\partial a_{r s}} + (λ_{1} λ_{2}) \frac{\partial λ_{3}}{\partial a_{r s}}$

2.6. Error Estimate

In the numerical examples to follow, we assess the performance of the Newton-Raphson (NR) method to accurately predict the steady-state values of $a_{i j}$ based on the relative absolute error between results of our NR approach and that of a reference method, in this case the explicit, single-step, Dormand-Prince 4th and 5th order Runge–Kutta adaptive-stepping method (RK45) with a set relative and absolute tolerance of $10^{- 12}$ , which ensures the solutions are accurate up to 12 digits. We define the percentage error as

(121) ${e r r}_{1} = \frac{a_{m n}^{N R} - a_{m n}^{R K 4}}{a_{m n}^{R K 4}} \times 100 %$

The NR iteration is said to converge to a unique solution, and the iteration is terminated once the NR error tolerance

{e r r}_{N R} = {‖{\underline{R}}_{i j}^{N R}‖}_{2}

falls below some predetermined tolerance value,

{e r r}_{d} = 10^{- 6}

, where

{\underline{R}}_{i j}^{N R}

is given by Equation (7).

3. Results

Example calculations are presented below to illustrate and validate our NR-based approach for computing steady-state values of $a_{i j}$ under various selected flow conditions. Prior to providing these examples, we first validate the explicit analytical expressions for the Jacobian introduced in Equation (4) for the various fiber orientation tensor evolution models and closures presented above. Results computed at steady state using our Newton Raphson approach are validated with computed results obtained using MATLAB’s 2025a inbuilt explicit, adaptive-stepping, Dormand-Prince 4th and 5th order, Runge–Kutta ode45 (RK45) solver (MathWorks, Natick, MA, USA). We note that the RK45 basis for validating our NR results are likewise numerical approximations, and as such, the relative error metric provided only serves as an indication of the agreement between these two outcomes. The validation of the NR approach is performed using three different orientation tensor model case studies which include (a) an anisotropic rotary diffusion (ARD) tensor model, (b) a slow orientation kinetic (SOK) tensor model, and (c) a composite/coupled tensor model. Subsequent sections are dedicated to assessing the performance of the NR method using various 4th order orientation tensor closure approximations and various homogenous flows. Besides the relative error tolerance, other metrics for assessing the NR method’s performance in the example calculations below include the rate of numerical convergence and the number of NR iterations to converge.

3.1. Validation of Exact Jacobians Based on Finite Difference Approximation

Exact Jacobians ${\underline{\underline{J}}}^{e x a c t}$ defined through the exact derivative expressions provided in Section 2 are validated here using approximate values of the Jacobian ${\underline{\underline{J}}}^{F D}$ calculated with the $O (h^{4})$ central finite difference method. The validity of our computed exact Jacobians is assessed here based on the error defined by the Euclidean norm, given as

(122) ${e r r}_{2} = {‖{\underline{\underline{J}}}^{e x a c t} - {\underline{\underline{J}}}^{F D}‖}_{2}$

Values of

{\underline{\underline{J}}}^{F D}

are computed in this work using the higher order central difference finite difference approximation

{C T}_{4}

written as

(123) $J_{m n i j}^{F D} = \frac{- R_{m n} (a_{i j} + 2 δ a_{i j}) + 8 R_{m n} (a_{i j} + δ a_{i j}) - 8 R_{m n} (a_{i j} - δ a_{i j}) + R_{m n} (a_{i j} - 2 δ a_{i j})}{12 δ a_{i j}} \dots + O (δ^{4})$

where

δ a_{i j}

is the finite perturbation of the ij-component of the 2nd order orientation tensor. A perturbation size of

δ = 10^{- 6}

is applied individually to each component of

a_{i j}

to define

δ a_{i j}

. We also consider the sensitivity of the result to the finite difference type including the nth order backward

({B W}_{n})

, forward

({F W}_{n})

and central

({C T}_{n})

differences, (see, e.g., [44], for finite difference formulas). The model parameters used here for the various 2nd order fiber orientation tensor models considered for this validation exercise appear in Table 2 [23].

We assume for this validation exercise a ‘randomly’ generated orientation state ${\underline{\underline{a}}}^{0}$ , given as

${\underline{\underline{a}}}^{0} = [\begin{matrix} 0.0622 & 0.0765 & 0.0398 \\ 0.0765 & 0.5521 & 0.0186 \\ 0.0398 & 0.0186 & 0.3857 \end{matrix}]$

Values of

{e r r}_{2}

from Equation (125) appear in Table 3, Table 4 and Table 5 for the various 2nd order orientation tensor models, identified in the left most column, and 4th order closure approximations, as identified in the table headers. The results show that with a perturbation size of

δ = 10^{- 6}

, the

{e r r}_{2}

magnitude values are below

10^{- 6}

which ensures the accuracy of our exact Jacobian derivations and numerical implementations over the various 2nd order orientation tensor models and 4th order closures approximations considered above. The

{e r r}_{2}

magnitude values appear lower for the Hinch and Leal closure approximations (cf. Table 3) on the order of

10^{- 9} - 10^{- 8}

as compared to the EBF and IBF closures (cf. Table 4), which is on the order of

10^{- 8}

. Alternatively, the

{e r r}_{2}

magnitude values for the IBOF and EBOF closures (cf. Table 5) are relatively higher on the order of about

10^{- 7}

. In Figure 2a, the

{e r r}_{2}

is shown to decrease monotonically with increasing perturbation size, showing the optimum result for

δ = 10^{- 3}

, below which

{e r r}_{2}

increases monotonically. Additionally, Figure 2b shows that the

{e r r}_{2}

depends on the finite difference type and decreases with increasing order of approximation. Values of

{e r r}_{2}

are relatively higher with lower order finite difference approximations (i.e., the 1st order backward and forward difference

{B W}_{1} a n d {F W}_{1}

, respectively) and lower with higher order approximations (i.e., the 4th order central difference

{C T}_{4}

3.2. Comparison of the NR Method and the Dormand-Prince Runge–Kutta (RK45) Method

In this section, computed results of the steady state orientation states obtained for various orientation evolution models and closure approximations using our Newton Raphson algorithm are compared to those obtained using the RK45 method. The RK45 method is implemented to provide transient solutions to the initial-value orientation ODE (cf. Equation (7)). A sufficiently small time step is selected for the RK45 simulations to minimize truncation errors [44] and the computations are performed until steady state is achieved, defined where the Euclidean norm of the orientation tensor material derivative is less than $10^{- 6}$ , i.e., at ${e r r}_{R K 45} (t) = {‖{\underline{R}}_{i j}^{R K 45} (t)‖}_{2}, {\underline{R}}_{i j}^{R K 45} (t) = \partial a_{i j}^{R K 45} / \partial t$ .

Three sample cases are included here; the first set of models is based on the study by Falvoro et al. [23], and two other sets are based on the study by Tseng et al. [10]. These cases are used to assess the performance of the NR method relative to the RK4 method for steady state predictions based on various 2nd order orientation tensor models. The first case includes various anisotropic rotary diffusion (ARD) tensor models presented in [23] that account for fiber-fiber and fiber-matrix interactions. The second case considers slow orientation kinetics (SOK) tensor-based models such as the SRF, RSC, and RPR that correct the over-prediction inaccuracies of the transient fiber orientation response computed with the AT model. These models account for the effect of both deformation and diffusion of the fibers and matrix on the transient fibers response [10]. The third case is based on coupled tensor models involving a combination of an ARD model and an SOK model that has been shown to yield improved orientation tensor accuracy including the iARD-RPR and pARD-RPR which have been incorporated into Moldex3D (CoreTech System Co. of Taiwan) [19] and iARD-RSC models incorporated into AutoDesk MoldFlow^® (Autodesk Inc., San Francisco, CA, USA). The models in the latter case study are typically used in injection molding simulations specifically for long-fiber concentrated suspensions [10,22]. The EBOF closure approximation of Verweyst [32] is employed for all analyses.

Case 1: The first study considers various anisotropic rotary diffusion (ARD) orientation tensor models with model parameters appearing in Table 2. The initial orientation state for the RK45 simulations is

(124) ${\underline{\underline{a}}}^{0} = [\begin{matrix} 0.30 & 0.00 & 0.00 \\ 0.00 & 0.60 & 0.10 \\ 0.00 & 0.10 & 0.10 \end{matrix}]$

which is the same as that used as the initial guess in the NR method calculations. The initial guess value is selected such that the target steady state orientation is physical with non-negative eigenvalues in all cases. This is required for some tensor models that can exhibit numerical convergence instability such as the PT model [10,23]. One purpose of this first case study is to show that the NR method adequately predicts the steady state orientation irrespective of the choice and complexity of the rotary diffusion model or transient simulation behavior. The RK45 calculated transient responses for select components of

a_{i j},

and the results appear in Figure 3.

The computed responses in Figure 3a,b show varying amplitudes of initial overshoot of $a_{11}$ and initial trough of $a_{22}$ depending on the orientation model that do not affect the final steady state orientation. Table 6 includes values of ${e r r}_{1}$ (cf. Equation (121)) which show that the NR-based steady state orientation predictions agree well with RK4 values indicated by ${e r r}_{1} < 0.005 %$ for both diagonal and off-diagonal tensor components. Although the NR method does not provide information about the transient orientation history such as the initial overshoot and/or the initial trough, steady state orientations appear to be independent of the transient behavior of the orientation tensor model.

The NR convergence history for the various diffusion tensor models appears in Figure 4, which indicates a relatively higher rate of convergence with a quadratic nature for the fundamental AT isotropic rotary diffusion (IRD) model, requiring only 5 iterations for convergence compared to the PT, WPT, Dz, ARD models, which all required 7 iterations to converge. The improved and principal ARD models (iARD and pARD) show the slowest convergence rate, requiring 8 iterations to converge. All models, however, exhibit stable convergence with a monotonically decreasing NR residual norm. With the RK45 method, the orientation tensor reached steady state faster for the AT model, at a time fraction of $\dot{γ} t = 35 s$ compared to the iARD, pARD, Dz, and MRD models, which required an average of about $\dot{γ} t = 126 s$ to converge. The PT and WPT models were the slowest, requiring at least $\dot{γ} t = 185 s$ for the orientation to approach steady state. Unlike the AT, PT, and WPT models, the iARD, pARD, Dz, and MRD model transient responses had undulations, which may be an effect of the numerical solver due to the stiffness of the ODE.

In terms of computational efficiency, computation time was obtained for all NR and RK5 simulations using a x64-based, Intel(R) Core(TM) i5-8365U CPU with a 1.90 GHz processing speed and a 16.0 GB RAM, operating on Windows 11 Pro. The average computation time for each technique appears in Table 7. Overall, the NR technique is seen to be, on average, 7.4 x faster than the RK45 method for the error tolerance given in Table 6.

Case 2: Unlike the first case study, Case 2 investigates the behavior of the NR method for SOK tensor models that involves corrections to slow down the transient fiber orientation kinematic response as compared to the AT model, which includes the SRF, RSC and RPR models. Although the various model corrections have no effect on the final steady state solution, the effect of the SOK model complexity on NR convergence is shown. The model parameters used for this case study are taken from Tseng et al. [10] and appear in Table 8.

The starting orientation for the RK45 analysis and the initial guess ${\underline{\underline{a}}}^{0}$ used for the Newton Raphson method is the same as that used in case study 1 (cf. Equation (124)). Two flow fields L₁ and L₂ are considered:

L₁: Simple shear flow in the 1–2 plane, $L_{12} = \dot{γ}$ .

L₂: Balanced shear/planar-elongation flow, simple shear in 1-2 plane superimposed on planar elongation in 1–2 plane. $L_{11} = - \dot{ε}, L_{22} = \dot{ε}, L_{12} = \dot{γ}$ given that $\dot{γ} / \dot{ε} = 10$ .

The time evolution of the components of the 2nd order orientation tensor based on the RK45 method are shown in Figure 5 below.

The percentage error between NR and RK45 steady state values are presented in Table 9. The results show a high level of accuracy in prediction based on the NR method, with a maximum error tolerance of only ${e r r}_{1} < 0.04 %$ among all tensor components for all models. The AT and SRF models exhibited relatively faster convergence rate requiring only 5 iterations to converge for both flow types (cf. Figure 6) compared to the RSC and RPR models which required an extra iteration to converge. The orientation tensor reached steady state faster with the AT tensor model compared to other tensor models in Figure 5 because the RK45 residual approached the set tolerance value of $10^{- 6}$ faster than other tensor models (cf. Figure 6b). The RSC and RPR models both require complex eigenvalue and eigenvector spectral decompositions which may influence the convergence rate. Moreover, the derivatives of the eigenvectors as shown in Appendix A may impact the convergence rate of these models. The convergence rate was seen to be independent of both flow types considered.

Figure 7a,b present the NR residual norm for the AT and RSC models and for different initial orientation states obtained by varying $a_{22}^{0}$ (i.e., $0 . \leq a_{22}^{0} \leq 0.7$ ) while retaining the initial states of other $a_{i j}^{0}$ components. As $a_{22}^{0}$ increases, the initial orientation states approach the final steady state orientation; thus, less iterations are required for the solutions to converge. The AT converged as expected for all $a_{22}^{0}$ values considered (cf. Figure 7a); however, the RSC model was very sensitive to the $a_{22}^{0}$ values (cf. Figure 7b). Only the initial orientation states with $a_{22}^{0} = 0.49, 0.56, 0.63, 0.70$ converged to physical solutions (as shown in Figure 7b), while the other $a_{22}^{0}$ values yielded non-physical solutions.

For the purpose of visualization, we present a typical 3D orientation distribution function (ODF) (cf. Figure 8) recovered from the deviatoric forms of the average steady state 2nd and 4th order fiber orientation tensors obtained from the NR solutions for two fiber interaction coefficients (C_I = 0.01 and C_I = 1.0) using the AT model and based on the method of Onat and Leckie [7].

Case 3: The third case study investigates the performance of the NR method to coupled iARD-RPR and pARD-RPR orientation tensor models and includes model parameters to simulate three different polymer composites materials in simple shear flow with $\dot{γ} = 1 s^{- 1}$ . The model parameters used for the three coupled models and three polymer composites (i.e., M1: 40 wt. % glass-fiber/PP, M2: 31 wt. % carbon-fiber/PP and M3: 40 wt. % glass-fiber/nylon) considered in this case study appear in Table 10 and Table 11. We assume the same RK45 initial orientation state and NR initial guess as in case study 1 (cf. Equation (124)). The simulated RK45 results for the coupled orientation models and polymer composites considered here appear in Figure 9a–c. The percentage error between the NR and RK45 steady state values are given in Table 12. The results show negligible discrepancy in the calculated steady state solutions, generally yielding an error ${e r r}_{1}$ below $0.3 %$ for all tensor components. Figure 10a, however, shows that the NR convergence rate is dependent on the model and composite material composition. The M2 (31 wt. % carbon-fiber/PP) with lower fiber content had a faster convergence rate for all three coupled tensor models compared to the M1 and M3 composites with higher fiber content which had a slower convergence rate. The coupled models showed slight instability through the initial NR iterations, especially for the iARD-RPR model with M1. The RK45 convergence profiles (cf. Figure 10b) is seen to have similar behavior as the NR convergence profiles in terms of the order of convergence for the various models and material types.

3.3. Performance of NR Method Using Various Closure Approximations

The performance of the NR method in computing steady state orientation tensor component values for various closure approximations of the 4th order orientation tensor in terms of accuracy and stability is assessed in this section. For this assessment, we consider the AT model with $C_{I} = 0.01$ . The same initial orientation state for the RK45 method and the NR method’s initial guess used in case 1 of the preceding section (cf. Equation (124)) is also employed, and simple shear flow with $\dot{γ} = 1 s^{- 1}$ is used for all analyses. The first set of results are presented for Hinch and Leal (HL) closures in Section 2.5.4 and for fitted closures in Section 2.5.5.

Figure 11a–d show the transient response of the $a_{11} a n d a_{12}$ from RK45 analysis through steady state. These results show that the choice of closure approximation affects the time taken to reach a steady state orientation. The model with the QDR and HYB₂ closures reached steady-state faster than other models, while models with the LIN, ISO and SF2 closures required more time to reach steady-state. As shown in Table 13, the NR method results tend to agree well with RK45 reference values, having a maximum ${e r r}_{1} = 0.004 %$ . Likewise, NR simulations using the fitted closure approximations performed well with a good accuracy in predicted results (cf. Table 13). Figure 12a shows that the LIN and ISO closures had the fastest convergence rate with a quadratic nature, requiring only 3 iterations to converge followed by the HL1, HL2, QDR and HYB₁ closures, which required 6 iterations to converge. The invariant-based HYB₂ closure required an extra iteration to converge, while the SF2 closure had the slowest convergence rate, requiring up to 10 iterations to converge. Surprisingly, all the fitted closure approximations had similar convergence behavior, with only 5 iterations needed to converge except the IBOF closure, which required an extra iteration to converge (cf. Figure 12b). Conversely, RK45 convergence using the Hinch and Leal closures (cf. Figure 12c) has a rather distinct behavior, with the LIN, ISO and SF2 closures having the slowest convergence rate and the HYB₂ and QDR closures having the fastest convergence rate. Similarly to the NR results that were obtained using fitted closures, the RK45 simulations converged over a narrower range (cf. Figure 12d) as compared to the Hinch and Leal closures.

3.4. Homogenous Flow Considerations

In this section we employ six homogenous flows to ensure the stability of our NR method beyond simple shear flow. The following flow conditions taken from [48] were considered:

SS: Simple shear flow in the 2–3 plane, i.e., $L_{23} = \dot{γ}$ .

SUA: Two stretching/shearing flows, including simple shear in 2–3 plane superimposed with uniaxial elongation in 3-direction, i.e., $L_{11} = - \dot{ε}, L_{22} = - \dot{ε}, L_{33} = 2 \dot{ε}, L_{23} = \dot{γ}$ . Two cases are considered, including balanced shear/stretch with $\dot{γ} / \dot{ε} = 10$ and dominant stretch with $\dot{γ} / \dot{ε} = 1$ .

UA: Uniaxial elongation flow in the 3—direction, i.e., $L_{11} = L_{22} = - \dot{ε}, L_{33} = 2 \dot{ε}$ .

BA: Biaxial elongation (BA) flow in the 2–3 plane, i.e., $L_{11} = - 2 \dot{ε}, L_{22} = L_{33} = \dot{ε} .$

PST: Two shear/planar-elongation flows, including simple shear in 2–3 plane superimposed on planar elongation in 1–3 plane, $L_{11} = - \dot{ε}, L_{33} = \dot{ε}, L_{23} = \dot{γ}$ . Two cases are considered, including balanced shear-planar elongation with $\dot{γ} / \dot{ε} = 10$ and dominant planar elongation with $\dot{γ} / \dot{ε} = 1$ .

SBA: Balanced shear/bi-axial elongation flow, including simple shear in 2–3 plane superimposed on biaxial elongation, i.e., $L_{11} = - 2 \dot{ε}, L_{22} = L_{33} = \dot{ε}, L_{23} = \dot{γ}$ . Two cases are considered for the shear to extension rate, $\dot{γ} / \dot{ε} = 1$ and $\dot{γ} / \dot{ε} = 10$ .

The natural orthotropic closures, i.e., exact midpoint fit, NAT₁, was employed for this study. The initial orientation state for the RK4 simulations and the initial guess ${\underline{\underline{a}}}^{0}$ for the NR iterations is

(125) ${\underline{\underline{a}}}^{0} = [\begin{matrix} 0.40 & 0.00 & 0.00 \\ 0.00 & 0.00 & 0.10 \\ 0.00 & 0.10 & 0.45 \end{matrix}]$

The NR method provided exceptional prediction accuracy of the orientation tensor components for all flow conditions considered where the percentage error for all orientation tensor components was zero to four decimals (i.e., ${e r r}_{1} = 0.0000 %$ ). Figure 13a,b show transient results for $a_{33} a n d a_{23}$ from RK45 simulations, respectively. The NR iteration history for the same (cf. Figure 13c) shows that a higher shear-to-extension rate ratio improves the NR rate of convergence. The SS, PST₂, SUA_2, and SBA₂ flows converged faster than other flow types. Conversely, the flow types with the highest extension-to-shear rate ratio (UA and BA) had a slower rate of convergence. Alternatively, the RK45 solution (cf. Figure 13d) showed that the higher extension-to-shear rate ratio enhanced the convergence rate, contrary to what was observed for the NR convergence history. The behavior of both NR and RK45 numerical techniques under the various flow and model considerations can provide useful insights when selecting a suitable numerical technique for solving a particular problem.

4. Conclusions

A Newton-Raphson (NR) method was derived and implemented for computing steady state 2nd order fiber orientation tensors in homogeneous flows using exact 4th order Jacobians. Analytical expressions are provided for the Jacobians which are derived from partial derivatives of 2nd order fiber orientation tensor material derivatives taken with respect to the 2nd order fiber orientation tensor itself. Various macroscopic fiber orientation moment-tensor models and closure approximations of the 4th order fiber orientation tensor are included in the derivations, and the performance of the NR method in various homogenous flows have been studied.

In summary, our results show that

The computed exact Jacobians aligned closely with finite difference approximations for all fiber orientation tensor models and 4th order fiber orientation tensor closure approximations considered, and the degree of alignment depended on the perturbation size, finite difference type and finite difference order of approximation.

The steady state orientation predictions from the NR solutions matched closely with the results predicted by the RK45 method for all fiber orientation tensor models, 4th order fiber orientation tensor closure approximations and homogenous flow types considered in this paper. The NR method was observed to be comparatively faster than the RK45 method for computing steady state values of the second order orientation tensor in all cases.

The NR convergence and stability behavior depended on the type and complexity of the fiber orientation tensor model, the initial guess of the orientation tensor and the material properties of the underlying composite material. A good initial guess resulted in faster and more stable convergence to physical steady state solutions.

The NR convergence behavior also depended on the flow type and extension-to-shear rate ratio. Higher shear rate dominance and flow symmetry increased the NR convergence rate and vice versa.

Overall, the NR method for determining the steady state orientation tensor performed well in terms of accuracy and computational efficiency for all fiber orientation tensor models, 4th order fiber orientation tensor closure approximations and homogenous flows considered in this article.

Author Contributions

Conceptualization, A.E.A. and D.E.S.; methodology, A.E.A. and D.E.S.; software, D.E.S.; validation, A.E.A. and D.E.S.; formal analysis, A.E.A.; investigation, A.E.A.; resources, D.E.S.; data curation, A.E.A.; writing—original draft preparation, A.E.A.; writing—review and editing, D.E.S.; visualization, A.E.A.; supervision, D.E.S.; project administration, D.E.S.; funding acquisition, D.E.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data presented in this study are available upon request from the corresponding author due to privacy restrictions.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Single ‘rigid’ ellipsoidal fiber orientation.

Figure 2 Finite difference error ${e r r}_{2}$ for various orientation models and closures (a) over perturbation size $δ$ for the 4th order central difference approximation and for (b) finite difference type and order of approximation at their corresponding optimum perturbation size $δ$ .

Figure 3 Time evolution of the 2nd order orientation tensor for calibrated FT, PT, iARD, pARD, WPT, Dz and MRD models for (a) $a_{11}$ component and (b) $a_{22}$ component.

Figure 4 Values of the residual norm for the (a) NR method and (b) RK45 method for the AT, PT, iARD, pARD, WPT, Dz and MRD models.

Figure 5 Time evolution of the $a_{11}, a_{22} a n d a_{13}$ components of the 2nd order orientation tensor for calibrated AT, SRF, RSC and FT-RPR models for (a) simple shear flow and (b) shearing/stretching combination flow.

Figure 6 History of the (a) NR residual norm and (b) RK45 absolute maximum material derivative tensor component for the FT, RSC, SRF, and RPR models for the 2 different flow fields (L1 and L2).

Figure 7 NR residual norm for different initial orientation states $a_{i j}^{0}$ and for (a) AT and (b) RSC orientation tensor models.

Figure 8 Reconstruction of the orientation distribution function (ODF) for the steady-state orientation for fiber interaction coefficients of (a) C_I = 1.0 and (b) C_I = 0.01.

Figure 9 Time evolution of $a_{11}, a_{22} a n d a_{13}$ for the iARD-RPR, iARD-RSC, and pARD-RPR models using materials (a) M1, (b) M2, and (c) M3.

Figure 10 Iteration history of the (a) NR residual norm and (b) RK45 residual norm for the various coupled models and polymer composites.

Figure 11 Transient simulation results for (a) $a_{11}$ and (b) $a_{12}$ for the various Hinch and Leal closure approximations and (c) $a_{11}$ and (d) $a_{12}$ for the fitted closures.

Figure 12 Iteration history of the NR residual norm for the (a) Hinch and Leal closures and (b) polynomial fitted closure approximations, and iteration history of the RK45 residual norm for the (c) Hinch and Leal closures and (d) polynomial fitted closure approximations.

Figure 13 Transient responses of the 2nd order orientation tensor evolution for (a) component $a_{33}$ and (b) component $a_{23}$ , and residual norm history profiles for the (c) NR method and (d) RK45 method. Results are presented for various flow considerations.

Table 2

Case Study 1 parameters for the FT, Dz, iARD, pARD, WPT, MRD and PT models [23].

	$C_{I}$	ARD Parameters
$F T$	$0.0311$	-
$D z$	$0.0258$	$D_{z} = 0.051$ , $\hat{n} = [\begin{matrix} 0 & 0 & 1 \end{matrix}]$
$i A R D$	$0.0562$	$C_{M} = 09977$
$p A R D$	$0.0169$	$Ω = 0.9868$
$W P T$	$0.0504$	$w = 0.9950$
$M R D$	$0.0198$	$D_{1} = 1.0, D_{2} = 0.7946, D_{3} = 0.0120$
$P T$	-	$b_{1}$	$b_{2}$	$b_{3}$	$b_{4}$	$b_{5}$
$P T$	-	$1.924 \times 10^{- 4}$	$5.839 \times 10^{- 3}$	$0.04$	$1.168 \times 10^{- 5}$	$0$

Table 3

Computed error $(\times 10^{- 8})$ for various orientation tensor models using Hinch and Leal closure approximations.

	HYB₁	HYB₂	ISO	LIN	QDR	SF2	HL₁	HL₂
AT	0.6436	0.9385	0.2220	0.4188	0.2691	2.0949	0.9618	4.3940
PT	0.8088	0.7549	0.5837	0.5003	0.4244	1.6776	0.8241	3.4809
iARD	0.5737	1.2712	0.3444	0.6336	0.5728	0.5100	0.8774	1.6148
pARD	0.7169	0.5475	0.2722	0.2438	0.4155	1.4805	0.9818	3.6185
WPT	0.8563	1.0386	0.3773	0.2525	0.2926	1.4543	0.9632	3.5284
Dz	0.5899	0.8248	0.2373	0.5484	0.3233	0.7137	1.0594	2.7732
NEM	0.6490	0.9306	0.4012	0.4314	0.1612	2.1012	0.9846	4.4062
pARD-RSC	1.0030	1.3343	1.3062	1.0506	1.3699	0.5378	1.3441	1.6482
iARD-RPR	0.5645	0.6900	0.3478	0.3687	0.5222	1.0731	1.0324	2.0512

Table 4

Computed error $(\times 10^{- 7})$ for various orientation tensor models using IBF and EBF closure approximations.

	NAT₁	NAT₂	ORS	ORT	ORW	ORW3
AT	0.3557	0.5812	0.4573	0.3351	0.5994	0.6496
PT	0.3076	0.3946	0.3327	0.2517	0.5693	0.5284
iARD	0.2512	0.2663	0.2129	0.1934	0.4762	0.4213
pARD	0.2975	0.3956	0.3216	0.2479	0.5354	0.5229
WPT	0.3040	0.4068	0.3430	0.2622	0.5720	0.5612
Dz	0.3532	0.3276	0.2805	0.2920	0.6704	0.6010
NEM	0.3649	0.5869	0.4615	0.3388	0.6062	0.6520
pARD-RSC	0.2529	0.2800	0.2074	0.1802	0.4772	0.4155
iARD-RPR	0.1964	0.1709	0.1601	0.1404	0.3687	0.2820

Table 5

Computed error $(\times 10^{- 7})$ for various orientation tensor models using IBOF and EBOF closure approximations.

	IBOF	WTZ	LAR32	VST	FFLAR4	LAR4
AT	6.1748	4.3147	5.0800	3.1567	4.3188	4.3101
PT	5.6437	4.0286	4.7162	2.9443	4.0435	3.9967
iARD	4.4942	3.0776	3.6782	2.2662	3.0115	3.0665
pARD	5.5458	3.8741	4.5415	2.8548	3.8500	3.8818
WPT	5.6412	4.0391	4.7234	2.9350	4.0486	3.9851
Dz	6.4226	4.8010	5.5454	3.3924	4.8016	4.6691
NEM	6.1978	4.3254	5.0936	3.1653	4.3271	4.3182
pARD-RSC	4.1821	2.9401	3.4878	2.1358	2.9236	2.9070
iARD-RPR	3.4882	2.4509	2.9082	1.7359	2.3590	2.4110

Table 6

Error estimates of the $a_{11}, a_{22} a n d a_{13}$ steady-state orientation tensor component values for FT, Dz, iARD, pARD, WPT, MRD and PT models.

	$a_{11}$	$a_{22}$	$a_{13}$
FT	0.0015	0.0003	0.0000
PT	0.0046	0.0021	0.0031
iARD	0.0009	0.0006	0.0000
pARD	0.0006	0.0003	0.0032
WPT	0.0043	0.0020	0.0030
Dz	0.0021	0.0011	0.0021
MRD	0.0015	0.0008	0.0046

Table 7

Average processor time taken (sec) to compute the steady state orientation values for both numerical techniques and for the various diffusion models.

	FT	PT	iARD	pARD	WPT	Dz	MRD
NR	0.7210	0.4075	0.5916	0.3784	0.3038	0.4615	0.7908
RK45	3.5829	3.9775	3.0067	2.6441	2.7791	3.4486	6.7519

Table 8

Case Study 2 parameters for the FT, SRF, RSC and RPR models [10].

	$A T$	$S R F$	$R S C$	$R P R$
$C_{I}$	$0.01$	$0.01$	$0.01$	$0.01$
$ᶄ$	$-$	$0.1$	$0.1$	$-$
$α$	$-$	$-$	$-$	$0.9$
$β$	$-$	$-$	$-$	$0$

Table 9

Calculated $a_{11}, a_{22} a n d a_{13}$ steady-state error values for RSC, AT, SRF, and RPR models and for the 2 flow fields L₁ and L₂.

	L₁			L₂
	$a_{11}$	$a_{22}$	$a_{13}$	$a_{11}$	$a_{22}$	$a_{13}$
RSC	0.0374	0.0058	0.0000	0.0178	0.0122	0.0017
FT	0.0036	0.0005	0.0000	0.0015	0.0011	0.0000
SRF	0.0374	0.0055	0.0024	0.0158	0.0109	0.0050
RPR	0.0374	0.0058	0.0000	0.0178	0.0122	0.0017

Table 10

ARD-RSC parameters [10] for (a) M1—40 wt. % glass-fiber/PP; (b) M2—31 wt. % carbon-fiber/PP; and (c) M3—40 wt. % glass-fiber/nylon.

	M1	M2	M3
$ᶄ$	1/30	1/30	1/20
$b_{1}$	$3.842 \times 10^{- 4}$	$3.728 \times 10^{- 3}$	$4.643 \times 10^{- 4}$
$b_{2}$	$- 1.786 \times 10^{- 3}$	$- 1.695 \times 10^{- 2}$	$- 6.169 \times 10^{- 4}$
$b_{3}$	$5.250 \times 10^{- 2}$	$1.750 \times 10^{- 1}$	$1.900 \times 10^{- 2}$
$b_{4}$	$1.168 \times 10^{- 5}$	$- 3.367 \times 10^{- 3}$	$9.650 \times 10^{- 4}$
$b_{5}$	$- 5.000 \times 10^{- 4}$	$- 1.000 \times 10^{- 2}$	$7.000 \times 10^{- 4}$

Table 11

Coupled iARD-RPR and pARD-RPR model parameters [10] for (a) model M1: 40 wt. % glass-fiber/PP, (b) model M2: 31 wt. % carbon-fiber/PP, and (c) model M3: 40 wt. % glass-fiber/nylon.

	M1	M2	M3
$C_{I}$	0.0165	0.0630	0.0060
$C_{M}$	0.9990	1.0100	0.9000
$Ω$	0.9880	0.9650	0.9000
$α$	0.9650	0.9650	0.9500
$β$	0.0000	0.0000	0.0000

Table 12

Calculated errors for $a_{11}, a_{22} a n d a_{13}$ steady-state values for iARD-RPR, iARD-RSC, and pARD-RPR models and materials M1, M2 and M3.

	iARD-RPR			pARD-RPR			iARD-RSC
	M1	M2	M3	M1	M2	M3	M1	M2	M3
$a_{11}$	0.1039	0.0757	0.0876	0.0878	0.0278	0.0971	0.1107	0.0395	0.3231
$a_{22}$	0.0550	0.0255	0.0219	0.0444	0.0131	0.0245	0.0553	0.0202	0.0678
$a_{13}$	0.1286	0.1703	0.0401	0.0656	0.0270	0.0415	0.0872	0.0450	0.0235

Table 13

Calculated error for $a_{11}, a_{22} a n d a_{12}$ steady-state values for (a) Hinch and Leal closure approximations and (b) fitted closure approximations.

(a)				(b)
	$a_{11}$	$a_{22}$	$a_{12}$		$a_{11}$	$a_{22}$	$a_{12}$
HYB₁	0.0042	0.0016	0.0000	IBOF	0.0068	0.0000	0.0000
HYB₂	0.0000	0.0000	0.0000	ORS	0.0000	0.0000	0.0000
ISO	0.0000	0.0000	0.0000	ORT	0.0061	0.0013	0.0000
LIN	0.0000	0.0000	0.0000	NAT₁	0.0000	0.0000	0.0000
QDR	0.0000	0.0000	0.0000	ORW	0.0067	0.0000	0.0000
SF2	0.0055	0.0000	0.0000	NAT₂	0.0072	0.0000	0.0000
HL1	0.0042	0.0014	0.0000	WTZ	0.0060	0.0000	0.0000
HL2	0.0038	0.0000	0.0000	LAR32	0.0066	0.0000	0.0000
				ORW3	0.0066	0.0000	0.0000
				VST	0.0066	0.0000	0.0000
				FFLAR4	0.0062	0.0013	0.0000
				LAR4	0.0066	0.0000	0.0000

Appendix A

Appendix A.1. Eigenvalues and Eigenvector Derivatives

The eigenvalue for a matrix $\underline{\underline{a}}$ is often given in terms of the eigenvalues $λ^{k}$ and corresponding eigenvectors ${\underline{Φ}}^{k}$ as [49](A1) $[\underline{\underline{a}} - λ_{k} \underline{\underline{I}}] {\underline{Φ}}^{k} = {\underline{F}}^{k}$ For the homogeneous system ${\underline{F}}^{k} = 0$ , requiring ${\underline{Φ}}^{k} \neq 0$ to yield non-trivial solutions, Equation (A1) becomes $|\underline{\underline{a}} - λ_{k} \underline{\underline{I}}| = 0$ , where the vertical bars indicate the matrix determinant. The resulting system matrix $[a_{m n} - λ^{k} δ_{m n}]$ is known to be rank deficient by an order of 1, requiring that we impose a normalization condition on the eigenvectors ${\underline{Φ}}^{k}$ which can be defined through a scalar function $G^{k} ({\underline{Φ}}^{k}) = 0$ . The eigenvectors are thus obtained by replacing the nth row of the residual column vector $Σ_{i}^{k} = (a_{i j} - λ^{k} δ_{i j}) Φ_{j}^{k} - {\underline{F}}_{i}^{k}$ with $G^{k} ({\underline{Φ}}^{k}) = 0$ and solving for ${\underline{Φ}}^{k}$ from the equation ${\underline{Σ}}^{k} ({\underline{Φ}}^{k}) = 0$ through any iterative algorithm or explicit solver. Here, we employ Newton-Raphson’s method to obtain ${\underline{Φ}}^{k}$ such that(A2) ${\underline{Φ}}^{k}^{+} = {\underline{Φ}}^{k}^{-} - {\underline{\underline{J}}}^{k}^{- 1} Σ^{k}$ (A3) $Σ_{i}^{k} = \{\begin{matrix} S_{i j}^{k} Φ_{j}^{k} - {\underline{F}}_{i}^{k} & i \neq n \\ G^{k} ({\underline{Φ}}^{k}) & i = n \end{matrix}$ (A4) $J_{i j}^{k} = \frac{\partial Σ_{i}^{k}}{\partial Φ_{j}^{k}} = \{\begin{matrix} S_{i j}^{k} & i \neq n \\ \frac{\partial G^{k}}{\partial Φ_{j}^{k}} & i = n \end{matrix}$ where(A5) $S_{i j}^{k} = a_{i j} - λ^{k} δ_{i j}$ The derivative of the eigenvalues with respect to components $a_{i j}$ can thus be obtained by differentiating Equation (A3), assuming symmetry of system matrix, i.e., $S_{i j}^{k} = S_{j i}^{k}$ , such that(A6) $Φ_{i}^{k} \frac{\partial S_{i j}^{k}}{\partial a_{r s}} Φ_{j}^{k} + {\underline{F}}_{i}^{k} \frac{\partial Φ_{i}^{k}}{\partial a_{r s}} - Φ_{i}^{k} \frac{\partial F_{i}^{k}}{\partial a_{r s}} = 0$ (A7) $Φ_{i}^{k} S_{i j}^{k} = S_{j i}^{k} Φ_{j}^{k} = S_{i j}^{k} Φ_{j}^{k} = {\underline{F}}_{i}^{k}$ where(A8) $\frac{\partial S_{i j}^{k}}{\partial a_{r s}} = δ_{i r} δ_{j s} - \frac{\partial λ^{k}}{\partial a_{r s}} δ_{i j}$ Assuming ${\underline{F}}_{i}^{k} = 0$ , we obtain(A9) $\frac{\partial λ^{k}}{\partial a_{r s}} = \frac{1}{(Φ_{m}^{k} Φ_{m}^{k})} \{Φ_{r}^{k} Φ_{s}^{k}\}$ Consequently, we can determine the corresponding derivatives $\partial Φ_{i}^{k} / \partial a_{r s}$ given $\partial λ^{k} / \partial a_{r s}$ from(A10) $S_{i j}^{k} \frac{\partial Φ_{i}^{k}}{\partial a_{r s}} = - \frac{\partial S_{j i}^{k}}{\partial a_{r s}} Φ_{i}^{k}$ Recalling that the system matrix $[\underline{\underline{a}} - λ^{k} \underline{\underline{I}}]$ is inherently singular, we can substitute an nth row of the above equation by adopting a normalization technique [49]: (A11) $\frac{\partial G^{k}}{\partial Φ_{j}^{k}} \frac{\partial Φ_{j}^{k}}{\partial a_{r s}} = \frac{\partial G^{k}}{\partial a_{r s}}$ The modified differential equation can be recast as given in Equation (A12) below, allowing the inversion of the modified system matrix to solve(A12) $J_{i j}^{k} \frac{\partial Φ_{i}^{k}}{\partial a_{r s}} = \frac{\partial Q_{i}^{k}}{\partial a_{r s}}$ where(A13) $\frac{\partial Q_{i}^{k}}{\partial a_{r s}} = \{\begin{matrix} - \frac{\partial S_{i j}^{k}}{\partial a_{r s}} Φ_{j}^{k} & i \neq n \\ - \frac{\partial G^{k}}{\partial a_{r s}} & i = n \end{matrix}$ Smith et al. [49] presents three common normalization techniques:

‘Mass’ Normalization

(A14) $G^{k} ({\underline{Φ}}^{k}) = Φ_{j}^{k} Φ_{j}^{k} - 1, \frac{\partial G^{k}}{\partial Φ_{j}^{k}} = 2 Φ_{j}^{k}, \frac{\partial G^{k}}{\partial a_{r s}} = 0$

Preassigning an mth Component of ${\underline{Φ}}^{k}$

(A15) $G^{k} ({\underline{Φ}}^{k}) = δ_{m j} Φ_{j}^{k} - α, \frac{\partial G^{k}}{\partial Φ_{j}^{k}} = δ_{m j}, \frac{\partial G^{k}}{\partial a_{r s}} = 0$

Predefining the Euclidean norm of ${\underline{Φ}}^{k}$

(A16) $G^{k} ({\underline{Φ}}^{k}) = \sqrt{Φ_{j}^{k} Φ_{j}^{k}} - β, \frac{\partial G^{k}}{\partial Φ_{j}^{k}} = Φ_{j}^{k}, \frac{\partial G^{k}}{\partial a_{r s}} = 0$

Appendix B

Appendix B.1. Optimal Fitted Closure Approximation Constants/Coefficients

Appendix B.1.1. Eigenvalue-Based Optimal Fitting Closure (EBOF) Approximation

In the main body of the paper, we presented the generalized formulations for the eigenvalue-based orthotropic binomial fitted closure approximation which is written in terms of a constant coefficient matrix $C_{i j}^{(n)}$ and an nth order permuted bivariate polynomial vector $⩓_{j}^{(n)}$ with $(n + 1) (n + 2) / 2$ number of parameters (cf. Equations (101)–(103)). Depending on the polynomial order $(n)$ the various constants and coefficient matrix ${\underline{\underline{C}}}^{(n)}$ for the closures considered in this paper taken from Kuzmin [38] are presented below.

Linear Orthotropic Fitted Closure $(n = 1)$

For the orthotropic linear (ORT) closure, the constant coefficient matrix ${\underline{\underline{C}}}^{(1)}$ is given as ${\underline{\underline{C}}}^{(1)} = \frac{1}{7} [\begin{matrix} - 3 / 5 & 6 & 0 \\ - 3 / 5 & 0 & 6 \\ 27 / 5 & - 6 & - 6 \end{matrix}]$ and for the orthotropic fitted closure by Cintra and Tucker (cf. [39]), the constant coefficient matrix ${\underline{\underline{C}}}^{(1)}$ is given as ${\underline{\underline{C}}}^{(1)} = [\begin{matrix} - 0.15 & 1.15 & - 0.10 \\ - 0.15 & 0.15 & 0.90 \\ 0.60 & - 0.60 & - 0.60 \end{matrix}]$

Quadratic Orthotropic Fitted Closures $(n = 2)$

The simple general orthotropic quadratic closure has a constant coefficient matrix ${\underline{\underline{C}}}^{(2)}$ given as ${\underline{\underline{C}}}^{(2)} = [\begin{matrix} 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 1 & - 2 & - 2 & 1 & 2 & 1 \end{matrix}]$ The orthotropic natural closure—exact midpoint fit (NAT₁) [38] has a constant coefficient matrix ${\underline{\underline{C}}}^{(2)}$ given as ${\underline{\underline{C}}}^{(2)} = [\begin{matrix} 0.0708 & 0.3236 & - 0.3776 & 0.6056 & 0.4124 & 0.3068 \\ 0.0708 & - 0.2792 & 0.2252 & 0.2084 & 0.4124 & 0.7040 \\ 1.1880 & - 2.0136 & - 2.1264 & 0.8256 & 1.7640 & 0.9384 \end{matrix}]$ The ORF independent coefficients are derived from a 2nd order polynomial fit of the principal axis data obtained from DFC via a minimization technique. For the orthotropic fitted closure by Cintra and Tucker (ORF) (cf. [39]), the constant coefficient matrix ${\underline{\underline{C}}}^{(2)}$ is given as ${\underline{\underline{C}}}^{(2)} = [\begin{matrix} 0.060964 & 0.371243 & - 0.369160 & 0.555301 & 0.371218 & 0.318266 \\ 0.124711 & - 0.389402 & 0.086169 & 0.258844 & 0.544992 & 0.796080 \\ 1.228982 & - 2.054116 & - 2.260574 & 0.821548 & 1.819756 & 1.053907 \end{matrix}]$ Chung and Kwon [56] improved the ORF and developed the 2nd order ORW closure for wide interaction coefficients that proved to be stable for all flow conditions considered. The improved orthotropic fitted closure (ORW2) by Chung and Kwon (cf. [56]) has a constant coefficient matrix ${\underline{\underline{C}}}^{(2)}$ given as ${\underline{\underline{C}}}^{(2)} = [\begin{matrix} 0.070055 & 0.339376 & - 0.396796 & 0.590331 & 0.411944 & 0.333693 \\ 0.115177 & - 0.368267 & 0.094820 & 0.252880 & 0.535224 & 0.800181 \\ 1.249811 & - 2.148297 & - 2.290157 & 0.898521 & 1.934914 & 1.044147 \end{matrix}]$ Kuzmin [38] presents details on derivations of some orthotropic fitted closures via a numerical bottom-up approach.

Cubic Orthotropic Fitted Closures $(n = 3)$

Recently, higher order polynomial fitted closures were developed for improved accuracy. Although, the orthotropic natural closure—extended quadratic fit (NAT₂) has a cubic order, it is essentially quadratic.

Non-rational Fitted Closure

The constant coefficient matrix ${\underline{\underline{C}}}^{(3)}$ for this closure approximation is given as ${\underline{\underline{C}}}^{(3)} = [\begin{matrix} 0 & 0.5 & 0 & 0.5 & - 0.6 & 0 & 0 & 0.6 & 0.6 & 0 \\ 0 & 0 & 0.5 & 0 & - 0.6 & 0.5 & 0 & 0.6 & 0.6 & 0 \\ 1 & - 1.5 & - 1.5 & 0.5 & 0.4 & 0.5 & 0 & 0.6 & 0.6 & 0 \end{matrix}]$ Chung and Kwon [56] also extended the 2nd order ORW to develop a 3^rd order polynomial ORW3 closure using new flow data from ODFC calculation. For the improved 3^rd order orthotropic fitted closure ORW3 by Chung and Kwon [56], the constant coefficient matrix is given as ${[{\underline{\underline{C}}}^{(3)}]}^{T} = [\begin{array}{r} - 0.1480648093 & - 0.2106349673 & 0.4868019601 \\ 0.8084618453 & 0.9092350296 & 0.5776328438 \\ 0.7765597096 & 1.1104441966 & 0.4605743789 \\ 0.3722003446 & - 1.2840654776 & - 2.2462007509 \\ - 1.7366749542 & - 2.5375632310 & - 4.8900459209 \\ - 1.3431772379 & 0.1260059291 & - 1.9088154281 \\ - 0.0324756095 & 0.5856304774 & 1.1817992322 \\ 0.8895946393 & 1.9988098293 & 4.0544348937 \\ 1.7367571741 & 1.4863151577 & 3.8542602127 \\ 0.6631716575 & - 0.0756740034 & 0.9512305286 \end{array}]$ Eigenvalue-based optimal fitting (EBOF) closures were developed, which are extensions to the ORF but with higher order polynomial fits for improved accuracy. These include the rational elliptical (RE) closure by Wetzel [54] and the 4th order polynomial fitted closure by Verweyst [53] (VST).

Rational Fitted Closure

The rational elliptical (RE) closure developed by Wetzel [41] is a higher order extension to the ORF using Carlson elliptical integrals. The rational ellipsoid fitted closure has two permutation vectors, with the denominator being one order less than the numerator, i.e., $m = n - 1$ , which is of cubic order, $n = 3$ . The corresponding constant coefficient matrix for the Wetzel numerator $(n = 3)$ is [54] ${[{\underline{\underline{C}}}^{(3)}]}^{T} = [\begin{array}{r} 0.1433751825 & 0.1433751825 & 0.9685744898 \\ - 0.6566650339 & - 0.5209453949 & - 2.5526857671 \\ - 0.5106016916 & - 0.6463213306 & - 2.5756669706 \\ 3.5295952199 & 0.6031924921 & 2.2044050704 \\ 4.4349137241 & 2.3303190917 & 4.4520903005 \\ 0.1229618909 & 5.1539592511 & 2.2485545147 \\ - 2.9144388828 & - 0.2256222796 & - 0.6202937932 \\ - 5.5556896198 & - 1.6481269200 & - 1.8811803355 \\ - 2.8284365891 & - 5.4494528976 & - 1.9023485762 \\ 0.2292109036 & - 3.7461520908 & - 0.6414620339 \end{array}]$ And for the denominator $(m = n - 1 = 2)$ , ${[{\underline{\underline{C}}}^{(2)}]}^{T} = [\begin{array}{r} 1.0000000000 & 1.0000000000 & 1.0000000000 \\ 0.7257989503 & 0.6916858207 & - 1.2134964928 \\ 3.0941511876 & 3.1282643172 & - 1.2128608265 \\ - 1.6239324646 & - 1.5898193351 & 0.2393747647 \\ - 4.7303686308 & - 4.7303686308 & 0.6004510415 \\ - 3.1742364608 & - 3.2083495904 & 0.2162486576 \end{array}]$ Mullens M. [40] developed several high order polynomial fitted closures for short fiber composites and introduced the time derivative fitted closures. For the LAR32 closure by Mullens M. [40], the corresponding constant coefficient matrix for the numerator $(n = 3)$ is ${[{\underline{\underline{C}}}^{(3)}]}^{T} = [\begin{array}{r} 0.087602233 & 0.156805152 & 1.072423739 \\ 0.028205550 & - 0.577818864 & - 2.803554028 \\ - 0.426784335 & - 0.514280920 & - 2.661576129 \\ 1.274677110 & 0.684250887 & 2.389379765 \\ 0.876469059 & 2.132305029 & 4.566728489 \\ 0.602031647 & 3.454835266 & 2.097523143 \\ - 1.066583115 & - 0.263237143 & - 0.658248930 \\ - 1.918931146 & - 1.614122610 & - 1.904704744 \\ - 0.934291306 & - 4.005261132 & - 1.754978355 \\ - 0.262854903 & - 2.228133231 & - 0.508282668 \end{array}]$ And for the denominator $(m = n - 1 = 2)$ , ${[{\underline{\underline{C}}}^{(2)}]}^{T} = [\begin{array}{r} 1.000000000 & 1.000000000 & 1.000000000 \\ - 0.244001948 & 0.365652907 & - 1.068512526 \\ - 0.574150861 & 1.385725477 & - 0.771356469 \\ - 0.432097367 & - 1.359687152 & 0.067386858 \\ - 0.895226091 & - 2.866357848 & 0.206908269 \\ - 0.462709527 & - 1.518996192 & - 0.248999874 \end{array}]$

Quartic Orthotropic Fitted Closures $(n = 4)$

The constant coefficient matrix ${\underline{\underline{C}}}^{'}$ based on a regression analysis by Verweyst [53] developed from Carlson elliptic integrals is ${[{\underline{\underline{C}}}^{(4)}]}^{T} = [\begin{array}{r} 0.6363 & 0.6363 & 2.7405 \\ - 1.8727 & - 3.3153 & - 9.1220 \\ - 4.4797 & - 3.0371 & - 12.2571 \\ 11.9590 & 11.8273 & 34.3199 \\ 3.8446 & 6.8815 & 13.8295 \\ 11.3421 & 8.4368 & 25.8685 \\ - 10.9583 & - 15.9121 & - 37.7029 \\ - 20.7278 & - 15.1516 & - 50.2756 \\ - 2.1162 & - 6.4873 & - 10.8802 \\ - 12.3876 & - 8.6389 & - 26.9637 \\ 9.8160 & 9.3252 & 27.3347 \\ 3.4790 & 7.7468 & 15.2651 \\ 11.7493 & 7.4815 & 26.1135 \\ 0.5080 & 2.2848 & 3.4321 \\ 4.8837 & 3.5977 & 10.6117 \end{array}]$ The constant coefficient matrix ${\underline{\underline{C}}}^{'}$ based on a regression analysis for the FFLAR4 closure by Mullens [40] is ${[{\underline{\underline{C}}}^{(4)}]}^{T} = [\begin{array}{r} 0.678225884 & 0.748226727 & 3.167356369 \\ - 3.834359034 & - 4.249612053 & - 13.288266400 \\ - 2.664862865 & - 2.987266447 & - 11.680179330 \\ 9.746185193 & 8.641488072 & 23.788431340 \\ 14.209962670 & 14.938209410 & 43.700607680 \\ 2.700369681 & 5.974489008 & 17.383121430 \\ - 8.013024236 & - 7.521216405 & - 19.959054610 \\ - 22.447252700 & - 21.757217160 & - 58.354308000 \\ - 13.078649640 & - 15.798676320 & - 49.513705640 \\ - 0.125467689 & - 3.616551654 & - 11.755525930 \\ 2.417857515 & 2.376441613 & 6.291273472 \\ 10.563248410 & 10.222185780 & 25.844317920 \\ 12.689484570 & 12.640352670 & 35.425354130 \\ 2.487386515 & 4.788201652 & 18.226443930 \\ - 0.328195677 & 1.056519961 & 2.925785795 \end{array}]$ The constant coefficient matrix ${\underline{\underline{C}}}^{'}$ based on a regression analysis for the LAR4 closure by Mullens [40] is ${[{\underline{\underline{C}}}^{(4)}]}^{T} = [\begin{array}{r} 0.813175172 & 1.768619587 & 4.525066937 \\ - 3.065410883 & - 9.826017151 & - 19.259137620 \\ - 4.659333003 & - 6.484058476 & - 17.650178090 \\ 6.329870878 & 19.986994700 & 33.901239610 \\ 14.747639770 & 28.905936750 & 61.543979540 \\ 9.739797775 & 10.759963010 & 28.467355970 \\ - 4.216519964 & - 17.715409270 & - 27.768082700 \\ - 15.922240910 & - 40.492387100 & - 76.738638810 \\ - 20.818571900 & - 27.442217500 & - 68.977583290 \\ - 8.993993112 & - 7.230748101 & - 22.399036130 \\ 1.138888034 & 5.785725498 & 8.600822308 \\ 5.834142985 & 18.709047480 & 32.480679940 \\ 11.470974520 & 19.729631240 & 43.875135630 \\ 9.874209286 & 8.882877701 & 26.928320210 \\ 3.100457733 & 2.224834058 & 7.101978254 \end{array}]$

Appendix B.1.2. Invariant-Based Optimal Fitting Closure (IBOF) Approximation

The unknown independent coefficients of the binomial expansion for the six parameters in the IBOF closure representation based on regression fitting by Chung et al. [34] of actual flow data obtained from the distribution function closure considering different flow types like EBF closures are given in Table A1 below.

Table A1

The 5th order binomial fitting coefficients for the IBOF closure approximation.

$k \ m$	1	2	3
0	2.49409081657860 × 10¹	−4.97217790110754 × 10⁻¹	2.34146291570999 × 10¹
1	−4.35101153160329 × 10²	2.34980797511405 × 10¹	−4.12048043372534 × 10²
2	7.03443657916476 × 10³	1.53965820593506 × 10²	5.73259594331015 × 10³
3	3.72389335663877 × 10³	−3.91044251397838 × 10²	3.19553200392089 × 10³
4	−1.33931929894245 × 10⁵	−2.13755248785646 × 10³	−6.05006113515592 × 10⁴
5	8.23995187366106 × 10⁵	1.52772950743819 × 10⁵	−4.85212803064813 × 10⁴
6	−1.59392396237307 × 10⁴	2.96004865275814 × 10³	−1.10656935176569 × 10⁴
7	8.80683515327916 × 10⁵	−4.00138947092812 × 10³	−4.77173740017567 × 10⁴
8	−9.91630690741981 × 10⁶	−1.85949305922308 × 10⁶	5.99066486689836 × 10⁶
9	8.00970026849796 × 10⁶	2.47717810054366 × 10⁶	−4.60543580680696 × 10⁷
10	3.22219416256417 × 10⁴	−1.04092072189767 × 10⁴	1.28967058686204 × 10⁴
11	−2.37010458689252 × 10⁶	1.01013983339062 × 10⁵	2.03042960322874 × 10⁶
12	3.79010599355267 × 10⁷	7.32341494213578 × 10⁶	−5.56606156734835 × 10⁷
13	−3.37010820273821 × 10⁷	−1.47919027644202 × 10⁷	5.67424911007837 × 10⁸
14	−2.57258805870567 × 10⁸	−6.35149929624336 × 10⁷	−1.52752854956514 × 10⁹
15	−2.32153488525298 × 10⁴	1.38088690964946 × 10⁴	4.66767581292985 × 10³
16	2.14419090344474 × 10⁶	−2.47435106210237 × 10⁵	−4.99321746092534 × 10⁶
17	−4.49275591851490 × 10⁷	−9.02980378929272 × 10⁶	1.32124828143333 × 10⁸
18	−2.13133920223355 × 10⁷	7.24969796807399 × 10⁶	−1.62359994620983 × 10⁹
19	1.57076702372204 × 10⁹	4.87093452892595 × 10⁸	7.92526849882218 × 10⁹
20	−3.95769398304473 × 10⁹	−1.60162178614234 × 10⁹	−1.28050778279459 × 10¹⁰

References

1. Jeffery, G.B. The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. Lond. Ser. Contain. Pap. Math. Phys. Character; 1922; 102, pp. 161-179. [DOI: https://dx.doi.org/10.1098/rspa.1922.0078]

2. Saffman, P.G. On the motion of small spheroidal particles in a viscous liquid. J. Fluid Mech.; 1956; 1, pp. 540-553. [DOI: https://dx.doi.org/10.1017/S0022112056000354]

3. Folgar, F.; Tucker, C.L. Orientation Behavior of Fibers in Concentrated Suspensions. J. Reinf. Plast. Compos.; 1984; 3, pp. 98-119. [DOI: https://dx.doi.org/10.1177/073168448400300201]

4. Folgar Portillo, F. Fiber Orientation Distribution in Concentrated Suspensions: A Predictive Model. Ph.D. Dissertation/Thesis; University of Illinois at Urbana: Champaign, IL, USA, 1983; Available online: http://ezproxy.baylor.edu/login?url=https://www.proquest.com/dissertations-theses/fiber-orientation-distribution-concentrated/docview/303264729/se-2?accountid=7014 (accessed on 9 March 2023).

5. Bay, R.S. Fiber Orientation in Injection-Molded Composites: A Comparison of Theory and Experiment. Ph.D. Dissertation/Thesis; University of Illinois at Urbana: Champaign, IL, USA, 1991; Available online: http://ezproxy.baylor.edu/login?url=https://www.proquest.com/dissertations-theses/fiber-orientation-injection-molded-composites/docview/303941981/se-2?accountid=7014 (accessed on 22 December 2021).

6. Montgomery-Smith, S.; Jack, D.; Smith, D.E. The Fast Exact Closure for Jeffery’s equation with diffusion. J. Non-Newton. Fluid Mech.; 2011; 166, pp. 343-353. [DOI: https://dx.doi.org/10.1016/j.jnnfm.2010.12.010]

7. Advani, S.G.; Tucker, C.L. The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites. J. Rheol.; 1987; 31, pp. 751-784. [DOI: https://dx.doi.org/10.1122/1.549945]

8. Huynh, H.M. Improved Fiber Orientation Predictions for Injection-Molded Composites. Master’s Thesis; University of Illinois at Urbana: Champaign, IL, USA, 2001.

9. Wang, J.; O’Gara, J.F.; Tucker, C.L. An objective model for slow orientation kinetics in concentrated fiber suspensions: Theory and rheological evidence. J. Rheol.; 2008; 52, pp. 1179-1200. [DOI: https://dx.doi.org/10.1122/1.2946437]

10. Tseng, H.-C.; Chang, R.-Y.; Hsu, C.-H. Phenomenological improvements to predictive models of fiber orientation in concentrated suspensions. J. Rheol.; 2013; 57, pp. 1597-1631. [DOI: https://dx.doi.org/10.1122/1.4821038]

11. Tseng, H.-C.; Chang, R.-Y.; Hsu, C.-H. Method and Computer Readable Media for Determining Orientation of Fibers in a Fluid. U.S. Patent; 8571828 B2, 29 August 2023.

12. Phelps, J.H.; Tucker, C.L. An anisotropic rotary diffusion model for fiber orientation in short- and long-fiber thermoplastics. J. Non-Newton. Fluid Mech.; 2009; 156, pp. 165-176. [DOI: https://dx.doi.org/10.1016/j.jnnfm.2008.08.002]

13. Ranganathan, S.; Advani, S.G. Fiber–fiber interactions in homogeneous flows of nondilute suspensions. J. Rheol.; 1991; 35, pp. 1499-1522. [DOI: https://dx.doi.org/10.1122/1.550244]

14. Fan, X.; Phan-Thien, N.; Zheng, R. A direct simulation of fibre suspensions. J. Non-Newton. Fluid Mech.; 1998; 74, pp. 113-135. [DOI: https://dx.doi.org/10.1016/S0377-0257(97)00050-5]

15. Phan-Thien, N.; Fan, X.-J.; Tanner, R.I.; Zheng, R. Folgar–Tucker constant for a fibre suspension in a Newtonian fluid. J. Non-Newton. Fluid Mech.; 2002; 103, pp. 251-260. [DOI: https://dx.doi.org/10.1016/S0377-0257(02)00006-X]

16. Koch, D.L. A model for orientational diffusion in fiber suspensions. Phys. Fluids; 1995; 7, pp. 2086-2088. [DOI: https://dx.doi.org/10.1063/1.868455]

17. Hand, G.L. A theory of anisotropic fluids. J. Fluid Mech.; 1962; 13, pp. 33-46. [DOI: https://dx.doi.org/10.1017/S0022112062000476]

18. Tseng, H.-C.; Chang, R.-Y.; Hsu, C.-H. An objective tensor to predict anisotropic fiber orientation in concentrated suspensions. J. Rheol.; 2016; 60, pp. 215-224. [DOI: https://dx.doi.org/10.1122/1.4939098]

19. Tseng, H.-C.; Chang, R.-Y.; Hsu, C.-H. The use of principal spatial tensor to predict anisotropic fiber orientation in concentrated fiber suspensions. J. Rheol.; 2018; 62, pp. 313-320. [DOI: https://dx.doi.org/10.1122/1.4998520]

20. Bakharev, A.Y.; Yu, H.; Ray, S.; Speight, R.; Wang, J. Using new anisotropic rotational diffusion model to improve prediction of short fibers in thermoplastic injection molding. Proceedings of the 2018 Society of Plastics Engineers (SPE) International Polyolefins Conference; Lubbock, TX, USA, 25–28 February 2018.

21. Latz, A.; Strautins, U.; Niedziela, D. Comparative numerical study of two concentrated fiber suspension models. J. Non-Newton. Fluid Mech.; 2010; 165, pp. 764-781. [DOI: https://dx.doi.org/10.1016/j.jnnfm.2010.04.001]

22. Kugler, S.K.; Kech, A.; Cruz, C.; Osswald, T. Fiber Orientation Predictions—A Review of Existing Models. J. Compos. Sci.; 2020; 4, 69. [DOI: https://dx.doi.org/10.3390/jcs4020069]

23. Favaloro, A.J.; Tucker, C.L. Analysis of anisotropic rotary diffusion models for fiber orientation. Compos. Part Appl. Sci. Manuf.; 2019; 126, 105605. [DOI: https://dx.doi.org/10.1016/j.compositesa.2019.105605]

24. Agboola, B.O.; Jack, D.A.; Montgomery-Smith, S. Effectiveness of recent fiber-interaction diffusion models for orientation and the part stiffness predictions in injection molded short-fiber reinforced composites. Compos. Part A Appl. Sci. Manuf.; 2012; 43, pp. 1959-1970. [DOI: https://dx.doi.org/10.1016/j.compositesa.2012.06.015]

25. Park, J.M.; Park, S.J. Modeling and Simulation of Fiber Orientation in Injection Molding of Polymer Composites. Math. Probl. Eng.; 2011; 2011, 105637. [DOI: https://dx.doi.org/10.1155/2011/105637]

26. Bay, R.S.; Tucker, C.L. Fiber orientation in simple injection moldings. Part I: Theory and numerical methods. Polym. Compos.; 1992; 13, pp. 317-331. [DOI: https://dx.doi.org/10.1002/pc.750130409]

27. Bay, R.S.; Tucker, C.L. Fiber orientation in simple injection moldings. Part II: Experimental results. Polym. Compos.; 1992; 13, pp. 332-341. [DOI: https://dx.doi.org/10.1002/pc.750130410]

28. Gupta, M.; Wang, K.K. Fiber orientation and mechanical properties of short-fiber-reinforced injection-molded composites: Simulated and experimental results. Polym. Compos.; 1993; 14, pp. 367-382. [DOI: https://dx.doi.org/10.1002/pc.750140503]

29. Chung, S.T.; Kwon, T.H. Coupled analysis of injection molding filling and fiber orientation, including in-plane velocity gradient effect. Polym. Compos.; 1996; 17, pp. 859-872. [DOI: https://dx.doi.org/10.1002/pc.10679]

30. Vincent, M.; Devilers, E.; Agassant, J.-F. Fibre orientation calculation in injection moulding of reinforced thermoplastics. J. Non-Newton. Fluid Mech.; 1997; 73, pp. 317-326. [DOI: https://dx.doi.org/10.1016/S0377-0257(97)00048-7]

31. Vincent, M.; Giroud, T.; Clarke, A.; Eberhardt, C. Description and modeling of fiber orientation in injection molding of fiber reinforced thermoplastics. Polymer; 2005; 46, pp. 6719-6725. [DOI: https://dx.doi.org/10.1016/j.polymer.2005.05.026]

32. VerWeyst, B.E.; Tucker, C.L.; Foss, P.H.; O’Gara, J.F. Fiber Orientation in 3-D Injection Molded Features: Prediction and Experiment. Int. Polym. Process.; 1999; 14, pp. 409-420. [DOI: https://dx.doi.org/10.3139/217.1568]

33. Hinch, E.J.; Leal, L.G. Constitutive equations in suspension mechanics. Part 2. Approximate forms for a suspension of rigid particles affected by Brownian rotations. J. Fluid Mech.; 1976; 76, pp. 187-208. [DOI: https://dx.doi.org/10.1017/S0022112076003200]

34. Chung, D.H.; Kwon, T.H. Invariant-based optimal fitting closure approximation for the numerical prediction of flow-induced fiber orientation. J. Rheol.; 2002; 46, pp. 169-194. [DOI: https://dx.doi.org/10.1122/1.1423312]

35. Jack, D.A.; Smith, D.E. The effect of fibre orientation closure approximations on mechanical property predictions. Compos. Part A Appl. Sci. Manuf.; 2007; 38, pp. 975-982. [DOI: https://dx.doi.org/10.1016/j.compositesa.2006.06.016]

36. Quintana, M.C.; Frontini, P.M.; Arriaga, A.; Plank, B.; Major, Z. Fiber Orientation Distribution Predictions for an Injection Molded Venturi-Shaped Part Validated Against Experimental Micro-Computed Tomography Characterization. Front. Mater.; 2020; 7, 169. [DOI: https://dx.doi.org/10.3389/fmats.2020.00169]

37. Jack, D.A. Advanced Analysis of Short-Fiber Polymer Composite Material Behavior. Ph.D. Thesis; University of Missouri: Columbia, MO, USA, 2006.

38. Kuzmin, D. Planar Orthotropic Closures for Orientation Tensors in Fiber Suspension Flow Models. SIAM J. Appl. Math.; 2018; 78, pp. 3040-3059. [DOI: https://dx.doi.org/10.1137/18M1175665]

39. Cintra, J.S.; Tucker, C.L. Orthotropic closure approximations for flow-induced fiber orientation. J. Rheol.; 1995; 39, pp. 1095-1122. [DOI: https://dx.doi.org/10.1122/1.550630]

40. Mullens, M. Developing New Fitted Closure Approximations for Short-Fiber Reinforced Polymer Composites. Master’s Thesis; University of Missouri: Columbia, MO, USA, 2010; Available online: https://core.ac.uk/download/pdf/62781195.pdf (accessed on 9 March 2023).

41. Jack, D.A.; Schache, B.; Smith, D.E. Neural network-based closure for modeling short-fiber suspensions. Polym. Compos.; 2010; 31, pp. 1125-1141. [DOI: https://dx.doi.org/10.1002/pc.20912]

42. Qadir, N.U.; Jack, D.A. Modeling fibre orientation in short fibre suspensions using the neural network-based orthotropic closure. Compos. Part A Appl. Sci. Manuf.; 2009; 40, pp. 1524-1533. [DOI: https://dx.doi.org/10.1016/j.compositesa.2009.06.010]

43. Jack, D.A. Investigating the Use of Tensors in Numerical Predictions for Short-Fiber Reinforced Polymer Composites. Ph.D. Thesis; University of Missouri: Columbia, MO, USA, 2003.

44. Chapra, S.C. Applied Numerical Methods with MATLAB for Engineers and Scientists; 3rd ed. McGraw-Hill: New York, NY, USA, 2012.

45. Tucker, C.L. Fundamentals of Fiber Orientation: Description, Measurement and Prediction; Hanser Publishers: Munich, Germany, 2022; [DOI: https://dx.doi.org/10.3139/9781569908761]

46. Eberle, A.P.R.; Ortman, K.; Baird, D.G. Structure and Rheology of Fiber Suspensions. Applied Polymer Rheology; 1st ed. Kontopoulou, M. Wiley: Hoboken, NJ, USA, 2011; pp. 113-151. [DOI: https://dx.doi.org/10.1002/9781118140611.ch4]

47. Zhang, D. Flow-Induced Micro- and Nano-Fiber Suspensions in Short-Fiber Reinforced Composite Materials Processing. Ph.D. Thesis/Dissertation; University of Missouri: Columbia, MO, USA, 2013; Available online: http://ezproxy.baylor.edu/login?url=https://www.proquest.com/dissertations-theses/flow-induced-micro-nano-fiber-suspensions-short/docview/2202998717/se-2?accountid=7014 (accessed on 8 December 2020).

48. Awenlimobor, A.; Smith, D.E. Effect of shear-thinning rheology on the dynamics and pressure distribution of a single rigid ellipsoidal particle in viscous fluid flow. Phys. Fluids; 2024; 36, 123113. [DOI: https://dx.doi.org/10.1063/5.0242953]

49. Smith, D.E.; Siddhi, V. Generalized Approach for Incorporating Normalization Conditions in Design Sensitivity Analysis of Eigenvectors. AIAA J.; 2006; 44, pp. 2552-2561. [DOI: https://dx.doi.org/10.2514/1.18760]

50. Wang, J. Moldflow Research & Development. Online. 2018; Available online: https://damassets.autodesk.net/content/dam/autodesk/www/campaigns/autodesk-moldflow-summit-recording/7_R&D-Update%20-%20Jin%20Wang.pdf (accessed on 9 March 2023).

51. Doi, M. Molecular dynamics rheological properties of concentrated solutions of rodlike polymers in isotropic liquid crystalline phases. J. Polym. Sci. Polym. Phys. Ed.; 1981; 19, pp. 229-243. [DOI: https://dx.doi.org/10.1002/pol.1981.180190205]

52. Lipscomb, G.G.; Denn, M.M.; Hur, D.U.; Boger, D.V. The flow of fiber suspensions in complex geometries. J. Non-Newton. Fluid Mech.; 1988; 26, pp. 297-325. [DOI: https://dx.doi.org/10.1016/0377-0257(88)80023-5]

53. VerWeyst, B.E. Numerical Predictions of Flow-Induced Fiber Orientation in Three-Dimensional Geometries. Ph.D. Dissertation/Thesis; University of Illinois at Urbana: Champaign, IL, USA, 1998; Available online: http://ezproxy.baylor.edu/login?url=https://www.proquest.com/dissertations-theses/numerical-predictions-flow-induced-fiber/docview/304440327/se-2?accountid=7014 (accessed on 18 February 2022).

54. Wetzel, E.D.; Tucker, C.L. Area tensors for modeling microstructure during laminar liquid-liquid mixing. Int. J. Multiph. Flow; 1999; 25, pp. 35-61. [DOI: https://dx.doi.org/10.1016/S0301-9322(98)00013-5]

55. Verleye, V.; Couniot, A.; Dupret, F. Numerical prediction of fibre orientation in complex injection-moulded parts. WIT Trans. Eng. Sci.; 1994; 4, 10. [DOI: https://dx.doi.org/10.2495/CP940341]

56. Chung, D.H.; Kwon, T.H. Improved model of orthotropic closure approximation for flow induced fiber orientation. Polym. Compos.; 2001; 22, pp. 636-649. [DOI: https://dx.doi.org/10.1002/pc.10566]

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Determination of the Steady State Fiber Orientation Tensor States in Homogeneous Flows with Newton–Raphson Iteration Using Exact Jacobians

Content area

Abstract

Full text

1. Introduction

2. Materials and Methods

2.1. Representation of Fiber Orientation State

2.2. Determination of Steady State Orientations

2.3. Fiber Orientation Modelling in the Dilute Regime

2.4. Fiber Orientation Modelling in Semi-Dilute and Concentrated Regime

2.4.1. The Folgar-Tucker Diffusion Model

2.4.2. Strain Reduction Factor (SRF) Model

2.4.3. Reduced Strain Closure (RSC) Model

2.4.4. Retarding Principal Rate (RPR) Model

2.4.5. Anisotropic Rotary Diffusion (ARD) Models

2.4.6. Nematic Potential (NEM) Model

2.5. Closure Approximations and Their Explicit Derivatives

2.5.1. Quadratic Closure Approximation

2.5.2. Linear Closure Approximation

2.5.3. Hybrid Closure Approximation

2.5.4. Hinch and Leal Closure Approximation

2.5.5. Eigenvalue-Based Fitted (EBF) Closure Approximations

2.5.6. Invariant-Based Fitted (IBF) Closure Approximations

2.6. Error Estimate

3. Results

3.1. Validation of Exact Jacobians Based on Finite Difference Approximation

3.2. Comparison of the NR Method and the Dormand-Prince Runge–Kutta (RK45) Method

3.3. Performance of NR Method Using Various Closure Approximations

3.4. Homogenous Flow Considerations

4. Conclusions