The Fast Exact Closure for Jeffery’s Equation with Diffusion

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.

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

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

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

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

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

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

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

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

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]).

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

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

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

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

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:

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

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

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

Then replace the first equation of equation (16) by

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

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

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

then decompose

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

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}$.

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.

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.

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.

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

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

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

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

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

Indeed if we define the six by six matrix

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

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

and