NUMERICAL RESEARCH ON THE THREE-DIMENSIONAL FIBER ORIENTATION DISTRIBUTION IN PLANAR SUSPENSION FLOWS

Abstract

To describe flow-induced fiber orientation, the Fokker-Planck equation is widely applied in the processing of composites and fiber suspensions. The analytical solution only exists when the Péclet number is infinite. So developing a numerical method covering a full range of Péclet number is of great significance. To accurately solve the Fokker-Planck equation, a numerical scheme based on the finite volume method is developed. Using spherical symmetry, the boundary is discretized and formulated into a cyclic tridiagonal matrix which is further solved by the CTDMA algorithm. To examine its validity, benchmark tests over a wide range of Péclet number are performed in a simple shear flow. For Pe=∞, the results agree well with the analytical solutions. For the other Pe numbers, the results are compared to results available in the literature. The tests show that this algorithm is accurate, stable, and globally conservative. Furthermore, this algorithm can be extended and used to predict the three-dimensional orientation distribution of complex suspension flows.

Keywords:
Fiber suspension; Orientation distribution; Fokker-Planck equation; Finite volume; Planar flow

INTRODUCTION

Fiber suspension flows are widely found in a variety of industrial applications, such as processing of thermoplastics, paper and pulp industries, etc. The final product is dependent on the microstructure of suspensions. Generally, a single fiber is affected by ambient fluid and nearby fibers. As a consequence, it changes the bulk flow rheological properties and leads to different macroscopic material behaviors. For instance, fiber orientation exerts potential effects on the physical properties of thermoplastic parts during molding, and distinctly determines the subsequent stiffness, thermal conductivity, and mechanical behavior (Givler et al., 1983; Folgar and Tucker, 1984; Advani and Tucker, 1987).

For molding parts and paper sheets, the geometrical configuration is always very slim. By using this feature, it is convenient to approximate the fiber suspension as a two-dimensional or planar flow. This approximation is widely recognized as the Hele-Shaw model, which was first adopted by Altan et al. (1990) to predict the three-dimensional orientation within channel and contraction flows. To improve computational efficiency, the quadratic closure approximation was used to close the orientation tensor or the moment of orientations originally proposed by Advani and Tucker (1987,1990). The closure model approximates the higher order orientation tensor by lower order tensors, such as assuming <pppp> ≈ <pp> <pp>. Similarly, a variety of candidates have been developed and tested by Cintra and Tucker (1995), Larson and Doi (1991), Chung and Kwon (2001), Bay (1991), and VerWeyst (1998). The quadratic closure, the orthogonal closure, the natural closure, and the hybrid closure are most commonly used by Larson (1999). Nevertheless, a closure model universally accurate for all types of flows is still not available.

In fact, fiber suspensions belong to a class of liquid type materials, which approach the equilibrium state slowly compared to the changes within their microstructural configuration. Thus, the macroscopic properties are commonly induced by the flow. To describe this type of field configuration, the evolution of the probability distribution function is commonly used and denoted by the Fokker-Planck equation (Öttinger, 1996). However, the transient analytical solution to this function only exists when the Péclet number is infinite (Pe=∞). Clearly, developing numerical techniques to solve this distribution function is of great importance.

The distribution function provides a complete description of the orientation state. Therefore, the memory requirement is huge and the computation is time-consuming, especially in the early 1990s (Bay, 1991; Akbar and Altan, 1992). Consequently, the reduced two-dimensional model has been applied for a large variety of flows (Olson et al., 2004; Krochak et al., 2009; Chinesta et al., 2003). For a three-dimensional solution, there are only a few researches that have concerned the numerical methods, where the finite difference (Advani and Tucker, 1987; Kamal and Mutel, 1989) and the finite element (Strand et al., 1987; Han and Im, 1999) have mainly been employed. However, the solving procedure and the Péclet number range was not mentioned in Advani and Tucker 1987). The maximum Péclet number reported in Kamal and Mutel (1989) is 2×10⁴ and the maximum numbers are 60 and 100 in Strand et al. (1987) and Han and Im (1999), respectively.

With modern computers, the numerical technique for directly solving the Fokker-Planck equation is feasible and necessary. In addition, in the past two decades, the finite volume method has been widely applied in molding process simulation and property analysis (Férec et al., 2008; Pei et al., 2012). In this paper, a numerical scheme based on the finite volume method is presented to solve the three-dimensional orientation distribution function. The numerical algorithm is tested by comparison with the analytical solution and data available in the literature, which covers the full range of Péclet number.

THEORETICAL FOUNDATION

Fiber Orientation Distribution Model

Using spherical coordinates, the fiber orientation is represented by a unit vector, i.e., p ={p₁, p₂, p₃ = (sinθcosΦ,sinθsinΦ,cosθ}, as shown in Figure 1.

Figure 1
Schematic of a unit vector p.

With millions of short fibers immersed in the suspensions, it is impractical to trace each fiber and to store its orientation state. An alternative is to adopt the stochastic process, where the fiber interactions are described by random collisions. On this basis, a statistical model was proposed by Folgar and Tucker (1984). Making use of this model, a probability distribution function Ψ(p,t) is defined to denote the probability to find a fiber located within the region [p,p+dp] at time t, and it can be written as

(1) $\begin{array}{c}\mathit{P}\left({\mathit{\theta }}_0\le \mathit{\theta }\le {\mathit{\theta }}_0+\mathit{d}\mathit{\theta },{\mathit{\phi }}_0\le \mathit{\phi }\le {\mathit{\phi }}_0+\mathit{d}\mathit{\phi }\right)\hfill \\ =\mathit{\psi }\left({\mathit{\theta }}_0,{\mathit{\phi }}_0\right)\sin {\mathit{\theta }}_0\mathit{d}\mathit{\theta }\mathit{d}\mathit{\phi }\hfill \end{array}$

The value of P is in [0,1]. Integrating the probability distribution function over the orientation space, a normalization condition is formed as:

(2) $∲\mathit{\psi }\left(\mathbf{p},\mathit{t}\right)\mathit{d}\mathbf{p}={\int }_0^{\mathit{\pi }}{\int }_0^{2\mathit{\pi }}\mathit{\psi }\left(\mathit{\theta },\mathit{\phi }\right)\sin \mathit{\theta }\mathit{d}\mathit{\theta }\mathit{d}\mathit{\phi }=1$

Using the conservation principle in the orientation space, the probability distribution function must satisfy the continuity equation as

(3) $\frac{\partial \mathit{\psi }}{\partial \mathit{t}}+\nabla ·\left(\dot{\mathbf{p}}\mathit{\psi }\right)=0$

where ▽ is the gradient operator which refers to $\frac{\partial }{\partial \mathbf{p}}$, ṗ is the change rate of p and it can be written as:

(4) $\dot{\mathbf{p}}=-\mathit{\omega }·\mathbf{p}+\mathit{\lambda }\left(\mathit{E}·\mathbf{p}-\mathit{E}:\mathbf{ppp}\right)$

where $\mathit{\omega }=\frac{1}{2}\left(\nabla \mathit{u}{}^{\mathit{T}}-\nabla \mathit{u}\right)$ is the vorticity tensor, $\mathit{E}=\frac{1}{2}\left(\nabla {\mathit{u}}^{\mathit{T}}+\nabla \mathit{u}\right)$ is the strain rate tensor λ = (r²_p - 1) / (r²_p + 1), is the fiber shape factor, and r_p = L / d is the fiber aspect ratio, i.e., the ratio of fiber's length L to its diameter d. Eq. (4) is the well-known Jeffery equation (Jeffery, 1922).

To model the fiber-fiber interactions or Brownian motion in fiber suspensions, a random term $\frac{{\mathit{D}}_{\mathit{r}}}{\mathit{\psi }}\nabla \mathit{\psi }$ is inserted into the right hand side of Eq. (4), and it becomes (Folgar and Tucker, 1984; Doi and Edwards, 1988; Dinh and Armstrong, 1984):

(5) ${\dot{\mathbf{p}}}_{\mathit{D}}=\mathit{\omega }·\mathbf{p}+\mathit{\lambda }\left(\mathit{E}·\mathbf{p}-\mathit{E}:\mathbf{ppp}\right)-\frac{{\mathit{D}}_{\mathit{r}}}{\mathit{\psi }}\nabla \mathit{\psi }$

where D_r is the rotary diffusion tensor and here is considered to be isotropic and replaced by a scalar D_r = C_i|ẏ| which is called the rotary diffusivity (Folgar and Tucker, 1984), C_i is a constant describing the random interaction between fibers, ẏ = ▽u^T + ▽u = 2E is the shear rate tensor, and $\left|\dot{\mathit{\gamma }}\right|=\sqrt{\frac{1}{2}\dot{\mathit{\gamma }}:\dot{\mathit{\gamma }}}$ is the effective shear rate.

Putting Eq. (5) into Eq. (3), the continuity equation can be written as:

(6) $\frac{\partial \mathit{\psi }}{\partial \mathit{t}}=-\nabla ·\left(\dot{\mathbf{p}}\mathit{\psi }\right)+{\mathit{D}}_{\mathit{r}}{\nabla }^2\mathit{\psi }$

where ${\nabla }^2=\frac{{\partial }^2}{\partial {\mathbf{p}}^2}$ is the Laplace operator. Eq. (6) is the Fokker-Planck equation in fiber suspensions.

Constitutive Model of Fiber Suspensions

Based on the probability distribution function, it is convenient to apply statistical mechanics to the rheology of fiber suspensions. The macroscopic quantities of suspensions are expressed as the ensemble average of the corresponding quantities. For instance, for a tensor A(p), one has:

(7) $⟨A\left(p\right)⟩=\oint A\left(p\right)\mathit{\psi }\left(p,t\right)dp$

where <A(p)> is the macroscopic value of A(p). The moment tensor of orientation is a valuable property for the fiber suspensions. When $A\left(p\right)\mathbf{=}\underset{n}{\underbrace{p\dots p}},$, the orientation tensor is:

(8) $⟨\underset{\mathit{n}}{\underbrace{\mathbf{p}\dots \mathbf{p}}}⟩=\oint \underset{\mathit{n}}{\underbrace{\mathbf{p}\dots \mathbf{p}}}\mathit{\psi }\left(\mathbf{p},\mathit{t}\right)\mathit{d}\mathbf{p}$

Based on the spherical symmetry, it is easy to write Ψ(p,t) = Ψ(p,t), that is, the function Ψ(p,t) is an even function. For odd numbers of n, one has

(9) $⟨\underset{\mathit{n}}{\underbrace{\mathbf{p}\dots \mathbf{p}}}⟩=0,\mathit{n}=1,3,5\dots$

In the viewpoint of continuum mechanics, the fiber suspensions can be treated as the transversely isotropic fluid. On this basis, the constitutive equation has been established by Ericksen (1960) and Hand (1962). Furthermore, the relationship between the stress and the rate of strain for dilute fiber suspensions has been described by Batchelor (1970) as:

(10) $\mathit{\sigma }=\mathit{\mu }\dot{\mathit{\gamma }}+\frac{1}{2}{\mathit{\mu }}_{\mathit{f}}\dot{\mathit{\gamma }}\left[⟨\mathbf{pppp}⟩-\frac{1}{3}\mathit{I}⟨\mathbf{pppp}⟩\right]$

where µ is the viscosity of the Newtonian fluid and µ_f is a function of fiber concentration, I is the second order unit tensor, and <pppp> and <pp> are, respectively, the fourth and second order orientation tensors.

Calculation of Orientation Tensors

Though the most complete description is presented in the probability distribution function, only the orientation tensors are required to be known for the constitutive relationship Eq. (10). Moreover, the orientation tensors are symmetric, i.e., there are five independent components for <pp> (using the normalization condition: p₁² + p₂² + p²₃ = 1). Similarly, there are forty four independent components for <pppp>.

Generally, the orientation tensors are calculated by the ensemble average with Eq. (8), where the probability distribution function is determined by solving the Fokker-Planck equation. An alternative method is to solve the evolution equation for the orientation tensors. Multiplying Eq. (6) by the tensor A(p) and integrating it:

(11) $\oint \mathbf{A}\left(\mathbf{p}\right)\frac{\partial \mathit{\psi }}{\partial \mathit{t}}\mathit{d}\mathbf{p}=-\oint \mathbf{A}\left(\mathbf{p}\right)\nabla ·\left(\left(-\mathit{\omega }·\mathbf{p}+\mathit{\lambda }\left(\mathit{E}·\mathbf{p}-\mathit{E}:\mathbf{ppp}\right)\right)\mathit{\psi }\right)\mathit{d}\mathbf{p}+\oint \mathbf{A}\left(\mathbf{p}\right){\mathit{D}}_{\mathit{r}}{\nabla }^2\mathit{\psi }\mathit{d}\mathbf{p}$

Taking the time partial derivative $\frac{\partial }{\partial \mathit{t}}$ over the integral term ∮A(p)Ψdp dp, and noting the moment equation Eq. (7), Eq. (10) becomes:

(12) $\frac{\partial }{\partial \mathit{t}}⟨\mathbf{A}\left(\mathbf{p}\right)⟩=-\oint \mathbf{A}\left(\mathbf{p}\right)\nabla ·\left(\left(-\mathit{\omega }·\mathbf{p}+\mathit{\lambda }\left(\mathit{E}·\mathbf{p}-\mathit{E}:\mathbf{ppp}\right)\right)\mathit{\psi }\right)\mathit{d}\mathbf{p}+\oint \mathbf{A}\left(\mathbf{p}\right){\mathit{D}}_{\mathit{r}}{\nabla }^2\mathit{\psi }\mathit{d}\mathbf{p}$

When A(p) = pp, and integrating the right hand side of Eq. (12) by parts, the evolution equation for the second order tensor is:

(13) $\frac{\partial }{\partial \mathit{t}}⟨\mathbf{pp}⟩=-\frac{1}{2}\left(\mathit{\omega }·⟨\mathbf{pp}⟩-⟨\mathbf{pp}⟩·\mathit{\omega }\right)+\frac{\mathit{\lambda }}{2}\left(\dot{\mathit{\gamma }}·⟨\mathbf{pp}⟩+⟨\mathbf{pp}⟩·\dot{\mathit{\gamma }}-2\dot{\mathit{\gamma }}:⟨\mathbf{pppp}⟩\right)+2{\mathit{D}}_{\mathit{r}}\left(\mathit{I}-3⟨\mathbf{pp}⟩\right)$

where the new term <pppp> appears, and an approximation is needed to close this equation. Using subscripts and denoting <pp> and <pppp> as a_ij and a_ijkl, respectively. The quadratic closure approximation is written as:

(14) $\widetilde{{\mathit{a}}_{\mathit{ijkl}}}={\mathit{a}}_{\mathit{ij}}{\mathit{a}}_{\mathit{kl}}$

This model is rather simple and only accurate for a fully aligned orientation distribution. In addition, the linear approximation is

(15) $\widetilde{{\mathit{a}}_{\mathit{ijkl}}}=-\frac{1}{35}\left({\mathit{\delta }}_{\mathit{ij}}{\mathit{\delta }}_{\mathit{kl}}+{\mathit{\delta }}_{\mathit{ik}}{\mathit{\delta }}_{\mathit{jl}}+{\mathit{\delta }}_{\mathit{il}}{\mathit{\delta }}_{\mathit{jk}}\right)+\frac{1}{7}\left({\mathit{a}}_{\mathit{ij}}{\mathit{\delta }}_{\mathit{kl}}+{\mathit{a}}_{\mathit{ik}}{\mathit{\delta }}_{\mathit{jl}}+{\mathit{a}}_{\mathit{il}}{\mathit{\delta }}_{\mathit{jk}}+{\mathit{a}}_{\mathit{kl}}{\mathit{\delta }}_{\mathit{ij}}+{\mathit{a}}_{\mathit{jl}}{\mathit{\delta }}_{\mathit{ik}}+{\mathit{a}}_{\mathit{jk}}{\mathit{\delta }}_{\mathit{il}}\right)$

The linear closure approximation is exact for an isotropic orientation distribution. To circumvent these limitations, many sophisticated models, such as orthogonal closure (Cintra and Tucker, 1995; Chung and Kwon, 2001), natural closure (Verleye and Dupret, 1993), etc., have been developed. The closure model description is compact, general and easy to calculate, which is more widely used for fiber orientation characterization (Advani and Tucker, 1987). As we have seen, the closures are formulated by polynomials and their coefficients are determined by fitting data obtained from a variety of flow fields. In practice, users should be aware of the model characteristics and its applicable parameter range.

NUMERICAL METHOD

Dimensionless Form of the Fokker-Planck Equation

Substituting the nondimensional terms ẏ^* = ẏ/|ẏ|, ω^* = ω/|y| and τ = tD_r into Eq. (6), the dimensionless continuity equation can be written as:

(16) $\frac{\partial \mathit{\psi }}{\partial \mathit{\tau }}=-\mathit{Pe}\nabla ·\left(\dot{\mathbf{p}}*\mathit{\psi }\right)+{\nabla }^2\mathit{\psi }$

where the Péclet number $\mathrm{Pe}=\frac{\left|\dot{\mathit{\gamma }}\right|}{{\mathit{D}}_{\mathit{r}}}$ is defined as the ratio of the effective shear rate to the rotary diffusivity.

Finite Volume Discretization

Integrating Eq. (16) over the orientation space, the integral equation is:

(17) ${\int }_{\mathit{V}}\frac{\partial \mathit{\psi }}{\partial \mathit{\tau }}\mathit{dV}+\mathit{Pe}{\int }_{\mathit{V}}\nabla ·\left(\dot{\mathbf{p}}*\mathit{\psi }\right)\mathit{dV}={\int }_{\mathit{V}}{\nabla }^2\mathit{\psi }\mathit{dV}$

The integral equation is further discretized by the finite volume method. In spherical coordinates, the basic control volume is depicted in Figure 2.

Figure 2
Spherical coordinates and the control volume.

Using Gauss theorem, the convective-diffusive type integral equation is simplified to the face integral equation, and can be written as:

(18) $\begin{array}{c}\frac{{\mathit{\psi }}^{\mathit{\tau }+\Delta \mathit{\tau }}-{\mathit{\psi }}^{\mathit{\tau }}}{\Delta \mathit{\tau }}\sin \mathit{\theta }\Delta \mathit{\theta }\Delta \mathit{\phi }+\mathit{Pe}{\oint }_{\mathit{s}}\left(\dot{\mathbf{p}}*\mathit{\psi }\right)·\mathit{ndS}\hfill \\ ={\oint }_{\mathit{s}}\frac{\partial \mathit{\psi }}{\partial \mathit{n}}\mathit{ds}\hfill \end{array}$

where n is the unit vector normal to the sphere surface, and its total differential is:

(19) $\mathit{dn}={\mathit{dre}}_1+\mathit{rd}\mathit{\theta }{\mathit{e}}_2+\mathit{r}\sin \mathit{\theta }\mathit{d}\mathit{\phi }{\mathit{e}}_3$

where e₁ (r), e₂ (θ),and e₃ (Φ) are the unit differential components in the spherical coordinates. For a unit spherical surface, Eq. (19) is reduced to

(20) $\mathit{dn}=\mathit{d}\mathit{\theta }{\mathit{e}}_2+\sin \mathit{\theta }\mathit{d}\mathit{\phi }{\mathit{e}}_3$

In addition, the gradient of the scalar Ψ is given by:

(21) $\frac{\partial \mathit{\psi }}{\partial \mathit{n}}=\frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}{\mathit{e}}_2+\frac{1}{\sin \mathit{\theta }}\frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}{\mathit{e}}_3$

The discretization is carried out on a unit sphere surface, and mapped onto the plane as shown in Figure 3. A cell-centered discretization scheme is adopted and the unknown variables are stored at the solid points as shown in Figure 3(a).

Figure 3
Schematic of the discretized control volume, and e, w, s and n represent the faces in east, west, south and north directions, respectively.

Making use of Eq. (19)-Eq. (21), the convective part in Eq. (18) is discretized as:

(22) $\begin{array}{c}{\oint }_{\mathit{s}}\left(\dot{\mathbf{p}}*\mathit{\psi }\right)·\mathit{nds}=\left(\dot{\mathit{\theta }\mathit{\psi }}\right)\sin \mathit{\theta }\Delta \mathit{\phi }{|}_{\mathit{e}}-\left(\dot{\mathit{\theta }\mathit{\psi }}\right)\sin \mathit{\theta }\Delta \mathit{\phi }{|}_{\mathit{w}}\\ +\left(\dot{\mathit{\phi }\mathit{\psi }}\sin {\mathit{\theta }}_{\mathit{p}}\right)\Delta \mathit{\theta }{|}_{\mathit{n}}\\ -\left(\dot{\mathit{\phi }\mathit{\psi }}\sin {\mathit{\theta }}_{\mathit{p}}\right)\Delta \mathit{\theta }{|}_{\mathit{s}}\end{array}$

Using the notations of convective flux, the internal flux through each face is:

(23) $\begin{array}{c}{\mathit{F}}_{\mathit{e}}=\left(\mathit{Pe}\right)\dot{\mathit{\theta }}\sin \mathit{\theta }\Delta \mathit{\phi }{|}_{\mathit{e}},{\mathit{F}}_{\mathit{w}}=\left(\mathit{Pe}\right)\dot{\mathit{\theta }}\sin \mathit{\theta }\Delta \mathit{\phi }{|}_{\mathit{w}},\hfill \\ {\mathit{F}}_{\mathit{n}}=\left(\mathit{Pe}\right)\dot{\mathit{\phi }}\sin {\mathit{\theta }}_{\mathit{p}}\Delta \mathit{\theta }{|}_{\mathit{n}},{\mathit{F}}_{\mathit{s}}=\left(\mathit{Pe}\right)\dot{\mathit{\phi }}\sin {\mathit{\theta }}_{\mathit{p}}\Delta \mathit{\theta }{|}_{\mathit{s}}\hfill \end{array}$

Deriving from Eq. (4), the velocities of θ and Φ are respectively written as:

(24) $\dot{\mathit{\theta }}=\frac{\mathit{\lambda }}{2}\sin 2\mathit{\theta }\left[\begin{array}{c}\left(\frac{\partial {\mathit{\nu }}_{\mathit{x}}}{\partial \mathit{x}}{\cos }^2\mathit{\phi }+\frac{\partial {\mathit{\nu }}_{\mathit{y}}}{\partial \mathit{y}}{\sin }^2\mathit{\phi }\right)\hfill \\ +\frac{1}{2}\sin 2\mathit{\phi }\left(\frac{\partial {\mathit{\nu }}_{\mathit{x}}}{\partial \mathit{y}}+\frac{\partial {\mathit{\nu }}_{\mathit{y}}}{\partial \mathit{x}}\right)\hfill \end{array}\right]$

(25) $\begin{array}{c}\dot{\mathit{\phi }}=\frac{\mathit{\lambda }}{2}\left[\sin 2\mathit{\phi }\left(\frac{\partial {\mathit{\nu }}_{\mathit{y}}}{\partial \mathit{y}}-\frac{\partial {\mathit{\nu }}_{\mathit{x}}}{\partial \mathit{x}}\right)+\cos 2\mathit{\phi }\left(\frac{\partial {\mathit{\nu }}_{\mathit{x}}}{\partial \mathit{y}}+\frac{\partial {\mathit{\nu }}_{\mathit{y}}}{\partial \mathit{x}}\right)\right]\hfill \\ +\frac{1}{2}\left(\frac{\partial {\mathit{\nu }}_{\mathit{y}}}{\partial \mathit{x}}-\frac{\partial {\mathit{\nu }}_{\mathit{x}}}{\partial \mathit{y}}\right)\hfill \end{array}$

The derivatives of the velocity in Eq. (24) and Eq. (25) are assumed to be known at this stage. And similar to the convective term, the diffusive part in Eq. (18) is discretized as:

(26) $\begin{array}{c}{\oint }_{\mathit{s}}\frac{\partial \mathit{\psi }}{\partial \mathit{n}}\mathit{dS}=\frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}\sin \mathit{\theta }\Delta \mathit{\phi }{|}_{\mathit{e}}-\frac{\partial \mathit{\psi }}{\partial \mathit{\theta }}\sin \mathit{\theta }\Delta \mathit{\phi }{|}_{\mathit{w}}\hfill \\ +\frac{1}{\sin {\mathit{\theta }}_{\mathit{p}}}\frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}\Delta \mathit{\theta }{|}_{\mathit{n}}-\frac{1}{\sin {\mathit{\theta }}_{\mathit{p}}}\frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}\Delta \mathit{\theta }{|}_{\mathit{s}}\hfill \end{array}$

where the internal diffusive flux through each face is:

(27) $\begin{array}{c}{\mathit{D}}_{\mathit{e}}=\sin \mathit{\theta }\frac{\Delta \mathit{\phi }}{\Delta \mathit{\theta }}{|}_{\mathit{e}},{\mathit{D}}_{\mathit{w}}=\sin \mathit{\theta }\frac{\Delta \mathit{\phi }}{\Delta \mathit{\theta }}{|}_{\mathit{w}},\hfill \\ {\mathit{D}}_{\mathit{n}}=\frac{1}{\sin {\mathit{\theta }}_{\mathit{p}}}\frac{\Delta \mathit{\theta }}{\Delta \mathit{\phi }}{|}_{\mathit{n}},{\mathit{D}}_{\mathit{s}}=\frac{1}{\sin {\mathit{\theta }}_{\mathit{p}}}\frac{\Delta \mathit{\theta }}{\Delta \mathit{\phi }}{|}_{\mathit{s}}\hfill \end{array}$

Then, putting Eq. (22)-Eq. (26) into Eq. (18), and the discrete algebraic equation can be written as:

(28) $\begin{array}{c}\left({\mathit{\psi }}_{\mathit{P}}^{\mathit{\tau }+\Delta \mathit{\tau }}-{\mathit{\psi }}_{\mathit{P}}^{\mathit{\tau }}\right){\mathit{a}}_{\mathit{\tau }}+\mathit{\alpha }{\mathit{a}}_{\mathit{P}}{\mathit{\psi }}_{\mathit{P}}^{\mathit{\tau }+\Delta \mathit{\tau }}+\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{P}}{\mathit{\psi }}_{\mathit{P}}^{\mathit{\tau }}\hfill \\ -\mathit{\alpha }{\left({\mathit{a}}_{\mathit{E}}{\mathit{\psi }}_{\mathit{E}}+{\mathit{a}}_{\mathit{W}}{\mathit{\psi }}_{\mathit{W}}+{\mathit{a}}_{\mathit{N}}{\mathit{\psi }}_{\mathit{N}}+{\mathit{a}}_{\mathit{S}}{\mathit{\psi }}_{\mathit{S}}\right)}^{\mathit{\tau }+\Delta \mathit{\tau }}\hfill \\ =\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{E}}{\mathit{\psi }}_{\mathit{E}}^{\mathit{\tau }}+\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{W}}{\mathit{\psi }}_{\mathit{W}}^{\mathit{\tau }}+\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{N}}{\mathit{\psi }}_{\mathit{N}}^{\mathit{\tau }}\hfill \\ +\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{S}}{\mathit{\psi }}_{\mathit{S}}^{\mathit{\tau }}\hfill \end{array}$

After rearrangement, Eq. (28) becomes:

(29) $\begin{array}{c}\left({\mathit{a}}_{\mathit{\tau }}+\mathit{\alpha }{\mathit{a}}_{\mathit{P}}\right){\mathit{\psi }}_{\mathit{P}}^{\mathit{\tau }+\Delta \mathit{\tau }}-\mathit{\alpha }{\mathit{a}}_{\mathit{E}}{\mathit{\psi }}_{\mathit{E}}^{\mathit{\tau }+\Delta \mathit{\tau }}\hfill \\ -\mathit{\alpha }{\mathit{a}}_{\mathit{W}}{\mathit{\psi }}_{\mathit{W}}^{\mathit{\tau }+\Delta \mathit{\tau }}-\mathit{\alpha }{\mathit{a}}_{\mathit{N}}{\mathit{\psi }}_{\mathit{N}}^{\mathit{\tau }+\Delta \mathit{\tau }}-\mathit{\alpha }{\mathit{a}}_{\mathit{s}}{\mathit{\psi }}_{\mathit{S}}^{\mathit{\tau }+\Delta \mathit{\tau }}\hfill \\ =\left(1-\mathit{\alpha }\right){\mathit{\alpha }}_{\mathit{E}}{\mathit{\psi }}_{\mathit{E}}^{\mathit{\tau }}+\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{W}}{\mathit{\psi }}_{\mathit{W}}^{\mathit{\tau }}+\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{N}}{\mathit{\psi }}_{\mathit{N}}^{\mathit{\tau }}\hfill \\ +\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{S}}{\mathit{\psi }}_{\mathit{S}}^{\mathit{\tau }}+\left[{\mathit{a}}_{\mathit{\tau }}-\left(1-\mathit{\alpha }\right){\mathit{a}}_{\mathit{P}}\right]{\mathit{\psi }}_{\mathit{P}}^{\mathit{\tau }}\hfill \end{array}$

where a_E, a_W, a_S, a_N and a_P are coefficients listed in the ^Appendix. When α = 1/2, the Crank-Nicolson scheme is formed. To accelerate the convergence, the dimensionless time step is multiplied by a coefficient α_τ, i.e., Δτ = α_τ × D_r, such that the time step can adapt to the varied Pe number.

Boundary Condition Implementation

On the sphere boundary, there are two pole points where singularity is always encountered because of the term $\frac{1}{\sin {\mathit{\theta }}_{\mathit{p}}}$, while in the finite volume discretization, this term is not directly evaluated on the boundary. Alternatively, the flux is evaluated on the boundary instead of its value. In addition, it is easy to see that the flow area sinθΔθΔΦ equals zero when θ = 0 and θ = π, thus there is no flux through the boundaries of the two pole points. In the finite difference discretization, the condition Ψ(θ = 0,Φ) = Ψ(θ = π,Φ) = 0 has been imposed on the boundary (Advani and Tucker, 1987; Kamal and Mutel, 1989), which means that there is no possibility that one fiber is located at these two pole points. Obviously, this condition is not always correct. Furthermore, the finite difference scheme does not guarantee the local flux conservation (Férec et al., 2008; Bay and Tucker, 1992). A finite volume scheme was provided by Bay (1991), Férec et al. (2008), Lin et al. (2010), Bay and Tucker (1992), but the implementation of the boundary condition was not detailed.

In practice, when Φ = 0 and Φ = 2π, the boundaries are joined together; thus, the spherical symmetrical boundary condition is written as:

(30) $\mathit{\psi }\left(\mathit{\theta },0\right)=\mathit{\psi }\left(\mathit{\theta },2\mathit{\pi }\right)$

Assuming that the initial fiber orientation is isotropic, and using Eq. (2), one has:

(31) $\mathit{\psi }\left(\mathit{\theta },\mathit{\phi }\right){\int }_0^{\mathit{\pi }}{\int }_0^{2\mathit{\pi }}\sin \mathit{\theta }\mathit{d}\mathit{\theta }\mathit{d}\mathit{\phi }=1$

That is, Ψ(θ,Φ) = 1/4π.

The CTDMA Algorithm

Furthermore, Eq. (29) is discretized into a tridiagonal matrix form, and at the control volume center (i, j), it can be written as:

(32) ${\mathit{A}}_{\mathit{i}}{,}_{\mathit{j}}{\mathit{\psi }}_{\mathit{i},\mathit{j}}={\mathit{B}}_{\mathit{i},\mathit{j}}{\mathit{\psi }}_{\mathit{i},\mathit{j}+1}+{\mathit{C}}_{\mathit{i},\mathit{j}}{\mathit{\psi }}_{\mathit{i},\mathit{j}-1}+{\mathit{D}}_{\mathit{i},\mathit{j}}$

where A_i,j = a_τ + αa_P, B_i,j = αa_N, C_i,j = αa_S, and the other terms in Eq. (29) enter into D_i,j. Similarly, at the cell center (i, j-1), assuming Ψ_i,j-1 = E_i,j-1Ψi,j + F_i,j-1Ψ_i,1+G_i,j-1 and inserting it into Eq. (32), and after rearrangement, the equation becomes:

(33) $\begin{array}{c}\left({\mathit{A}}_{\mathit{i},\mathit{j}}-{\mathit{C}}_{\mathit{i},\mathit{j}}{\mathit{E}}_{\mathit{i},\mathit{j}-1}\right){\mathit{\psi }}_{\mathit{i},\mathit{j}}={\mathit{B}}_{\mathit{i},\mathit{j}}{\mathit{\psi }}_{\mathit{i},\mathit{j}+1}+{\mathit{C}}_{\mathit{i},\mathit{j}}{\mathit{F}}_{\mathit{i},\mathit{j}-1}{\mathit{\psi }}_{\mathit{i},1}\\ +{\mathit{C}}_{\mathit{i},\mathit{j}}{\mathit{G}}_{\mathit{i},\mathit{j}-1}+{\mathit{D}}_{\mathit{i},\mathit{j}}\end{array}$

And rewriting Eq. (33) as:

(34) ${\mathit{\psi }}_{\mathit{i},\mathit{j}}={\mathit{E}}_{\mathit{i},\mathit{j}}{\mathit{\psi }}_{\mathit{i},\mathit{j}+1}+{\mathit{F}}_{\mathit{i},\mathit{j}}{\mathit{\psi }}_{\mathit{i},1}+{\mathit{G}}_{\mathit{i},\mathit{j}}$

Taking A_i,j - C_i,jE_i,j-1 = U_i,j, an iterative relationship is obtained as:

(35) $\begin{array}{c}{\mathit{E}}_{\mathit{i},\mathit{j}}=\frac{{\mathit{B}}_{\mathit{i},\mathit{j}}}{{\mathit{U}}_{\mathit{i},\mathit{j}}},{\mathit{E}}_{\mathit{i},\mathit{j}}=\frac{{\mathit{C}}_{\mathit{i},\mathit{j}}{\mathit{F}}_{\mathit{i},\mathit{j}-1}}{{\mathit{U}}_{\mathit{i},\mathit{j}}},\hfill \\ {\mathit{E}}_{\mathit{i},\mathit{j}}=\frac{{\mathit{C}}_{\mathit{i},\mathit{j}}{\mathit{G}}_{\mathit{i},\mathit{j}-1}+{\mathit{D}}_{\mathit{i},\mathit{j}}}{{\mathit{U}}_{\mathit{i},\mathit{j}}}\hfill \end{array}$

At the boundary cell center (i, N), Eq. (32) becomes:

(36) ${\mathit{A}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},\mathit{N}}={\mathit{B}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{iN}+1}+{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},\mathit{N}-1}+{\mathit{D}}_{\mathit{i},\mathit{N}}$

And at the internal cell center (i, N-1), Eq. (34) becomes:

(37) ${\mathit{\psi }}_{\mathit{i},\mathit{N}-1}={\mathit{E}}_{\mathit{i},\mathit{N}-1}{\mathit{\psi }}_{\mathit{i},\mathit{N}}+{\mathit{F}}_{\mathit{i},\mathit{N}-1}{\mathit{\psi }}_{\mathit{i},1}+{\mathit{G}}_{\mathit{i},\mathit{N}-1}$

Then inserting Eq. (37) into Eq. (36), leads to:

(38) $\begin{array}{c}\hfill \left({\mathit{A}}_{\mathit{i},\mathit{N}}-{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{E}}_{\mathit{i},\mathit{N}-1}\right){\mathit{\psi }}_{\mathit{i},\mathit{N}}={\mathit{B}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},\mathit{N}+1}\\ \hfill +{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{F}}_{\mathit{i},\mathit{N}-1}{\mathit{\psi }}_{\mathit{i},1}\\ \hfill +{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{G}}_{\mathit{i},\mathit{N}-1}+{\mathit{D}}_{\mathit{i},\mathit{N}}\end{array}$

And similar to Eq. (37), Eq. (38) can be written as

(39) ${\mathit{\psi }}_{\mathit{i},\mathit{N}}={\mathit{E}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},\mathit{N}+1}+{\mathit{F}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},1}+{\mathit{G}}_{\mathit{i},\mathit{N}}$

Substituting Eq. (39) into Eq. (38), the expression for Ψ_i,N+1 becomes:

(40) $\begin{array}{c}\underset{{\mathit{P}}_{\mathit{i},\mathit{N}}}{\underbrace{\left[\left({\mathit{A}}_{\mathit{i},\mathit{N}}-{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{E}}_{\mathit{i},\mathit{N}-1}\right){\mathit{E}}_{\mathit{i},\mathit{N}}-{\mathit{B}}_{\mathit{i},\mathit{N}}\right]}}{\mathit{\psi }}_{\mathit{i},\mathit{N}+1}\hfill \\ =\underset{{\mathit{q}}_{\mathit{i},\mathit{N}}}{\underbrace{\left[{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{F}}_{\mathit{i},\mathit{N}-1}-\left({\mathit{A}}_{\mathit{i},\mathit{N}}-{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{E}}_{\mathit{i},\mathit{N}-1}\right){\mathit{F}}_{\mathit{i},\mathit{N}}\right]{\mathit{\psi }}_{\mathit{i},1}}}\hfill \\ \underset{{\mathit{m}}_{\mathit{i},\mathit{N}}}{\underbrace{+{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{G}}_{\mathit{i},\mathit{N}-1}+{\mathit{D}}_{\mathit{i},\mathit{N}}-\left({\mathit{A}}_{\mathit{i},\mathit{N}}-{\mathit{C}}_{\mathit{i},\mathit{N}}{\mathit{E}}_{\mathit{i},\mathit{N}-1}\right){\mathit{G}}_{\mathit{i},\mathit{N}}}}\hfill \end{array}$

Then taking N+2-i=k, and substituting the subscript i with k, Eq. (40) becomes:

(41) $\left\{\begin{array}{c}{\mathit{p}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},\mathit{N}+1}={\mathit{q}}_{\mathit{i},\mathit{N}}{\mathit{\psi }}_{\mathit{i},1}+{\mathit{m}}_{\mathit{i},\mathit{N}}\hfill \\ {\mathit{p}}_{\mathit{k},\mathit{N}}{\mathit{\psi }}_{\mathit{k},\mathit{N}+1}={\mathit{q}}_{\mathit{k},\mathit{N}}{\mathit{\psi }}_{\mathit{k},1}+{\mathit{m}}_{\mathit{k},\mathit{N}}\hfill \end{array}\right.$

On the boundary surface, Eq. (30) is:

(42) $\left\{\begin{array}{c}{\mathit{\psi }}_{\mathit{i},\mathit{N}+1}={\mathit{\psi }}_{\mathit{k},1}\hfill \\ {\mathit{\psi }}_{\mathit{k},\mathit{N}+1}={\mathit{\psi }}_{\mathit{i},1}\hfill \end{array}\right.$

Multiplying the first equation in Eq. (41) by p_k,N and noting p_k,N Ψ_k,N+1, it leads to:

(43) $\begin{array}{c}{\mathit{p}}_{\mathit{i},\mathit{n}}{\mathit{p}}_{\mathit{k},\mathit{N}}{\mathit{\psi }}_{\mathit{i},\mathit{N}+1}={\mathit{q}}_{\mathit{i},\mathit{N}}{\mathit{q}}_{\mathit{k},\mathit{N}}{\mathit{\psi }}_{\mathit{k},1}+{\mathit{m}}_{\mathit{i},\mathit{N}}{\mathit{p}}_{\mathit{k},\mathit{N}}\\ +{\mathit{m}}_{\mathit{k},\mathit{N}}{\mathit{p}}_{\mathit{i},\mathit{N}}\end{array}$

Thus, the boundary value is written as

(44) ${\mathit{\psi }}_{\mathit{i},\mathit{N}+1}=\frac{{\mathit{m}}_{\mathit{i},\mathit{N}}{\mathit{p}}_{\mathit{k},\mathit{N}}+{\mathit{m}}_{\mathit{k},\mathit{N}}{\mathit{q}}_{\mathit{i},\mathit{N}}}{{\mathit{p}}_{\mathit{i},\mathit{N}}{\mathit{p}}_{\mathit{k},\mathit{N}}-{\mathit{q}}_{\mathit{i},\mathit{N}}{\mathit{q}}_{\mathit{k},\mathit{N}}}$

WithEq. (32), Eq. (34), Eq. (35), and using Eq. (44), the algebraic equation system can be iteratively solved by the LU decomposition, which has also been denoted as the cyclic tridiagonal-matrix algorithm (C-TDMA) (Férec et al., 2008; Bay and Tucker, 1992).

NUMERICAL CALCULATION AND VALIDATION

Test of Grid Relevance

To solve the time evolution of the probability distribution function, there were 57 grids equally spaced in Φ for the planar orientation (Advani and Tucker, 1987). The 39×77 grids were used in Advani and Tucker (1990), and the 40×40 grids were used in Férec et al. (2008) for the three-dimensional orientation. But the accuracy validation, especially the grid relevance, has not been tested. In addition, for the transient Fokker-Planck, the calculation is computationally-intensive using the high grid resolution. Alternatively, the steady Fokker-Planck equation is solved. Using Eq. (29) and noting α = 1, the steady discretization equation system is formed as:

(45) ${\mathit{a}}_{\mathit{P}}{\mathit{\psi }}_{\mathit{P}}-{\mathit{a}}_{\mathit{E}}{\mathit{\psi }}_{\mathit{E}}-{\mathit{a}}_{\mathit{W}}{\mathit{\psi }}_{\mathit{W}}-{\mathit{a}}_{\mathit{N}}{\mathit{\psi }}_{\mathit{N}}-{\mathit{a}}_{\mathit{S}}{\mathit{\psi }}_{\mathit{S}}=0$

The calculation is tested for different grid resolutions when Pe=100, and the results are listed in Table 1.

Table 1
Grid relevance test when Pe=100.

Grids/ a2	a 11	a 12	a 22
36×72	0.730462	0.080243	0.094655
72×144	0.728071	0.076282	0.092953
144×288	0.727373	0.075166	0.092479
288×576	0.727188	0.074868	0.092354
1000×2000	0.727162	0.074659	0.091948

These tests were carried out on a laptop with Intel T7300 CPU (2.0GHz) and 2GB RAM. The calculations take less than 1 minute for 36×72 grids, and several minutes to tens of minutes for 72×144, 144×288, 288×576 grids, but exceed ten hours for 1000×2000 grids. From Table 1, the error for 36×72 grids is rather big, so if the computational resources are enough, it is recommended to adopt high resolution grids.

Table 1: Grid relevance test when Pe=100.

Analytical Solution When Pe=∞

When the Pe number is infinite, and the diffusion term is negligible in the Fokker-Planck equation, we obtain

(46) $\frac{\partial \mathit{\psi }}{\partial \mathit{t}}=\nabla ·\left(\dot{\mathbf{p}}\mathit{\psi }\right)$

where the convection is dominant, and the analytical solution for the Fokker-Planck equation is (Dinh and Armstrong, 1984):

(47) $\mathit{\psi }\left(\mathbf{p},\mathit{t}\right)=\frac{1}{4\mathit{\pi }}{\left({\Delta }^{\mathit{t}}·\Delta :\mathbf{pp}\right)}^{-3/2}$

where Δ = E^-1, and E is the deformation gradient, calculated as

(48) $\frac{\partial \mathit{E}}{\partial \mathit{t}}=\left(-\mathit{\omega }·\mathbf{p}+\mathit{\lambda }\mathit{E}·\mathbf{p}\right)·\mathit{E}$

When the fiber shape factor λ is unity, Eq. (48) becomes:

(49) $\frac{\partial \mathit{E}}{\partial \mathit{t}}=\Delta \mathit{u}·\mathit{E}$

In a simple shear flow, {u = ẏy,v = 0, w=0}, and the shear rate is y=ẏt, E can be written as:

$\mathit{E}=\left[\begin{array}{ccc}1& \mathit{\gamma }& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],\mathit{and}\Delta ={\mathit{E}}^{-1}=\left[\begin{array}{ccc}1& -\mathit{\gamma }& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right].$

Using Eq. (47)-Eq. (49), and noting Eq. (8), the time evolution for the orientation tensor is calculated as

(50) ${\mathit{a}}_{\mathit{ij}}=⟨\mathbf{pp}⟩=\frac{1}{4\mathit{\pi }}\oint \mathbf{pp}{\left({\Delta }^{\mathit{t}}·\Delta :\mathbf{pp}\right)}^{-3/2}\mathit{d}\mathbf{p}$

To test the numerical code with high Pe number, the main components a₁₁, a₁₂ and a₂₂ are calculated with 1000×2000 grids and compared to the analytical solutions (Pe=10⁸), as shown in Figure 4.

Figure 4
Comparison of the analytical solutions with Pe=10⁸ (-a₁₁, ......a₁₂, and............a₂₂) and the numerical results (▲ a₁₁, ■ a₁₂, and ● a₂₂) for a simple shear flow.

Comparison with Literature Results

For transient evolution of the orientation distribution function, the calculation time becomes very long for high grid resolution. In the existing literature, the grid resolutions are 39×77 in Advani and Tucker (1990), 40×40 in Férec et al. (2008), or were not mentioned in Lin et al. (2010). So the calculation is carried on 36×72 grids and the results are listed and compared in Figure 5.

Figure 5
Comparison of numerical results adapted from the literature (Férec et al., 2008) (▲ Pe=1000, ■ Pe=100, ● Pe=10, and ♦ Pe=1) and the results in this study (-Pe=1000, ......Pe=100,............Pe=10, and-...-...-Pe=1) for a simple shear flow.

By comparison, the main orientation tensors are in good agreement with the literature. The small difference mainly comes from the grid resolution. For high Pe number, much more time is needed to reach the steady solution, so it manifests a big difference.

CONCLUSION

The objective of this research is to directly solve the orientation distribution function, i.e., the Fokker-Planck equation, which has been widely used in the prediction of fiber orientation distribution in fiber-filled polymers, pulp fiber suspensions, etc. To achieve this objective, a finite volume scheme is proposed. The scheme is characterized as globally conservative, accurate, and stable. Based on the face flux in the FV discretization, the singularity at pole points encountered in the FD method is naturally solved. In addition, the global flux conservation is guaranteed. Then the spherical symmetry boundary condition is numerically formulated, which enhances the scheme's accuracy. Subsequently, the numerical validation was conducted for a simple shear flow, where the results for P_e=∞ are compared with analytical solutions, and the others are compared with data in the literature. Furthermore, this scheme is computationally inexpensive and suitable for the evaluation of the rheological properties of complex fiber suspensions.

ACKNOWLEDGMENT

The authors are grateful for the financial support by the National Natural Science Foundation of China (No. 51309118), the Natural Science Foundation of Jiangsu Province (No. BK20130527), and the Postdoctoral Science foundation of China (No. 2013M531282).

APPENDIX

$\left\{\begin{array}{c}{\mathit{a}}_{\mathit{E}}={\mathit{D}}_{\mathit{e}}·\max \left\{0,{\left(1-0.1\left|{\mathit{Pe}}_{\mathit{e}}\right|\right)}^5\right\}+\max \left(-{\mathit{F}}_{\mathit{e}},0\right),\hfill \\ {\mathit{a}}_{\mathit{W}}={\mathit{D}}_{\mathit{e}}·\max \left\{0,{\left(1-0.1\left|{\mathit{Pe}}_{\mathit{W}}\right|\right)}^5\right\}+\max \left(-{\mathit{F}}_{\mathit{w}},0\right)\hfill \\ {\mathit{a}}_{\mathit{N}}={\mathit{D}}_{\mathit{n}}·\max \left\{0,{\left(1-0.1\left|{\mathit{Pe}}_{\mathit{n}}\right|\right)}^5\right\}+\max \left(-{\mathit{F}}_{\mathit{n}},0\right)\hfill \\ {\mathit{a}}_{\mathit{S}}={\mathit{D}}_{\mathit{s}}·\max \left\{0,{\left(1-0.1\left|{\mathit{Pe}}_{\mathit{s}}\right|\right)}^5\right\}+\max \left(-{\mathit{F}}_{\mathit{s}},0\right)\hfill \\ {\mathit{a}}_{\mathit{P}}={\mathit{a}}_{\mathit{E}}+{\mathit{a}}_{\mathit{W}}+{\mathit{a}}_{\mathit{N}}+{\mathit{a}}_{\mathit{S}}+{\mathit{F}}_{\mathit{e}}-{\mathit{F}}_{\mathit{w}}+{\mathit{F}}_{\mathit{n}}-{\mathit{F}}_{\mathit{s}},\hfill \\ {\mathit{a}}_{\mathit{\tau }}=\frac{\sin \mathit{\theta }\Delta \mathit{\theta }\Delta \mathit{\phi }}{\Delta \mathit{\tau }}\hfill \end{array}\right.\left\{\begin{array}{c}\frac{{\mathit{F}}_{\mathit{e}}}{{\mathit{D}}_{\mathit{e}}}\hfill \\ \frac{{\mathit{F}}_{\mathit{w}}}{{\mathit{D}}_{\mathit{w}}}\hfill \\ \frac{{\mathit{F}}_{\mathit{n}}}{{\mathit{D}}_{\mathit{n}}}\hfill \\ \frac{{\mathit{F}}_{\mathit{s}}}{{\mathit{D}}_{\mathit{s}}}\hfill \end{array}\right.$

REFERENCES

• Advani, S. G. and Tucker, C. L., Closure approximations for three-dimensional structure tensors. J. Rheol., 34, 367 (1990).
• Advani, S. G. and Tucker, C. L., The use of tensors to describe and predict fiber orientation in short fiber composites. J. Rheol., 31, 751 (1987).
• Akbar, S. and Altan, M. C., On the solution of fiber orientation in two-dimensional homogeneous flows. Polym. Eng. Sci., 32, 810 (1992).
• Altan, M. C., Subbiah, S., Güçeri, S. I. and Pipes, R. B., Numerical prediction of three-dimensional fiber orientation in Hele-Shaw flows. Poly. Eng. Sci., 30, 848 (1990).
• Batchelor, G. K., The stress system in a suspension of force-free particles. J. Fluid Mech., 41, 545 (1970).
• Bay, R. S., Fiber orientation in injection molded composites: A comparison of theory and experiment. PhD Thesis, University of Illinois at Urbana-Champaign (1991).
• Bay, R. R. and Tucker, C. L., Fiber orientation in simple injection moldings. 1. Theory and numerical methods. Poly. Com., 13, 317 (1992).
• Chinesta, F., Chaidron, G. and Poitou, A., On the solution of Fokker-Planck equations in steady recirculating flows involving short fiber suspensions. J. Non-Newt. Fluid Mech., 113, 97 (2003).
• Chung, D. H. and Kwon, T. H., Improved model of orthotropic closure approximation for flow induced fiber orientation. Polym. Compos., 22, 636 (2001).
• Cintra, J. S. and Tucker, C. L., Orthotropic closure approximation for flow-induced fiber orientation. J. Rheol., 39, 1095 (1995).
• Dinh, S. M., Armstrong, R. C., A rheological equation of state for semiconcentrated fiber suspensions. J. Rheol., 28, 207 (1984).
• Doi, M., Edwards, S. F., Theory of Polymer Dynamics. Oxford University Press, England (1988).
• Ericksen, J. L., Anisotropic Fluids, Arch. Rat. Mech. Anal., 4, 231 (1960).
• Férec, J., Heniche, M., Heuzey, M. C., Ausias, G. and Carreau, P. J., Numerical solution of the Fokker-Planck equation for fiber suspensions: Application to the Folgar-Tucker-Lipscomb model. J. Non-Newt. Fluid Mech., 155, 20 (2008).
• Folgar, F. and Tucker, C. L., A phonomenological model for concentrated fibre suspensions. J. Rei. Plas. Com., 3, 98 (1984).
• Givler, R. C., Crochet, M. J. and Pipes, R. B., Numerical Prediction of fibre orientation in dilute suspensions. J. Com. Mat., 17, 330 (1983).
• Han, K. H. and Im, Y. T., Modified hybrid closure approximation for prediction of flow-induced fiber orientation. J. Rheol., 43, 569 (1999).
• Hand, G. L., A theory of anisotropic fluids. J. Fluid Mech., 13, 33 (1962).
• Jeffery, G. B., The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. London, A102, 161 (1922).
• Kamal, M. R. and Mutel, A. T., The Prediction of flow and orientation behavior of short fiber reinforced melts in simple flow systems. Poly. Com., 10, 337 (1989).
• Krochak, P. J., Olson, J. A. and Martinez, D. M., Fiber suspension flow in a tapered channel: The effect of flow/fiber coupling. Int. J. Multiphase Flow, 35, 676 (2009).
• Larson, R. G., The Structure and Rheology of Complex Fluids. Oxford University Press, New York (1999).
• Larson, R. G. and Doi, M., Mesoscopic domain theory for textured liquid crystalline polymers. J. Rheol., 35, 539 (1991).
• Lin, J. Z., Zhang, Q. H., Zhang, K., Rheological properties of fiber suspensions flowing through a curved expansion duct. Poly. Eng. Sci., 50, 1994 (2010).
• Olson, J. A., Frigaard, I., Candice, C. and Hämaläine, J. P., Modeling a turbulent fibre suspension flowing in a planar contraction: The one-dimensional headbox. Int. J. Multiphase Flow, 30, 51 (2004).
• Öttinger, H. C., Stochastic Processes in Polymeric Fluids, Tools and Examples for Developing Simulation Algorithms. Springer, Berlin (1996).
• Pei, Z., Hu, B., Diao, C. and Yu, C., Investigation on the motion of different types of fibers in the vortex spinning nozzle. Polym. Eng. Sci., 52, 856 (2012).
• Strand, S. R., Kim, S. and Karrila, S. J., Computation of rheological properties of suspensions of rigid rods: stress growth after inception of steady shear flow. J. Non-Newtonian Fluid Mech., 24, 311 (1987).
• Verleye, V. and Dupret, F., Proceedings of ASME Winter Annual Meeting. ASME, New York, 139 (1993).
• VerWeyst, B. E., Numerical Predictions of flow induced fiber orientation in three-dimensional geometries. PhD Thesis, University of Illinois at Urbana-Champaign (1998).

Publication Dates

• Publication in this collection
  Jan-Mar 2017

History

• Received
  09 Dec 2014
• Reviewed
  06 July 2015
• Accepted
  19 Sept 2015