# The Fast Exact Closure for Jeffery’s Equation with Diffusion

Stephen Montgomery-Smith  David Jack111Corresponding Author. Email:
<EMAIL_ADDRESS>Phone: 254-710-3347. Douglas E. Smith  Department of
Mathematics, University of Missouri, Columbia MO 65211, U.S.A.  Department of
Mechanical Engineering, Baylor University, Waco TX 76798, U.S.A.  Department
of Mechanical and Aerospace Engineering, University of Missouri, Columbia MO
65211, U.S.A.

###### Abstract

Jeffery’s equation with diffusion is widely used to predict the motion of
concentrated fiber suspensions in flows with low Reynold’s numbers.
Unfortunately, the evaluation of the fiber orientation distribution can
require excessive computation, which is often avoided by solving the related
second order moment tensor equation. This approach requires a ‘closure’ that
approximates the distribution function’s fourth order moment tensor from its
second order moment tensor. This paper presents the _Fast Exact Closure_ (FEC)
which uses conversion tensors to obtain a pair of related ordinary
differential equations; avoiding approximations of the higher order moment
tensors altogether. The FEC is exact in that when there are no fiber
interactions, it exactly solves Jeffery’s equation. Numerical examples for
dense fiber suspensions are provided with both a Folgar-Tucker (1984)
diffusion term and the recent anisotropic rotary diffusion term proposed by
Phelps and Tucker (2009). Computations demonstrate that the FEC exhibits
improved accuracy with computational speeds equivalent to or better than
existing closure approximations.

###### keywords:

B. Directional orientation, B. Rheological properties, D. Injection molding,
Jeffery’s equation with rotary diffusion

##  1 Introduction

The industrial demand has continued to increase for high-strength, low-weight,
rapid production parts such as those made of short discontinuous fiber
composites with injection molding processes. For effective design, it is
essential to understand the dependance of the final part performance of short-
fiber injection molded composites with the variations in the microstructure
due to the processing (see e.g. [1, 2]). The Folgar and Tucker model of
isotropic diffusion [3] for fiber interactions within a suspension has been
used for several decades to compute fiber orientation and has been implemented
to some extent within most related industrial and research computer
simulations. Unfortunately, direct computations of the isotropic diffusion
model are computationally prohibitive, and most implementations employ the
orientation tensor approach of Advani and Tucker [4] where the moments of the
fiber orientation are solved, thus indirectly quantifying the fiber
orientation distribution. The orientation tensor approach requires knowledge
of the next higher-order moment tensor, thus requiring some form of a closure.
The hybrid closure of Advani and Tucker [4] has been used extensively due to
its computational efficiencies, but in implementation it will overpredict the
alignment state in simple shear flow [5]. Cintra and Tucker [6] introduced the
class of the orthotropic closures, which result in significant accuracy
improvements when compared to the hybrid closure, but at an increase in
computational costs.

With recent advances in part repeatability, the limitation of the isotropic
diffusion model has become apparent [7]. Recent anisotropic diffusion models
[8, 9, 10, 11] propose new forms with greater accuracies for modeling fiber
collisions, but these anisotropic diffusion models pose a new set of
computational complications. In particular is the concern that nearly all of
the fitted orthotropic closures are obtained by fitting orientation
information based on direct numerical solutions of the Folgar-Tucker diffusion
model. The exception is the orthotropic closures of Wetzel [12] and VerWeyst
[13] which were both constructed on distributions formed through the elliptic
integral form for orientations encompassing the eigenspace [6].

The Exact Closure of Montgomery-Smith et al. [14] presents an alternative to
the classical closure form, and provides an exact solution for pure Jeffery’s
motion (i.e., the dilute regime). The Exact Closure avoids the curve fitting
process required to define fitted closures, by solving a set of related ODEs
of the fiber orientation. In the present paper, we extend the Exact Closure
form to systems of concentrated suspensions that are more relevant to modeling
the processing of short-fiber composites. Furthermore, we introduce the new
_Fast Exact Closure_ (FEC) that defines conversion tensors that lead to a
coupled system of ordinary differential equations that avoid costly closure
computations. The FEC form is derived for fiber collision models for both the
isotropic diffusion model of Folgar and Tucker and the recent anisotropic
diffusion model of Phelps and Tucker [9]. Results presented will demonstrate
the effectiveness of this alternative approach for modeling fiber orientation,
both for accuracy and for computational speed.

##  2 Fiber Motion Basics

Jeffery’s equation [15] has been used to predict the motion of the direction
of axi-symmetric fibers under the influence of a low Reynold’s number flow of
a Newtonian fluid, whose velocity field is
${\mathbf{u}}={\mathbf{u}}({\mathbf{x}},t)$. The directions of the fibers is
represented by the fiber orientation distribution
$\psi=\psi({\mathbf{x}},{\mathbf{p}},t)$, where ${\mathbf{p}}$ is an element
of the orientation space, that is, the 2-dimensional sphere
$S=\\{{\mathbf{p}}=(p_{1},p_{2},p_{3}):p_{1}^{2}+p_{2}^{2}+p_{3}^{2}=1\\}$.
Thus given a subset $E$ of $S$, the proportion of fibers whose direction is in
$E$ is given by $\int_{E}\psi({\mathbf{x}},{\mathbf{p}},t)\,d{\mathbf{p}}$,
where $d{\mathbf{p}}$ represents the usual integration over $S$. In
particular, an isotropic distribution is represented by $\psi=1/4\pi$. The
Jeffery’s equation for the fiber orientation distribution is

$\frac{D\psi}{Dt}=-\tfrac{1}{2}{\boldsymbol{\nabla}}_{\mathbf{p}}\cdot((\Omega\cdot{\mathbf{p}}+\lambda\Gamma\cdot{\mathbf{p}}-\lambda\Gamma:{\mathbf{p}}{\mathbf{p}}{\mathbf{p}})\psi)$
(1)

Here $\Omega$ is the vorticity, that is, the anti-symmetric part
${\boldsymbol{\nabla}}{\mathbf{u}}-({\boldsymbol{\nabla}}{\mathbf{u}})^{T}$ of
the Jacobian of the velocity field
${\boldsymbol{\nabla}}{\mathbf{u}}=(\partial u_{i}/\partial x_{j})_{1\leq
i,j\leq 3}$, and $\Gamma$ is the rate of strain tensor, that is, the symmetric
part
${\boldsymbol{\nabla}}{\mathbf{u}}+({\boldsymbol{\nabla}}{\mathbf{u}})^{T}$ of
the Jacobean of the velocity field. Also, $D/Dt=\partial/\partial
t+{\mathbf{u}}\cdot{\boldsymbol{\nabla}}$ represents the material derivative,
and
${\boldsymbol{\nabla}}_{\mathbf{p}}=(I-{\mathbf{p}}{\mathbf{p}})\cdot\left(\tfrac{\partial}{\partial
p_{1}},\tfrac{\partial}{\partial p_{2}},\tfrac{\partial}{\partial
p_{3}}\right)$ is the gradient operator restricted to the sphere.

Equation (1) is modified to incorporate the rotary diffusion expressed by Bird
et al. [16], occasionally referred to as the generalized Fokker-Planck or the
Smoluchowski equation [17], as

$\frac{D\psi}{Dt}=-\tfrac{1}{2}{\boldsymbol{\nabla}}_{\mathbf{p}}\cdot((\Omega\cdot{\mathbf{p}}+\lambda\Gamma\cdot{\mathbf{p}}-\lambda\Gamma:{\mathbf{p}}{\mathbf{p}}{\mathbf{p}})\psi)+\Delta_{\mathbf{p}}(D_{r}\psi),$
(2)

where $D_{r}$ captures the effect of fiber interaction and depends upon the
flow kinetics. Here
$\Delta_{\mathbf{p}}={\boldsymbol{\nabla}}_{\mathbf{p}}\cdot{\boldsymbol{\nabla}}_{\mathbf{p}}$
represents the Beltrami-Laplace operator on the sphere. Folgar and Tucker [3]
selected $D_{r}=C_{I}\dot{\gamma}$ where
$\dot{\gamma}=\left(\frac{1}{2}\Gamma:\Gamma\right)^{1/2}$ and $C_{I}$ is a
constant that depends upon the volume fraction and aspect ratio of the fibers.

Other authors have considered a wider class of diffusion terms. For example,
Koch [10], and Phelps and Tucker [9] considered anisotropic diffusion

$\frac{D\psi}{Dt}=-\tfrac{1}{2}{\boldsymbol{\nabla}}_{\mathbf{p}}\cdot((\Omega\cdot{\mathbf{p}}+\lambda\Gamma\cdot{\mathbf{p}}-\lambda\Gamma:{\mathbf{p}}{\mathbf{p}}{\mathbf{p}})\psi)+{\boldsymbol{\nabla}}_{\mathbf{p}}\cdot(I-{\mathbf{p}}{\mathbf{p}})\cdot
D_{r}\cdot{\boldsymbol{\nabla}}_{\mathbf{p}}\psi$ (3)

where $D_{r}$ is the anisotropic diffusion matrix, calculated as a function of
$\psi$ and ${\boldsymbol{\nabla}}u$ (see, e.g., [9, 10]).

Since these are, in effect, partial differential equations in 5-spacial
dimensions (3 for space and 2 for the orientation defined on a unit sphere),
numerically calculating solutions can be rather daunting with solutions taking
days to weeks for simple flows. Hence Hinch and Leal [18] suggested to recast
the equation in terms of moment tensors. For example, the second and fourth
moment tensors are defined by

$\displaystyle
A=\int_{S}{\mathbf{p}}{\mathbf{p}}\psi\,d{\mathbf{p}},\qquad\qquad\mathbb{A}=\int_{S}{\mathbf{p}}{\mathbf{p}}{\mathbf{p}}{\mathbf{p}}\psi\,d{\mathbf{p}}$
(4)

Then Jeffery’s equation (1) for the second order moment tensor can be
expressed as

$\frac{DA}{Dt}=\tfrac{1}{2}(\Omega\cdot A-A\cdot\Omega+\lambda(\Gamma\cdot
A+A\cdot\Gamma)-2\lambda\mathbb{A}:\Gamma)$ (5)

and the equations (2) and (3) with diffusion terms become

$\frac{DA}{Dt}=\tfrac{1}{2}(\Omega\cdot A-A\cdot\Omega+\lambda(\Gamma\cdot
A+A\cdot\Gamma)-2\lambda\mathbb{A}:\Gamma)+\mathcal{D}[A]$ (6)

where $\mathcal{D}[A]$ for isotropic diffusion as expressed in equation (2)
becomes

$\mathcal{D}[A]=D_{r}(2I-6A)$ (7)

and subsequently the anisotropic diffusion of equation (3) (see [9]) is

$\mathcal{D}[A]=2D_{r}-2(\mathop{\text{tr}}D_{r})A-5(A\cdot D_{r}+D_{r}\cdot
A)+10\mathbb{A}:D_{r}$ (8)

The difficulty with equations (5) and (6) is that they explicitly include the
fourth order moment tensor, and implicitly the higher order diffusion models
of equation (8) include moments higher than the second-moment. To circumvent
this problem, various authors (for example, [18, 19, 20, 21, 6, 1, 22, 23,
24]) have proposed _closures_ , that is, formulae to calculate the fourth
order moment tensor $\mathbb{A}$ from the second order moment tensor $A$. The
mapping from $A$ to ${\mathbb{A}}$ is not unique, thus closures are only able
to approximately obtain a higher order moment from the lower order moments.
Most closures are often constructed by obtaining the best-fit coefficients of
for a polynomial by fitting numerical data obtained by directly evaluating
equation (2) using a finite element method to solve equation (2) (for example,
Bay [25]).

##  3 The Fast Exact Closure

Verleye and Dupret [21] (see also [12, 13, 26, 27, 28]) noted that there is an
exact closure for Jeffery’s equation when the diffusion terms are not present,
_in the particular case that the fiber orientation distribution is at some
time isotropic_. This exact closure is stated explicitly in [14] for the
scenario when the suspension is dilute. For the sake of labeling, the present
closure retains the reference _Exact Closure_ , as it is exact for Jeffery’s
equation without diffusion terms.

The Exact Closure may be directly computed by solving the elliptic integral
forms presented in equation (6), where $\mathbb{A}$ is computed from $A$ using
equations (38) and (39) as derived in [14]. This approach only gives the exact
answer to equations (2) and (3) when $D_{r}=0$ and when the orientation is
isotropic at some time. Nevertheless it is reasonable to suppose that the
exact closure should give a reasonable approximation in general, even when
$D_{r}\neq 0$ as in Verweyst et al. [1, 13]. Their ORT closure is a polynomial
approximation to the Exact Closure, and as we demonstrate below, gives answers
that are virtually indistinguishable from that of the Exact Closure.

The _Fast Exact Closure_ (FEC) performs the Exact Closure in a computationally
efficient manner. A version of FEC is described in [14], but only when the
diffusion terms are absent. In this section we describe the FEC from an
implementation perspective, and leave the full derivation to the appendix.

The idea behind the FEC is the computation of two rank 4 tensors $\mathbb{C}$
and $\mathbb{D}$, defined in equations (40) and (43), respectively, which we
define as _conversion tensors_. These tensors convert between $DA/Dt$ and
$DB/Dt$ according to the formulae

$\frac{DA}{Dt}=-\mathbb{C}:\frac{DB}{Dt},\qquad\qquad\frac{DB}{Dt}=-\mathbb{D}:\frac{DA}{Dt}$
(9)

as derived in equations (51) - (53). The orientation tensor $A$ retains the
classical meaning as described in [4] and the tensor $B$ turns out to be
extremely useful for computations. $B$ appears to be a more abstract quantity
to describe the degree of orientation much like the orientation tensor. For
example, when the orientation parameter $B$ is given as $B_{ij}=\delta_{ij}$
this is analogous to saying that the orientation is isotropic, whereas when
one of the diagonal terms of $B$ goes to $0$, it indicates that the
orientation is perfectly aligned along the corresponding coordinate axis.
Montgomery-Smith et al. [14] provide a further discussion as to the meaning of
the orientation parameter $B$

What makes everything work is the formula, proven in the appendix by equation
(54), that for any matrix $M$, we have

$\mathbb{C}:(B\cdot M+M^{T}\cdot B)=(\mathop{\text{tr}}M)A+M\cdot A+A\cdot
M^{T}-2\mathbb{A}:M$ (10)

where $\mathbb{A}$ and $A$ satisfy equations (38) and (39).

The FEC present in this paper will be of the form:

$\frac{DA}{Dt}=-\mathbb{C}:F(B)+G(A),\qquad\qquad\frac{DB}{Dt}=F(B)-\mathbb{D}:G(A)$
(11)

where $F(B)$ and $G(A)$ will be given explicitly below. This is a general form
that can be applied to a the known diffusion models that fit the form of
equation (2) or (3). The conversion tensors $\mathbb{C}$ and $\mathbb{D}$ are
defined later in this section, and in the appendix we provide a more
mathematical formula for them along with a proof of the above properties. It
is important to note that $\mathbb{C}$ and $\mathbb{D}$ may be computed
directly from $A$ and $B$ in a rather fast manner, involving nothing more than
the diagonalization and inversion of three by three symmetric matrices,
general simple arithmetic, and where appropriate invoking inverse
trigonometric or inverse hyperbolic functions.

The FEC solves the coupled ODEs of (11) simultaneously. If the initial fiber
orientation is isotropic, then $A=\tfrac{1}{3}I$ and $B=I$ at $t=0$. When the
initial fiber orientation is not isotropic, then one can compute the initial
condition for $B$ from $A$ by inverting equation (38), as described in [14].

It can be shown that the matrices $A$ and $B$ remain positive definite,
simultaneously diagonalizable, and satisfy the equations
$\mathop{\text{tr}}A=\det B=1$ for all time.

For example, the FEC for the Jeffery’s equation with isotropic diffusion given
in equation (2) is given by:

$\displaystyle\frac{DA}{Dt}=\tfrac{1}{2}\mathbb{C}:[B\cdot(\Omega+\lambda\Gamma)+(-\Omega+\lambda\Gamma)\cdot
B]+D_{r}(2I-6A)$ (12)
$\displaystyle\frac{DB}{Dt}=-\tfrac{1}{2}(B\cdot(\Omega+\lambda\Gamma)+(-\Omega+\lambda\Gamma)\cdot
B)-D_{r}\mathbb{D}:(2I-6A)$ (13)

and the FEC for Jeffery’s equation with anisotropic diffusion as shown in
equation (3) is given by

$\displaystyle\frac{DA}{Dt}=\tfrac{1}{2}\mathbb{C}:[B\cdot(\Omega+\lambda\Gamma)+(-\Omega+\lambda\Gamma)\cdot
B]+2D_{r}+3(\mathop{\text{tr}}D_{r})A-5\mathbb{C}:(B\cdot D_{r}+D_{r}\cdot B)$
(14)
$\displaystyle\frac{DB}{Dt}=-\tfrac{1}{2}(B\cdot(\Omega+\lambda\Gamma)+(-\Omega+\lambda\Gamma)\cdot
B)-\mathbb{D}:(2D_{r}+3(\mathop{\text{tr}}D_{r})A)+5(B\cdot D_{r}+D_{r}\cdot
B)$ (15)

Using equation (10) it can be seen that equation (12) comes directly from
equations (6) and (7), and equation (13) comes from applying equation (43) to
equation (12). Similarly for the anisotropic diffusion model, this can be
observed for equations (14) and (15).

Notice, for equations (12) and (13) and for equations (14) and (15), that the
fourth-order orientation tensor ${\mathbb{A}}$ does not appear. The equation
of motion for the orientation is now reduced to developing the relationship
between $A$ and $B$ with that of ${\mathbb{C}}$ and ${\mathbb{D}}$. The
conversion tensors ${\mathbb{C}}$ and ${\mathbb{D}}$ are both computed with
respect to the basis of orthonormal eigenvectors of $B$. With respect to this
basis, the matrix $B$ is diagonal with entries $b_{1}$, $b_{2}$ and $b_{3}$,
and $A$ is diagonal with entries $a_{1}$, $a_{2}$ and $a_{3}$ where we
constrain $b_{1}\leq b_{2}\leq b_{3}$ which implies that $a_{1}\geq a_{2}\geq
a_{3}$.

If the eigenvalues $b_{1}$, $b_{2}$ and $b_{3}$ are not close to each other,
then $\mathbb{C}$ is the symmetric tensor calculated using the formulae from
equations (48) and (49) from the appendix

$\begin{array}[]{lll}\mathbb{C}_{1122}=\frac{a_{1}-a_{2}}{2(b_{2}-b_{1})}&&\mathbb{C}_{1111}=\tfrac{1}{2}b_{1}^{-1}-\mathbb{C}_{1122}-\mathbb{C}_{1133}\\\
\mathbb{C}_{1133}=\frac{a_{1}-a_{3}}{2(b_{3}-b_{1})}&&\mathbb{C}_{2222}=\tfrac{1}{2}b_{2}^{-1}-\mathbb{C}_{1122}-\mathbb{C}_{2233}\\\
\mathbb{C}_{2233}=\frac{a_{2}-a_{3}}{2(b_{3}-b_{2})}&&\mathbb{C}_{3333}=\tfrac{1}{2}b_{3}^{-1}-\mathbb{C}_{1133}-\mathbb{C}_{2233}\\\
\mathbb{C}_{ijkk}=0\text{ if $i\neq j\neq k$}\end{array}$ (16)

If two or more of the eigenvalues are close to each other, then these
equations can give rise to large numerical errors, or even ‘divide by zero’
exceptions. So in this situation, we use different formulae to compute
$\mathbb{C}$.

Suppose two of the eigenvalues are close to each other, for example,
$b_{1}=b_{0}+\epsilon$ and $b_{2}=b_{0}-\epsilon$, where $\epsilon$ is small.
Thus $b_{0}=\tfrac{1}{2}(b_{1}+b_{2})$ and
$\epsilon=\tfrac{1}{2}(b_{1}-b_{2})$. Define the quantity $\mathcal{I}_{n}$
from equation (50) and with equations (57) and (58) this quantity can be
expressed as

$\begin{split}\mathcal{I}_{n+1}=\frac{2n-1}{2n(b_{0}-b_{3})}\mathcal{I}_{n}-\frac{\sqrt{b_{3}}}{nb_{0}^{n}(b_{0}-b_{3})}\text{
if $n\geq 1$}\\\
\mathcal{I}_{1}=\frac{2}{\sqrt{b_{0}-b_{3}}}\cos^{-1}\left(\sqrt{\frac{b_{3}}{b_{0}}}\right)\text{
if $b_{0}>b_{3}$}\\\
\mathcal{I}_{1}=\frac{2}{\sqrt{b_{3}-b_{0}}}\cosh^{-1}\left(\sqrt{\frac{b_{3}}{b_{0}}}\right)\text{
if $b_{0}<b_{3}$}\end{split}$ (17)

Then replace the first equation of equation (16) by

$\mathbb{C}_{1122}=\tfrac{1}{4}\mathcal{I}_{3}+\tfrac{3}{8}\mathcal{I}_{5}\epsilon^{2}+O(\epsilon^{4})$
(18)

If all three of the eigenvalues are almost equal, that is $b_{1}=1+c_{1}$,
$b_{2}=1+c_{2}$, $b_{3}=1+c_{3}$ with $|c_{1}|,|c_{2}|,|c_{3}|\leq\epsilon$,
then it can be similarly shown that

$\begin{split}\mathbb{C}_{1122}&=\textstyle\frac{1}{10}-\frac{3}{28}c_{1}-\frac{3}{28}c_{2}-\frac{1}{28}c_{3}+\frac{5}{48}c_{1}^{2}+\frac{1}{8}c_{1}c_{2}+\frac{1}{24}c_{1}c_{3}+\frac{5}{48}c_{2}^{2}+\frac{1}{24}c_{2}c_{3}+\frac{1}{48}c_{3}^{2}\\\
&\phantom{={}}\textstyle-\frac{35}{352}c_{1}^{3}-\frac{45}{352}c_{1}^{2}c_{2}-\frac{15}{352}c_{1}^{2}c_{3}-\frac{45}{352}c_{1}c_{2}^{2}-\frac{9}{176}c_{1}c_{2}c_{3}\\\
&\phantom{={}}\textstyle-\frac{9}{352}c_{1}c_{3}^{2}-\frac{35}{352}c_{2}^{3}-\frac{15}{352}c_{2}^{2}c_{3}-\frac{9}{352}c_{2}c_{3}^{2}-\frac{5}{352}c_{3}^{3}+O(\epsilon^{4})\end{split}$
(19)

with similar formulae for $\mathbb{C}_{1133}$ and $\mathbb{C}_{2233}$. The
remaining entries of $\mathbb{C}$ are computed using the last four equations
from (16).

The rank 4 conversion tensor $\mathbb{D}$ given in equation (9) is defined
through equation (43) with respect to the basis of orthonormal eigenvectors of
$B$, and can be simplified to

$\begin{split}\left[\begin{smallmatrix}\mathbb{D}_{1111}&\mathbb{D}_{1122}&\mathbb{D}_{1133}\\\
\mathbb{D}_{2211}&\mathbb{D}_{2222}&\mathbb{D}_{2233}\\\
\mathbb{D}_{3311}&\mathbb{D}_{3322}&\mathbb{D}_{3333}\end{smallmatrix}\right]\quad&=\quad\left[\begin{smallmatrix}\mathbb{C}_{1111}&\mathbb{C}_{1122}&\mathbb{C}_{1133}\\\
\mathbb{C}_{2211}&\mathbb{C}_{2222}&\mathbb{C}_{2233}\\\
\mathbb{C}_{3311}&\mathbb{C}_{3322}&\mathbb{C}_{3333}\end{smallmatrix}\right]^{-1}\\\
\mathbb{D}_{ijij}=\mathbb{D}_{ijji}=\frac{1}{4\mathbb{C}_{ijij}}\text{ if
$i\neq j$}&\qquad\quad\mathbb{D}_{ijkk}=0\text{ if $i\neq j\neq
k$}\end{split}$ (20)

Note that there is no reason to suppose that $\mathbb{D}$ is completely
symmetric because in general $\mathbb{D}_{ijij}$ will not be the same as
$\mathbb{D}_{iijj}$.

In performing the numerical calculations, it is more efficient when forming
$DA/Dt$ and $DB/Dt$ from equation (11) to calculate the right hand side in the
coordinate system of the orthonormal eigenvectors of $B$, and then convert
back to the standard coordinate system when solving for $A$ and $B$.

For example, suppose $\mathbb{B}$ is any rank four tensor such that
$\mathbb{B}_{ijkk}=0$ if $i\neq j\neq k$, and
$\mathbb{B}_{ijkl}=\mathbb{B}_{jikl}=\mathbb{B}_{klij}$. Suppose also that $N$
is a symmetric matrix. Then $\mathbb{B}:N$ can be calculated by first defining
the matrices $M_{\mathbb{B}}$ and $\tilde{M}_{\mathbb{B}}$ as

$M_{\mathbb{B}}=\left[\begin{smallmatrix}\mathbb{B}_{1111}&\mathbb{B}_{1122}&\mathbb{B}_{1133}\\\
\mathbb{B}_{1122}&\mathbb{B}_{2222}&\mathbb{B}_{2233}\\\
\mathbb{B}_{1133}&\mathbb{B}_{2233}&\mathbb{B}_{3333}\end{smallmatrix}\right],\qquad\tilde{M}_{\mathbb{B}}=\left[\begin{smallmatrix}0&\mathbb{B}_{1212}&\mathbb{B}_{1313}\\\
\mathbb{B}_{1212}&0&\mathbb{B}_{2323}\\\
\mathbb{B}_{1313}&\mathbb{B}_{2323}&0\end{smallmatrix}\right]$ (21)

then decompose

$N=\text{diag}(\mathbf{n})+\tilde{N}$ (22)

where $\mathbf{n}=(N_{11},N_{22},N_{33})$, and $\tilde{N}$ is the matrix of
the off-diagonal elements of $N$. It follows that

$\mathbb{B}:N=\text{diag}(M_{\mathbb{B}}\cdot\mathbf{n})+2\tilde{M}_{\mathbb{B}}\circ\tilde{N}$
(23)

where for any matrices $U$ and $V$ we define the entrywise product (also known
as the Hadamard or Schur product) by $(U\circ V)_{ij}=U_{ij}V_{ij}$.

###  3.1 The Reduced Strain Closure

Wang et al. [8] described a method that slows down the rate of alignment of
the fibers, which the paper calls the reduced strain closure model (RSC). The
method is implemented by selecting a number $0<\kappa\leq 1$, which is
identified as the rate of reduction. The authors [8] define the tensor

$\mathbb{M}=\sum_{i=1}^{3}{\mathbf{e}}_{i}{\mathbf{e}}_{i}{\mathbf{e}}_{i}{\mathbf{e}}_{i}$
(24)

where ${\mathbf{e}}_{1}$, ${\mathbf{e}}_{2}$, ${\mathbf{e}}_{3}$ are the
orthonormal eigenvectors for $A$. The RSC replaces equations of the form

$\frac{DA}{Dt}=F(A)$ (25)

by

$\frac{DA}{Dt}=F(A)-(1-\kappa)\mathbb{M}:F(A)$ (26)

It turns out this form is simple to reproduce for the FEC. If equation (25) is
represented by the FEC

$\frac{DA}{Dt}=F(A,B),\qquad\frac{DB}{Dt}=G(A,B)$ (27)

then the effect of equation (26) is precisely modeled by the new FEC

$\frac{DA}{Dt}=F(A,B)-(1-\kappa)\mathbb{M}:F(A,B),\qquad\frac{DB}{Dt}=G(A,B)-(1-\kappa)\mathbb{M}:G(A,B)$
(28)

Finally, from a computational point of view, it should be noticed that if we
are working in the basis of orthonormal eigenvectors of $B$, then for any
symmetric matrix $N$ we have that $\mathbb{M}:N$ is simply the diagonal part
of $N$, that is, $\text{diag}(N_{11},N_{22},N_{33})$.

###  3.2 Is the solution to FEC always physical?

By the phrase “the solutions stay physical” we mean that $A$ stays positive
definite with trace one, that is, there exists a fiber orientation
distribution $\psi$ that satisfies equation (4). In fact, if $A$ ever ceases
to become positive definite, then not only is the Exact Closure going to give
the wrong answer, it even ceases to have a meaning in that equation (38) which
is used to define $A$ in terms of $B$ cannot be solved. Thus another way to
state “the solutions stay physical” is that $B$ stays positive definite and
finite, that is, none of the eigenvalues of $B$ become zero, and none of them
become infinite.

###### Theorem 1

The FEC solution to the isotropic diffusion equations (12) and (13) have
global in time physical solutions if $\Omega$, $\Gamma$ and $D_{r}$ are
bounded.

###### Theorem 2

The FEC solution to the anisotropic diffusion equations (14) and (15) have
global in time physical solutions if $D_{r}$ is positive definite, and
$\Omega$, $\Gamma$, $D(D_{r})/Dt$, $D_{r}$ and $1/\|D_{r}^{-1}\|$ are bounded.

where the proofs for both theorems are given in the Appendix beginning with
equation (65). Unfortunately Theorem 2 will not necessarily apply to the Koch
model [10] nor to the Phelps-Tucker ARD model [9], as there is no guarantee
that $1/\|D_{r}^{-1}\|$ is bounded nor, in the ARD case, that $D_{r}$ is
positive definite, unless extra hypotheses are applied.

###  3.3 Algorithm Summary

The algorithm to solve the FEC closure for the second-order orientation tensor
$A$ and the second-order tensor $B$ can be summarized as:

  1. 1.

Initialize $A$ and $B$, and define $\lambda$ along with any constants needed
for the diffusion model ${\mathcal{D}}\left[A\right]$

  2. 2.

At time $t_{i}$, rotate the tensors $A$ and $B$ into the principal frame of
$B$

  3. 3.

When the eigenvalues are distinct, use equation (16) for ${\mathbb{C}}$.
Otherwise when two eigenvalues are repeated, use equation (17) along with
equation (18), or in the case when three eigenvalues are repeated, use
equation (19).

  4. 4.

From ${\mathbb{C}}$, compute ${\mathbb{D}}$ using equation (20) in the
principal frame of $B$

  5. 5.

Compute $DA/Dt$ and $DB/Dt$ using either equations (12) and (13) for isotropic
diffusion or equations (14), (15) and (28) for the anisotropic diffusion
model, ARD-RSC. For the symmetric rank four tensor contractions with rank two
tensors, use equation (23) to reduce the number of redundant multiplication
operations.

  6. 6.

Rotate $DA/Dt$ and $DB/Dt$ into the flow reference frame, and extrapolate
$A\left(t_{i+1}\right)$ and $B\left(t_{i+1}\right)$ from time $t_{i}$ using
any standard ODE solver.

There are a number of coding issues we encountered, and we feel it will be
helpful to share as it will aid others in their computational implementations.

  * 1.

There is a choice to compute the basis of orthonormal eigenvectors from either
$A$ or $B$, where in theory these should be identical. We compute the basis
from $B$, arguing that the quantity $B$ is somehow more ‘fundamental’ and $A$
is ‘derived’ from $B$, which is true in the absence of diffusion.

  * 2.

We solve a ten dimensional set of ODEs, five for $A$, and five for $B$, where
one of the components of both $A$ and $B$ can be obtained, respectively, from
the relationships $\mathop{\text{tr}}A=1$ and $\det B=1$.

  * 3.

When computing $A$ from the orthonormal eigenvector basis of $B$, it is
important to force the off diagonal entries to be non-zero to limit numerical
drifting. In our studies, we found that failing to do this could cause an
adaptive ODE solver to completely freeze in select scenarios.

  * 4.

We set the ODE solver to work with a relative tolerance of $10^{-5}$, and
choose to use equations (18) or (19) when the eigenvalues were within
$10^{-4}$ of each other. This should cause $\mathbb{C}$ to be computed with an
accuracy of about $10^{-8}$ when using equations (16), and nearly machine
precision when using equations (18) or (19).

##  4 Numerical Results

Results are presented to demonstrate the accuracy improvements from employing
the FEC closure, and just as important to demonstrate the computational speed
advances over the similarly accurate orthotropic closures. In the present
examples, all flows have an initial isotropic orientation state designated by
$A_{11}=A_{22}=A_{33}=1/3$ and $B_{11}=B_{22}=B_{33}=1$, with all other
components of $A$ and $B$ being zero. The accuracy of the closure does not
depend on the initial orientation state, the isotropic orientation state is
chosen for uniformity. The equations of motion are solved using the FEC
closure for $A$ and $B$ from equations (12) and (13) for isotropic diffusion
or from equations (14), (15) and (28) for the anisotropic rotary diffusion
model with the reduced strain closure ARD-RSC from Phelps and Tucker [9]. For
comparison, the classical equations of motion for the second-order orientation
tensor $A$ requiring a curve-fitted closure for the fourth-order orientation
tensor ${\mathbb{A}}$, are solved using equations (6) and (7) for Folgar-
Tucker diffusion and equations (6), (8) and (25) for the ARD-RSC diffusion
model. Results are compared to solutions obtained using the Spherical Harmonic
approach [29] for solving the full distribution function equations (2) and
(3). It has been demonstrated in [29] that solutions using the Spherical
Harmonic approach are only limited in their accuracy by machine precision and
require considerably less computational effort than solutions using the
control volume approach of Bay [25]. Although a great reduction in speed and
an advancement in accuracy, the Spherical Harmonic approach still requires
more effort than the orientation tensor approach, nor does it readily lend
itself to an applicable form for coupling with commercial FEA solvers. We
select three commonly employed closures for comparisons. The first is the
classical Hybrid closure of Advani and Tucker [4] is selected as it is
regularly used in commercial and research codes due to its computational
efficiency and ease of implementation. The second is an orthotropic closure,
whose class of closures has found increasing use due to their considerable
accuracy improvements over the Hybrid closure. In our study we select the ORT
closure presented by VerWeyst and Tucker [1] based on the Wetzel closure [12].
Our third closure is that of the IBOF from Chung and Kwon [22] which is
claimed to be a more computationally efficient orthotropic closure as it uses
the invariants of $A$ as opposed to the eigenvalues of $A$ thus avoiding
costly tensor rotations.

###  4.1 Results: Simple Shear Flow

The first example is that of a pure shearing flow, given by $v_{1}=Gx_{3}$ and
$v_{2}=v_{3}=0$. Pure shearing flow is commonly employed (see e.g., [6, 22,
30]) to demonstrate a particular closure problem due to the oscillatory nature
of alignment inherent to the Jeffery fiber orbits. Two scenarios are
presented, the first of the Folgar-Tucker isotropic diffusion model in
equation (2) where $D_{r}=C_{I}\dot{\gamma}$, and the second scenario for the
ARD-RSC anisotropic diffusion model.

####  4.1.1 Simple Shear Flow Orientation

In industrial simulations, the Folgar-Tucker isotropic diffusion model
typically has interaction coefficients that range from $C_{I}=10^{-3}$ to
$C_{I}=10^{-2}$. The effective fiber aspect ratio ranges from 5 to 30
($a_{e}\simeq 1.4\times a_{r}$, where $a_{r}$ is the aspect ratio of
cylindrical fibers), which corresponds to a shape correction factor ranging
from $\lambda=0.96$ to $\lambda=0.999$. Two simulation results using isotropic
diffusion are presented in Figures 1(a) and (b), the first is for
$C_{I}=10^{-3}$ with $\lambda=0.99$ and the later for $C_{I}=10^{-2}$ with
$\lambda=0.95$. Results for the IBOF closure are not shown as they are nearly
graphically indistinguishable from the ORT closure results. It is important to
observe that the ORT and the FEC closure yield results that are graphically
indistinguishable and reasonably close to the orientation state predicted from
the numerically exact Spherical Harmonic solution. Conversely, the orientation
results from the Hybrid closure tend to over predict the the true orientation
state. It is important to point out the apparent oscillatory nature of the
transient solution for the Spherical Harmonic results when $C_{I}=10^{-3}$
with $\lambda=0.99$, which occurs to a lesser extent for $C_{I}=10^{-2}$.
These oscillations are expected due to the low amount of diffusion present.
Equally important is to notice that the oscillations from the FEC closure, as
well as the ORT, both damp out to the same steady state value. Note also that
the FEC does not oscillate excessively for either of the isotropic flow
conditions presented, which was a problem that plagued the early orthotropic
closures (see e.g., [6] and [31]) and the early neural network closures [32].
There remains room for further accuracy improvements (see e.g., [33] for
several preliminary higher accuracy closures). However, it is speculated based
upon the discussion in Jack and Smith [34] that such improvements will be
slight when solving the second-order moment equations, and higher order moment
simulations, such as those that use sixth-order closures (see e.g., [24]) may
need to be considered for significant accuracy improvements.

The Folgar-Tucker model has been used for decades, but tends to overstate the
rate of alignment during the transient solution (see e.g., [7]). The ARD-RSC
model [9] seeks to address these limitations, but few studies have focused on
this new diffusion model and the dependance of computed results on the choice
of closure. In the ARD-RSC model, the rotary diffusion coefficient of Folgar
and Tucker isotropic diffusion model ($D_{r}=C_{I}\dot{\gamma}$ where
$\dot{\gamma}=\left(\frac{1}{2}\Gamma:\Gamma\right)^{1/2}$) is replaced by an
anisotropic diffusion coefficient expressed by

$D_{r}=b_{1}\dot{\gamma}I+b_{2}\dot{\gamma}A+b_{3}\dot{\gamma}A^{2}+\tfrac{1}{2}b_{4}\Gamma+\tfrac{1}{4}b_{5}\dot{\gamma}^{-1}\Gamma^{2}$
(29)

where

$(b_{1},b_{2},b_{3},b_{4},b_{5})=(1.924\times 10^{-4},5.839\times
10^{-3},4.0\times 10^{-2},1.168\times 10^{-5},0)$ (30)

The ARD-RSC model serves as an excellent example of the effectiveness of the
FEC approach for solving the tensor form of orientation as the ARD-RSC model
will yield orientation states that are considerably different than that of the
Folgar-Tucker model. Results from the various closures and the spherical
harmonic results are presented in Figure 2 for the ARD-RSC flow with
$\kappa=1/30$. The value of $\kappa=1/30$ is taken from the results presented
in Phelps and Tucker [9], which was based on their experimental observations.
For a fiber aspect ratio of $\sim 5$, corresponding to $\lambda=0.95$, each of
the investigated closures produces graphically similar results. During the
initial flow stages, the Hybrid tends to over predict alignment, whereas the
ORT and the FEC tend to under predict alignment. As steady state is attained,
the FEC and the ORT yield nearly identical results, both of which over predict
$A_{11}$ in the final orientation state whereas the Hybrid yields a reasonable
representation of the orientation. For a long fiber, corresponding to
$\lambda\rightarrow 1$, the trends are similar to those of the lower aspect
ratio fibers, but in this case the FEC and the ORT better represent the final
orientation state relative to the Hybrid.

####  4.1.2 Orthotropic Closure Errors

The ORT is a polynomial approximation to the Exact Closure, as demonstrated in
the preceding section, and it is not surprising that the two approaches yield
graphically indistinguishable results for many of the flows investigated. On
closer inspection of the transient solution of the ARD-RSC model for
$\kappa=1/30$ and $\lambda=1$ there is a slight difference. This difference is
shown in Figure 3(a) where a closeup view is provided of the $A_{11}$
component for the flow times of $800\leq Gt\leq 1,200$. These results indicate
how well the fitting was performed in the construction of the ORT. As the ORT
is an approximation of the Exact Closure of Montgomery-Smith et al. [14] for
pure Jeffery’s flow, it is of interest to determine whether the slight
deviation comes from the Jeffery’s component or the diffusion component of
equation (6). To this end, we performed a comparison for the derivative of $A$
computed in two different ways. First, for each point in time $t$, we computed
$A(t)$ and $B(t)$ using the FEC method. Then we computed four quantities:
$\frac{DA^{\text{\tiny FEC, Diff}}}{Dt}$ which contains the terms from the
right hand side of equation (14) that explicitly include $D_{r}$,
$\frac{DA^{\text{\tiny FEC, Jeff}}}{Dt}$ which contains the terms from the
right hand side of equation (14) that do not involve $D_{r}$,
$\frac{DA^{\text{\tiny ORT, Diff}}}{Dt}$ the right hand side of equation (8),
and $\frac{DA^{\text{\tiny ORT, Jeff}}}{Dt}$ the right hand side of equation
(6) when $\mathcal{D}(A)$ is set to zero. In the latter two cases $\mathbb{A}$
is computed using the ORT closure. The error is then defined as

$\displaystyle E_{\text{\tiny
Diffusion}}=\sqrt{\sum_{i=1}^{3}\sum_{j=1}^{3}\left(\frac{DA_{ij}^{\text{\tiny
FEC, Diff}}}{Dt}-\frac{DA_{ij}^{\text{\tiny ORT, Diff}}}{Dt}\right)^{2}}$ (31)
$\displaystyle E_{\text{\tiny
Jeffery}}=\sqrt{\sum_{i=1}^{3}\sum_{j=1}^{3}\left(\frac{DA_{ij}^{\text{\tiny
FEC, Jeff}}}{Dt}-\frac{DA_{ij}^{\text{\tiny ORT, Jeff}}}{Dt}\right)^{2}}$ (32)

Each of the two errors are plotted in Figure 3(b). It is clear from the figure
that although the ORT’s derivative calculation from the diffusion component is
not zero, it is minor in comparison to the error from the Jeffery’s part of
the orientation tensor equation of motion. This error is only a rough
indication of the sources of error, but values of 0.04% at a given moment in
flow time can account for an error as large as 40% for $A$ for the flow times
on the order 1,000. Since the errors from each of the possible sources
probably do not drive the error in the solution toward the same direction, the
total error would be expected to be less than the upper bound of 40%, where in
reality the error is closer to 0.9% as steady state is approached.

Since the ORT and FEC differ by about 0.9%, it begs the question as to which
is more accurate in computing the true exact closure. While the FEC in theory
should exactly compute the exact closure, it is possible that numerical errors
creep into the FEC. To test for this, we performed a consistency check. After
finding the solution $A(t)$ and $B(t)$ using the FEC, we calculated

$E_{\text{\tiny
Exact}}=\sqrt{\sum_{i=1}^{2}\sum_{j=1}^{3}\left(A(B)_{ij}-A_{ij}\right)^{2}}$
(33)

where $A(B)$ was computed using equation (38). This calculation was performed
by diagonalizing $B$, applying the elliptic integrals in equation set (47)
using the software package [35], and then performing the reverse change of
basis. The results for the ARD-RSC model with $\kappa=1/30$ and $\lambda=1.00$
show an error of less than $10^{-8}$ throughout the transient solution, thus
suggesting the implementation as presented in this paper for the FEC is quite
accurate.

###  4.2 Results: Orientation Error Summary

To quantify the errors observed in Figures 1(a) and (b) for the isotropic
diffusion models, a series of fourteen flows are studied as outlined in table
1 where $\lambda=1$ for each of the flows. The solution is obtained using the
classical closure methods and the FEC closure results are compared to
solutions obtained from the Spherical Harmonic approach. To quantify the
error, the time average of the Frobenius Norm of the difference between the
true solution $A_{ij}^{\mbox{\tiny Spherical}}(t)$ and the approximate
solution obtained from a closure $A_{ij}^{\mbox{\tiny Closure}}(t)$ is
computed as

$\displaystyle\overline{E}_{\mbox{\tiny
Closure}}=\frac{1}{t_{f}-t_{0}}\int_{t_{0}}^{t_{f}}\sqrt{\sum_{i=1}^{3}\sum_{j=1}^{3}\left|A_{ij}^{\mbox{\tiny
Spherical}}(t)-A_{ij}^{\mbox{\tiny Closure}}(t)\right|^{2}}dt$ (34)

where $t_{0}$ is the initial time where the fiber orientation is isotropic and
$t_{f}$ is the time when the steady state is attained, which in this example
will be defined when the magnitude of the largest derivative of the
eigenvalues of $A$ is less than $G\times 10^{-4}$. This can be expressed as
the smallest moment in time when the following is satisfied
$\left(\max_{i\in\\{1,2,3\\}}|\frac{DA_{(i)}}{Dt}(t)|\right)\leq G\times
10^{-4}$. The quantitative error metric in equation (34) yields a value for
the simple shear flow of Figure 1(b) for the FEC, ORT and Hybrid closures of,
respectively, $4.74\times 10^{-2}$, $4.85\times 10^{-2}$ and $1.75\times
10^{-1}$. As the objective is to compare the relative accuracy improvements
between the FEC closure and the existing closures we will normalize the error
metric in equation (34) as

$\displaystyle\overline{\varepsilon}_{\mbox{\tiny Closure}}$
$\displaystyle\equiv$ $\displaystyle\frac{\overline{E}_{\mbox{\tiny
Closure}}}{\min\limits_{\mbox{\tiny Closure}}\left(\overline{E}_{\mbox{\tiny
Closure}}\right)}$ (35)

where the closure with the greatest accuracy will have a value of
$\overline{\varepsilon}_{\mbox{\tiny Closure}}=1$, and the remaining closures
will have a value of $\overline{\varepsilon}_{\mbox{\tiny Closure}}$ in excess
of 1. For each of the flows studied, the normalized error of equation (35) is
tabulated in Table 1 for the FEC, ORT, IBOF and the Hybrid closures. In each
of the flows considered, the FEC performs as well as or better than the
orthotropic closures.

###  4.3 Results: Combined Flow

A classical flow to demonstrate the effectiveness and robustness of a closure
is that of the combined flow presented in Cintra and Tucker [6]. This flow is
often selected as the orientation state crisscrosses the eigenspace of
possible orientations. The combined flow begins with pure shear in the
$x_{1}-x_{2}$ direction for $0\leq Gt<10$ defined by the velocity field
$v_{1}=Gx_{2}$, $v_{2}=v_{3}=0$. The flow then transitions to shearing flow in
the $x_{2}-x_{3}$ plane with stretching in the $x_{3}$ direction during the
time $10\leq Gt<20$ defined by the velocity field $v_{1}=-1/20Gx_{1}$,
$v_{2}=-1/20Gx_{2}+Gx_{3}$ and $v_{3}=1/10Gx_{3}$. The flow then transitions
to a flow with a considerable amount of stretching in the $x_{1}$ direction
with a reduced amount of shearing in the $x_{2}-x_{3}$ plane for $20\leq Gt$
defined by the velocity field $v_{1}=Gx_{1}$, $v_{2}=-1/2Gx_{2}+Gx_{1}$ and
$v_{3}=-1/2x_{3}$. The times where the flow transitions are chosen to prevent
the orientation from attaining steady state, thus any error in the transient
solution will be propagated to the next flow state. As observed in Figure 4
for flow results from the Folgar-Tucker model with $C_{I}=10^{-2}$ and
$\lambda=1$, the ORT and the FEC again yield similar results. This is
significant as it further demonstrates the robustness and the accuracy of the
FEC.

###  4.4 Results: Center-gated Disk Flow

The final flow investigated is that of the center-gated disk, a typical flow
condition in industrial processes [1, 36]. The flow enters the mold through
the pin gate and flows radially outward, where the velocity is a function of
both the gap height $2b$ and the radial distance from the gate $r$. The
velocity gradient for a Newtonian fluid can be represented by [6]

$\displaystyle v_{r}=\frac{3Q}{8\pi
rb}\left(1-\left(\frac{z}{b}\right)^{2}\right),\>\>\>\>v_{\theta}=v_{z}=0$
(36) $\displaystyle\frac{\partial v_{i}}{\partial x_{j}}=\frac{3Q}{8\pi
rb}\left[\begin{smallmatrix}-\frac{1}{r}\left(1-\frac{z^{2}}{b^{2}}\right)&0&-\frac{2}{b}\frac{z}{b}\\\
0&\frac{1}{r}\left(1-\frac{z^{2}}{b^{2}}\right)&0\\\ 0&0&0\\\
\end{smallmatrix}\right]$ (37)

where $z$ is the gap height location between the mold walls, $b$ is half the
gap height thickness, and $Q$ is the flow rate. Orientation results are
presented in Figure 5 for a gap height of $z/b=4/10$ for isotropic diffusion
with $C_{I}=10^{-2}$ and $\lambda=1$. Again, the Hybrid overshoots the actual
orientation state, whereas the ORT and the FEC behave in a graphically
identical fashion. This last result further demonstrates the robustness of the
FEC approach. Similar tests were performed for gap heights of
$z/b=0,1/10,2/10,\ldots,9/10$ and similar conclusions were observed at all gap
heights.

###  4.5 Results: Computational Time Enhancement

An additional goal for any new closure is that of reducing the computational
requirements for numerical solutions. Simulations are performed using in-house
developed single threaded code using Intel’s FORTRAN 90 compiler version 11.1.
Computations are solved on a standard desktop with an Intel i7 processor with
8 GB of Ram. The solution of the ORT has been studied by the investigators for
several years, and a reasonably efficient algorithm has been developed.
Solutions for the IBOF were made using the FORTRAN 90 code discussed in Jack
et al. [30].

Notice from Equations (12) and (13) that the operations
${\mathbb{C}}:\left[\cdots\right]$ and ${\mathbb{D}}:\left[\cdots\right]$ are
independent of coordinate frame. As we explained in equation (23), in the
principal frame there are a considerable number of terms in both
${\mathbb{C}}$ and ${\mathbb{D}}$ that are zero that are known prior to any
calculations, and thus operations involving $0$ can be avoided in the coding.
In addition, computing $DA/Dt$ and $DB/Dt$ in the principle reference frame
and then rotating the resulting $3\times 3$ tensors into the local reference
frame will be more efficient than rotating the $3\times 3\times 3\times 3$
tensors ${\mathbb{C}}$ and ${\mathbb{D}}$ into the local reference frame and
then computing $DA/Dt$ and $DB/Dt$. All computations of the FEC utilize this
characteristic, and thus greatly reduce the computational efforts. In
addition, redundant calculations from Equations (12) and (13) are closely
followed and performed only once. These computations are particularly frequent
in the double contractions of the fourth-order tensors with the second-order
tensors.

In the first study, computations were performed for the previous closure
operations for the ORT and the Hybrid using algorithms similar to
implementations discussed in the literature. In studies using an adaptive step
size solver, solutions for the IBOF took nearly 10 times that of the ORT,
whereas for the fixed step size the two closures required similar
computational efforts. To avoid any computational comparisons introduced by an
adaptive step size solver, computations were performed using a fixed step-size
fourth-order Runge-Kutta (R-K) solver with a very small step size of $\Delta
Gt=10^{-4}$. Computational times are tabulated in Table 2 for both CPU time
and normalized time. Normalized time is defined based off of the often
employed Hybrid closure using the standard implementation for the Hybrid
closure with the very small step size. The ORT required nearly 770 seconds, a
factor of 31 times greater than that of the Hybrid. Conversely, the FEC
required only 26 seconds, a slight increase in effort beyond the Hybrid, which
required 25 seconds. This is very striking as the Hybrid closure is often
selected in research and industrial codes due to its computational efficiency,
while recognizing the sacrifice in computational accuracy. This is no longer
the case with the FEC as it has the same accuracy of the orthotropic closures
while providing computational speeds nearly identical to that of the Hybrid
closure.

In the process of developing the FEC algorithm, it was observed that many
redundant operations existed in the implementation of the ORT and the Hybrid
closures. For existing implementation of the classical closures, no special
consideration was given to the ${\mathbb{A}}:\Gamma$ term, but since the rank
four tensor ${\mathbb{A}}$ is symmetric, equation (23) can be used to reduce
the number operations of the double contraction to that of a simple rank two
tensor operations for both the hybrid closure and the ORT closure
implementations. For the ORT, the computational problem can be further
simplified by constructing the second-order tensor $DA/Dt$ in the principal
frame, and then performing the tensor rotation back into the reference frame.
Thus the costly rotations of the fourth-order tensor ${\mathbb{A}}$ are
avoided. These optimized results for the Hybrid and the ORT are shown in Table
2, and it is clear that the computational times were greatly reduced. The
optimized Hybrid implementation reduced the computational time to 30% of the
original time, whereas the ORT implementation improved by over an order of
magnitude. With these additional computational advances the ORT appears to be
a more viable alternative to the Hybrid, but the FEC still has similar
computational requirements. It is expected that with further studies, the FEC
algorithm could be improved to further reduce its computational times.

##  5 Conclusion

The Fast Exact Closure is a robust, computationally efficient, approach to
solve the fiber orientation equations of motion for the orientation tensors.
This unique approach does not require any form of curve fitting based on
orientation data obtained from numerical solutions of the full fiber
orientation distribution. The results presented demonstrate that the FEC is as
accurate and robust as the existing industrially accepted closures, while
enjoying computational speeds equivalent to the industrial form of the hybrid
closure.

##  6 Acknowledgments

The authors gratefully acknowledge support from the N.S.F. via grant C.M.M.I.
0727399 and Baylor University through their financial support through their
faculty member start-up package.

## Appendix: Justification and Proofs

By [14], the Exact Closure is this. Given $A$, compute the symmetric matrix
$B$ by solving

$A=A(B)=\tfrac{1}{2}\int_{0}^{\infty}\frac{(B+sI)^{-1}\,ds}{\sqrt{\text{det}(B+sI)}}$
(38)

It was shown in [14] that $B$ is unique with this property. Then compute
$\mathbb{A}$ using the formula

$\mathbb{A}=\tfrac{3}{4}\int_{0}^{\infty}\frac{s\,\mathcal{S}((B+sI)^{-1}\otimes(B+sI)^{-1})\,ds}{\sqrt{\text{det}(B+sI)}}$
(39)

Here $\mathcal{S}$ represents the symmetrization of a rank 4 tensor, that is,
$\mathcal{S}(\mathbb{B})_{ijkl}$ is the average of $\mathbb{B}_{mnpq}$ over
all permutations $(m,n,p,q)$ of $(i,j,k,l)$.

It can be shown that the following two statements are equivalent:

  1. 1.

Equation (38) holds for all time.

  2. 2.

Equation (38) holds at $t=0$, and equation (9) holds for all time, where

$\mathbb{C}=\tfrac{3}{4}\int_{0}^{\infty}\frac{\mathcal{S}((B+sI)^{-1}\otimes(B+sI)^{-1})\,ds}{\sqrt{\text{det}(B+sI)}}$
(40)

Furthermore, it can be shown for every symmetric matrix $M$ that

$\mathop{\text{tr}}(B^{-1}\cdot M)=2\mathop{\text{tr}}(\mathbb{C}:M)$ (41)

and hence it can be seen that $\mathop{\text{tr}}(DA/Dt)=0$ if and only if
$\mathop{\text{tr}}(B^{-1}\cdot(DB/Dt))=0$, that is, $\mathop{\text{tr}}A$
stays constant if and only if $\det B$ stays constant.

Next, we have

The linear map on symmetric matrices $M\mapsto\mathbb{C}:M$ is invertible (42)

that is, there exists a rank 4 tensor $\mathbb{D}$ such that

$\mathbb{C}:\mathbb{D}:M=\mathbb{D}:\mathbb{C}:M=M\text{ for any symmetric
matrix $M$}$ (43)

Indeed if we define the six by six matrix

$\mathcal{C}=\left[\begin{smallmatrix}\mathbb{C}_{1111}&\mathbb{C}_{1122}&\mathbb{C}_{1133}&2\mathbb{C}_{1112}&2\mathbb{C}_{1113}&2\mathbb{C}_{1123}\\\
\mathbb{C}_{2211}&\mathbb{C}_{2222}&\mathbb{C}_{2233}&2\mathbb{C}_{2212}&2\mathbb{C}_{2213}&2\mathbb{C}_{2223}\\\
\mathbb{C}_{3311}&\mathbb{C}_{3322}&\mathbb{C}_{3333}&2\mathbb{C}_{3312}&2\mathbb{C}_{3313}&2\mathbb{C}_{3323}\\\
2\mathbb{C}_{1211}&2\mathbb{C}_{1222}&2\mathbb{C}_{1233}&4\mathbb{C}_{1212}&4\mathbb{C}_{1213}&4\mathbb{C}_{1223}\\\
2\mathbb{C}_{1311}&2\mathbb{C}_{1322}&2\mathbb{C}_{1333}&4\mathbb{C}_{1312}&4\mathbb{C}_{1313}&4\mathbb{C}_{1323}\\\
2\mathbb{C}_{2311}&2\mathbb{C}_{2322}&2\mathbb{C}_{2333}&4\mathbb{C}_{2312}&4\mathbb{C}_{2313}&4\mathbb{C}_{2323}\end{smallmatrix}\right]$
(44)

then $\mathbb{D}$ can be calculated using the formula

$\displaystyle\left[\begin{smallmatrix}\mathbb{D}_{1111}&\mathbb{D}_{1122}&\mathbb{D}_{1133}&\mathbb{D}_{1112}&\mathbb{D}_{1113}&\mathbb{D}_{1123}\\\
\mathbb{D}_{2211}&\mathbb{D}_{2222}&\mathbb{D}_{2233}&\mathbb{D}_{2212}&\mathbb{D}_{2213}&\mathbb{D}_{2223}\\\
\mathbb{D}_{3311}&\mathbb{D}_{3322}&\mathbb{D}_{3333}&\mathbb{D}_{3312}&\mathbb{D}_{3313}&\mathbb{D}_{3323}\\\
\mathbb{D}_{1211}&\mathbb{D}_{1222}&\mathbb{D}_{1233}&\mathbb{D}_{1212}&\mathbb{D}_{1213}&\mathbb{D}_{1223}\\\
\mathbb{D}_{1311}&\mathbb{D}_{1322}&\mathbb{D}_{1333}&\mathbb{D}_{1312}&\mathbb{D}_{1313}&\mathbb{D}_{1323}\\\
\mathbb{D}_{2311}&\mathbb{D}_{2322}&\mathbb{D}_{2333}&\mathbb{D}_{2312}&\mathbb{D}_{2313}&\mathbb{D}_{2323}\end{smallmatrix}\right]=\mathcal{C}^{-1}$
(45) $\displaystyle\mathbb{D}_{ijkl}=\mathbb{D}_{jikl}=\mathbb{D}_{ijlk}$ (46)

In the basis of orthonormal eigenvectors of $B$, since $\mathbb{C}_{ijkk}=0$
whenever $i\neq j\neq k$, this reduces to equation (20).

Next, if $B$ is diagonal, then $A$ is diagonal with entries

$\begin{split}a_{1}=\tfrac{1}{2}\int_{0}^{\infty}\frac{ds}{(b_{1}+s)^{3/2}\sqrt{b_{2}+s}\sqrt{b_{3}+s}}\\\
a_{2}=\tfrac{1}{2}\int_{0}^{\infty}\frac{ds}{\sqrt{b_{1}+s}(b_{2}+s)^{3/2}\sqrt{b_{3}+s}}\\\
a_{3}=\tfrac{1}{2}\int_{0}^{\infty}\frac{ds}{\sqrt{b_{1}+s}\sqrt{b_{2}+s}(b_{3}+t)^{3/2}}\end{split}$
(47)

and

$\begin{split}\mathbb{C}_{1111}=\tfrac{3}{4}\int_{0}^{\infty}\frac{ds}{(b_{1}+s)^{5/2}\sqrt{b_{2}+s}\sqrt{b_{3}+s}}&\qquad\quad\mathbb{C}_{1122}=\tfrac{1}{4}\int_{0}^{\infty}\frac{ds}{(b_{1}+s)^{3/2}(b_{2}+s)^{3/2}\sqrt{b_{3}+s}}\\\
\mathbb{C}_{2222}=\tfrac{3}{4}\int_{0}^{\infty}\frac{ds}{\sqrt{b_{1}+s}(b_{2}+s)^{5/2}\sqrt{b_{3}+s}}&\qquad\quad\mathbb{C}_{2233}=\tfrac{1}{4}\int_{0}^{\infty}\frac{ds}{\sqrt{b_{1}+s}(b_{2}+s)^{3/2}(b_{3}+s)^{3/2}}\\\
\mathbb{C}_{3333}=\tfrac{3}{4}\int_{0}^{\infty}\frac{ds}{\sqrt{b_{1}+s}\sqrt{b_{2}+s}(b_{3}+t)^{5/2}}&\qquad\quad\mathbb{C}_{1133}=\tfrac{1}{4}\int_{0}^{\infty}\frac{ds}{(b_{1}+s)^{3/2}\sqrt{b_{2}+s}(b_{3}+s)^{3/2}}\end{split}$
(48)

all the other coefficients of $\mathbb{C}$ being zero. Furthermore

$\mathbb{C}:I=\tfrac{1}{2}B^{-1}$ (49)

Equations (16) now follow by an easy calculation. Equations (18) and (19) are
obtained by expanding the fourth equation of (48) using Taylor’s series, where

$\mathcal{I}_{n}=\int_{0}^{\infty}\frac{ds}{(b_{0}+s)^{n}\sqrt{b_{3}+s}}$ (50)

The proofs of various details now follow.

Proof of equations (9) and (40): Write $\dot{A}$ and $\dot{B}$ for
$\frac{DA}{Dt}$ and $\frac{DB}{Dt}$ respectively. Use the formulae

$\frac{D}{Dt}B^{-1}=-B^{-1}\cdot\dot{B}\cdot
B^{-1}\quad\text{and}\quad\frac{D}{Dt}\det
B=\mathop{\text{tr}}(B^{-1}\cdot\dot{B})\det B$ (51)

to obtain

$\begin{split}\dot{A}&=-\tfrac{1}{2}\int_{0}^{\infty}\frac{(B+sI)^{-1}\cdot\dot{B}\cdot(B+sI)^{-1}\,ds}{\sqrt{\text{det}(B+sI)}}-\tfrac{1}{4}\int_{0}^{\infty}\frac{[(B+sI)^{-1}:\dot{B}]\,(B+sI)^{-1}\,ds}{\sqrt{\text{det}(B+sI)}}\\\
&=-\mathbb{C}:\dot{B}\end{split}$ (52)

since for any symmetric matrix $K$ we have

$\mathcal{S}(K\otimes K):\dot{B}=\tfrac{1}{3}K\cdot\dot{B}\cdot
K+\tfrac{2}{3}(K:\dot{B})K$ (53)

Proof of equation (10): For any invertible symmetric matrix $K$

$\mathcal{S}(K^{-1}K^{-1}):(K\cdot M)=\mathcal{S}(K^{-1}K^{-1}):(M^{T}\cdot
K)=\tfrac{1}{3}((\mathop{\text{tr}}M)K^{-1}+M\cdot K^{-1}+K^{-1}\cdot M^{T})$
(54)

Setting $K=B+sI$, we multiply both sides by $3/(4\sqrt{\det(B+sI)})$, and
integrate with respect to $s$ from zero to infinity, to obtain equation (10).

Proof of equations (17) from (50): To compute $\mathcal{I}_{1}$, use the
formulae

$\displaystyle\frac{d}{ds}\cos^{-1}\left(\sqrt{\frac{b_{3}+s}{b_{0}+s}}\right)=-\frac{\sqrt{b_{0}-b_{3}}}{2(b_{0}+s)\sqrt{b_{3}+s}}$
(55)
$\displaystyle\frac{d}{ds}\cosh^{-1}\left(\sqrt{\frac{b_{3}+s}{b_{0}+s}}\right)=-\frac{\sqrt{b_{3}-b_{0}}}{2(b_{0}+s)\sqrt{b_{3}+s}}$
(56)

Next, integrating by parts, we obtain

$\mathcal{I}_{n}=-\frac{\sqrt{b_{3}}}{2b_{0}^{n}}+\frac{n}{2}\int_{0}^{\infty}\frac{\sqrt{b_{3}+s}\,ds}{(b_{0}+s)^{n+1}}$
(57)

and simple algebra gives

$\int_{0}^{\infty}\frac{\sqrt{b_{3}+s}\,ds}{(b_{0}+s)^{n+1}}=\mathcal{I}_{n}+(b_{3}-b_{0})\mathcal{I}_{n+1}$
(58)

Proof of equation (41): For any positive definite matrix $X$, if

$A(X)=\tfrac{1}{2}\int_{0}^{\infty}\frac{(X+sI)^{-1}\,ds}{\sqrt{\text{det}(X+sI)}}$
(59)

then

$\mathop{\text{tr}}(A(X))=\tfrac{1}{2}\int_{0}^{\infty}\frac{\mathop{\text{tr}}((X+sI)^{-1})\,ds}{\sqrt{\text{det}(X+sI)}}=-\int_{0}^{\infty}\frac{d}{ds}\left(\frac{1}{\sqrt{\text{det}(X+sI)}}\right)\,ds=\frac{1}{\sqrt{\det
X}}$ (60)

If $\mathop{\text{tr}}(B^{-1}\cdot M)=\alpha$, then (remembering that $\det
B=1$) we have $\det(B+\epsilon M)=1+\epsilon\alpha+O(\epsilon^{2})$ as
$\epsilon\to 0$. Hence
$1-\tfrac{1}{2}\epsilon\alpha+O(\epsilon^{2})=\mathop{\text{tr}}(A(B+\epsilon
M))=1-\epsilon\mathop{\text{tr}}(\mathbb{C}:M)+O(\epsilon^{2})$. Therefore
$\mathop{\text{tr}}(\mathbb{C}:M)=\tfrac{1}{2}\alpha$.

Proof of equation (42): This follows because

$\text{$M$ is a symmetric non-zero matrix}\Rightarrow M:\mathbb{C}:M>0$ (61)

and hence $M\neq 0\Rightarrow\mathbb{C}:M\neq 0$.

To see this, suppose that $K$ is a positive definite three by three matrix,
and let $k_{1}$, $k_{2}$ and $k_{3}$ be its eigenvalues. Then in the basis of
corresponding orthonormal eigenvalues of $K$, we have that for any non-zero
symmetric $M$

$M:\mathcal{S}(K\otimes
K):M=\tfrac{1}{3}\left(\sum_{i=1}^{3}k_{i}M_{ii}\right)^{2}+\tfrac{2}{3}\sum_{i,j=1}^{3}k_{i}k_{j}M_{ij}^{2}>0$
(62)

Apply this to $K=(B+sI)^{-1}$, multiply by $(\det(B+sI))^{-1/2}$, and then
integrate over $s$ to obtain $M:\mathbb{C}:M>0$.

Proof of equation (49): Without loss of generality $B$ is diagonal. Hence we
need to prove statements such as
$\mathbb{C}_{1111}+\mathbb{C}_{1122}+\mathbb{C}_{1133}=\tfrac{1}{2}b_{1}^{-1}$
when $\mathbb{C}$ satisfies equation (48). But

$\begin{split}\mathbb{C}_{1111}+\mathbb{C}_{1122}+\mathbb{C}_{1133}&=-\tfrac{1}{2}\int_{0}^{\infty}\frac{d}{ds}\left(\frac{1}{(b_{1}+s)^{3/2}\sqrt{b_{2}+s}\sqrt{b_{3}+s}}\right)\,ds\\\
&=\tfrac{1}{2}b_{1}^{-3/2}b_{2}^{-1/2}b_{3}^{-1/2}\end{split}$ (63)

The result follows since $b_{1}b_{2}b_{3}=1$.

Proof of equation (28): From equation (9), we see that the RSC version of
equation (27) is equation (28) and

$\frac{DB}{Dt}=G(A,B)-(1-\kappa)\mathbb{D}:\mathbb{M}:F(A,B)$ (64)

Since $\mathbb{D}:F(A,B)=G(A,B)$, it follows that all we need to show is
$\mathbb{D}:\mathbb{M}:F(A,B)=\mathbb{M}:\mathbb{D}:F(A,B)$. This is easily
seen by working in the basis of orthonormal eigenvectors of $B$, noticing that
then $\mathbb{M}:N$ is simply the diagonal part of $N$, and applying equation
(23).

Proof of Theorem 1: It follows from $\det B=1$ that the only way that the
solutions can become non-physical is if $B$ ‘blows up,’ that is, if one or
more of the eigenvalues of $B$ become infinite in finite time. (Also, [37,
Theorem 1.4] can be used to show that the finiteness of the eigenvalues of $B$
imply the differential equations have a unique solution.)

Substituting $M=I$ into equation (10), we obtain $\mathbb{C}:B=\tfrac{3}{2}A$,
that is,

$\mathbb{D}:A=\tfrac{2}{3}B$ (65)

Take the trace of equation (13) and use equation (65), to obtain

$\frac{D}{Dt}\mathop{\text{tr}}B\leq
c(\|\Omega\|+\|\Gamma\|+D_{r})(\mathop{\text{tr}}B)-2D_{r}(I:\mathbb{D}:I)$
(66)

for some universal constant $c>0$. Here $\|\cdot\|$ denotes the spectral norm
of a matrix, and we have used the inequality $\mathop{\text{tr}}(X\cdot
Y)\leq\|X\|\mathop{\text{tr}}Y$ whenever $Y$ is positive definite. By equation
(61), we have $I:\mathbb{D}:I=M:\mathbb{C}:M\geq 0$, where $M=\mathbb{D}:I$,
and hence

$\frac{D}{Dt}\mathop{\text{tr}}B\leq
c(\|\Omega\|+\|\Gamma\|+D_{r})(\mathop{\text{tr}}B)$ (67)

Now we can apply Gronwall’s inequality [37, Chapter 2.1.1] (in Lagrangian
coordinates) to obtain

$\mathop{\text{tr}}B\leq(\mathop{\text{tr}}B_{0})e^{ctL}$ (68)

where $L$ is an upper bound for $\|\Omega\|+\|\Gamma\|+D_{r}$, and $B_{0}$ is
the value of $B$ at $t=0$. Therefore $\mathop{\text{tr}}B$ remains finite, and
since $B$ is positive definite, no eigenvalue of $B$ blows up to infinity in
finite time.

Proof of Theorem 2: Note that the positive definiteness of $D_{r}$, and the
boundedness of $D_{r}$ and $1/\|D_{r}^{-1}\|$ guarantee that the ratio of
$\mathop{\text{tr}}B$ and $D_{r}:B$ is bounded from above and below. From
equation (15), and using equation (65)

$\frac{D}{Dt}(D_{r}:B)\leq{}c\left(\|\Omega\|+\|\Gamma\|+\|D_{r}\|+\left\|\frac{D(D_{r})}{Dt}\right\|\right)(D_{r}:B)-2(D_{r}:\mathbb{D}:D_{r})$
(69)

The rest of the proof proceeds by a similar argument as above.

##  7 References

## References

  * [1] VerWeyst, B.E., C. Tucker, P. Foss, J. O’Gara, Fiber Orientation in 3-D Injection Molded Features: Prediction and Experiment, International Polymer Processing 14 (1999) 409–420.
  * [2] Fan, X., N. Phan-Thien, R. Zheng, A Direct Simulation of Fibre Suspensions, Jn. of Non-Newtonian Fluid Mechanics 74 (1998) 113–135.
  * [3] Folgar, F.P., C. Tucker, Orientation Behavior of Fibers in Concentrated Suspensions, Jn. of Reinforced Plastics and Composites 3 (1984) 98–119.
  * [4] Advani, S.G., C. Tucker, The Use of Tensors to Describe and Predict Fiber Orientation in Short Fiber Composites, Jn. of Rheology 31 (8) (1987) 751–784.
  * [5] Jack, D.A., D. Smith, The Effect of Fiber Orientation Closure Approximations on Mechanical Property Predictions, Composites, Part A 38 (2007) 975–982.
  * [6] Cintra, J. S., C. Tucker, Orthotropic Closure Approximations for Flow-Induced Fiber Orientation, Jn. of Rheology 39 (6) (1995) 1095–1122.
  * [7] Tucker, C.L., J. Wang, J. O’Gara, G. DeBarr, Improved Fiber Orientation Predictions for Injection Molded Composites, in: NSF/DOE/APC Workshop: The Future of Modeling in Composites Molding Processes, Washington, D.C., 2004.
  * [8] Wang, J., J. O’Gara, C. Tucker, An Objective Model for Slow Orientation Kinetics in Concentrated Fiber Suspensions: Theory and Rheological Evidence, Journal of Rheology 52 (5) (2008) 1179–1200.
  * [9] Phelps, J.H., C. Tucker, An Anisotropic Rotary Diffusion Model for Fiber Orienation in Short- and LongFiber Thermoplastics, Journal of Non-Newtonian Fluid Mechanics 156 (2009) 165–176.
  * [10] Koch, D.L., A Model for Orientational Diffusion in Fiber Suspensions, Physics of Fluids 7 (8) (1995) 2086–2088.
  * [11] Jack, D.A., S. Montgomery-Smith, D. Smith, Anisotropic Diffusion Model for Suspensions of Short-Fiber Composite Processes., in: The XVth International Congress on Rheology, the Society of Rheology 80th Annual Meeting, The Society of Rheology, Monterey, CA, 2008.
  * [12] Wetzel, E.D., Modeling Flow-Induced Microstructure of Inhomogeneous Liquid-Liquid Mixtures, Ph.D. thesis, University of Illinois at Urbana Champaign (1999).
  * [13] VerWeyst, B.E., Numerical Predictions of Flow Induced Fiber Orientation in Three-Dimensional Geometries, Ph.D. thesis, University of Illinois at Urbana Champaign (1998).
  * [14] Montgomery-Smith, S.J. W. He, D.A. Jack, D.E. Smith, Exact Tensor Closures for the Three Dimensional Jeffery’s Equation, Under Review Journal of Fluid Mechanics, draft available at http://www.math.missouri.edu/~stephen/preprints/exact-closure.html (2010).
  * [15] Jeffery, G.B., The Motion of Ellipsoidal Particles Immersed in a Viscous Fluid, Proceedings of the Royal Society of London A 102 (1923) 161–179.
  * [16] Bird, R. B., C. Curtiss, R. C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, 2nd Edition, Vol. 2: Kinetic Theory, John Wiley & Sons, Inc., New York, NY, 1987.
  * [17] Petrie, C.J.S., The Rheology of Fibre Suspensions, Journal of Non-Newtonian Fluid Mechanics 87 (1999) 369–402.
  * [18] Hinch, E.J., L. Leal, Time-Dependent Shear Flows of a Suspension of Particles with Weak Brownian Rotations, Journal of Fluid Mechanics 57 (1973) 753–767.
  * [19] Altan, M.C., S. Advani, S. Güçeri, R. Pipes, On the Description of the Orientation State for Fiber Suspensions in Homogeneous Flows, Jn. of Rheology 33 (7) (1989) 1129–1155.
  * [20] Altan, M.C., S. Subbiah, S. Guceri, R. Pipes, Numerical Prediction of Three-Dimensional Fiber Orientation in Hele-Shaw Flows, Polymer Engineering and Science 30 (14) (1990) 848–859.
  * [21] Verleye, V., F. Dupret, Prediction of Fiber Orientation in Complex Injection Molded Parts, in: Developments in Non-Newtonian Flows, 1993, pp. 139–163.
  * [22] Chung, D.H., T. Kwon, Invariant-Based Optimal Fitting Closure Approximation for the Numerical Prediction of Flow-Induced Fiber Orientation, Jn. of Rheology 46 (1) (2002) 169–194.
  * [23] Han, K.-H., Y.-T. Im, Numerical Simulation of Three-Dimensional Fiber Orientation in Short-Fiber-Reinforced Injection-Molded Parts, Jn. of Materials Processing Technology 124 (2002) 366–371.
  * [24] Jack, D.A., D. Smith, An Invariant Based Fitted Closure of the Sixth-order Orientation Tensor for Modeling Short-Fiber Suspensions, Jn. of Rheology 49 (5) (2005) 1091–1116.
  * [25] Bay, R.S., Fiber Orientation in Injection Molded Composites: A Comparison of Theory and Experiment, Ph.D. thesis, University of Illinois at Urbana-Champaign (August 1991).
  * [26] Dinh, S.M., R. Armstrong, A Rheological Equation of State for Semiconcentrated Fiber Suspensions, Jn. of Rheology 28 (3) (1984) 207–227.
  * [27] Lipscomb, G G II., M. Denn, D. Hur, D. Boger, Flow of Fiber Suspensions in Complex Geometries, Jn. of Non-Newtonian Fluid Mechanics 26 (1988) 297–325.
  * [28] Altan, M.C., L. Tang, Orientation tensors in simple flows of dilute suspensions of non-brownian rigid ellipsoids, comparison of analytical and approximate solutions, Rheologica Acta 32 (1993) 227–244.
  * [29] Montgomery-Smith, S.J., D. Jack, D. Smith, A Systematic Approach to Obtaining Numerical Solutions of Jeffery’s Type Equations using Spherical Harmonics, Composites Part A 41 (2010) 827–835.
  * [30] D. Jack, B. Schache, D. Smith, Neural Network Based Closure for Modeling Short-Fiber Suspensions, Polymer CompositesAccepted for Publication.
  * [31] Chung, D.H., T. Kwon, Improved Model of Orthotropic Closure Approximation for Flow Induced Fiber Orientation, Polymer Composites 22 (5) (2001) 636–649.
  * [32] Qadir, N., D. Jack, Modeling Fibre Orientation in Short Fibre Suspensions Using the Neural Network-Based Orthotropic Closure, Composites, Part A.
  * [33] Mullens, M., Developing New Fitted Closure Approximations for Short-Fiber Reinforced Polymer Composites, Master’s thesis, University of Missouri - Columbia (July 2010).
  * [34] Jack, D.A., D. Smith, Assessing the Use of Tensor Closure Methods With Orientation Distribution Reconstruction Functions, Jn. of Composite Materials 38 (21) (2004) 1851–1872.
  * [35] GNU Scientific Library, http://www.gnu.org/software/gsl.
  * [36] Jack D.A., D. Smith, Elastic Properties of Short-Fiber Polymer Composites, Derivation and Demonstration of Analytical Forms for Expectation and Variance from Orientation Tensors, Journal of Composite Materials 42 (3) (2008) 277–308.
  * [37] Chicone, C., Ordinary Differential Equations with Applications, 2nd Edition, Springer-Verlag, New York, 2006.

Flow # | $v_{1}$ | $v_{2}$ | $v_{3}$ | $C_{I}$ | $\lambda$ | $\overline{\varepsilon}_{\mbox{\tiny FEC}}$ | $\overline{\varepsilon}_{\mbox{\tiny ORT}}$ | $\overline{\varepsilon}_{\mbox{\tiny IBOF}}$ | $\overline{\varepsilon}_{\mbox{\tiny Hybrid}}$
---|---|---|---|---|---|---|---|---|---
1 | $Gx_{1}$ | $Gx_{2}$ | $-2Gx_{3}$ | $10^{-3}$ | 1 | 1 | 2.06 | 1.28 | 76.2
2 | $2Gx_{1}$ | $-Gx_{2}$ | $-Gx_{3}$ | $10^{-3}$ | 1 | 1.03 | 1.05 | 1 | 25.8
3a | $Gx_{3}$ | $0$ | $0$ | $10^{-3}$ | 0.99 | 1 | 1.02 | 1.02 | 3.69
3b | $Gx_{3}$ | $0$ | $0$ | $10^{-3}$ | 1 | 1 | 1.02 | 1.01 | 3.28
4 | $-Gx_{1}+10Gx_{2}$ | $-Gx_{2}$ | $2Gx_{3}$ | $10^{-3}$ | 1 | 1 | 2.25 | 1.36 | 12.9
5 | $-Gx_{1}+Gx_{2}$ | $-Gx_{2}$ | $2Gx_{3}$ | $10^{-3}$ | 1 | 1.02 | 1 | 1.23 | 22.6
6 | $Gx_{1}+2Gx_{3}$ | $Gx_{2}$ | $-2Gx_{3}$ | $10^{-2}$ | 1 | 1 | 1.01 | 1.08 | 3.57
7 | $Gx_{1}+2.75Gx_{3}$ | $Gx_{2}$ | $-2Gx_{3}$ | $10^{-2}$ | 1 | 1 | 1.02 | 1.05 | 2.98
8 | $Gx_{1}+1.25Gx_{3}$ | $Gx_{2}$ | $-2Gx_{3}$ | $10^{-2}$ | 1 | 1.02 | 1 | 1.12 | 3.85
9 | $-Gx_{1}+10Gx_{3}$ | $Gx_{2}$ | $0$ | $10^{-2}$ | 1 | 1.03 | 1 | 1.03 | 1.65
10 | $-Gx_{1}+Gx_{3}$ | $Gx_{2}$ | $0$ | $10^{-2}$ | 1 | 1.01 | 1 | 1.04 | 2.29
11 | $2Gx_{1}+3Gx_{3}$ | $-Gx_{2}$ | $-Gx_{3}$ | $10^{-2}$ | 1 | 1.04 | 1.04 | 1 | 2.52
12 | $-Gx_{1}+3.75Gx_{2}$ | $Gx_{2}$ | $2Gx_{3}$ | $10^{-2}$ | 1 | 1 | 1.03 | 1.06 | 2.03
13 | $-Gx_{1}+1.5Gx_{2}$ | $-Gx_{2}$ | $2Gx_{3}$ | $10^{-2}$ | 1 | 1.00 | 1 | 1.01 | 2.34
14a | $Gx_{3}$ | $0$ | $0$ | $10^{-2}$ | 0.99 | 1 | 1.00 | 1.03 | 4.14
14b | $Gx_{3}$ | $0$ | $0$ | $10^{-2}$ | 1 | 1 | 1.00 | 1.02 | 3.90
Table 1: Flows used, and the resulting error computation in computing the second-order orientation tensor $A$. Closure | CPU Time | Normalized Time
---|---|---
Hybrid - Original | 25 | 1
Hybrid - Optimized | 6.9 | 0.3
ORT - Original | 770 | 31
ORT - Optimized | 21 | 0.8
FEC | 26 | 1.0
Table 2: Normalized Computational Times

Figure 1: Transient Solution for selected components of $A$ for simple shear
flow under isotropic diffusion (a) $C_{I}=10^{-3}$ and $\lambda=0.99$ and (b)
$C_{I}=10^{-2}$ and $\lambda=0.95$.

Figure 2: Transient Solution for selected components of $A$ for simple shear
flow under anisotropic rotary diffusion (ARD-RSC) (a) $\kappa=1/30$ and
$\lambda=0.95$ and (b) $\kappa=1/30$ and $\lambda=1.0$.

Figure 3: Anisotropic rotary diffusion results, simple shear $\kappa=1/30$ and
$\lambda=1.0$ (a) Selected time range of for $A_{11}$ (b) Transient error in
derivative computation for the fitted orthotropic closure ORT compared to FEC.

Figure 4: Transient solution for selected components of $A$ for mixed flow
from the Folgar-Tucker model with $C_{I}=10^{-2}$ and $\lambda=1.0$.

Figure 5: Transient solution for selected components of $A$ for center-gated
disk flow from the Folgar-Tucker model with $C_{I}=10^{-2}$ and $\lambda=1.0$
for $r/b=4/10$.