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^{4} 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*_{1}, *p*_{2}, *p*_{3} = (sinθcosΦ,sinθsinΦ,cosθ}, as shown in Figure 1.

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**+*d***p**] at time *t*, and it can be written as

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

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

where ▽ is the gradient operator which refers to
**ṗ** is the change rate of **p** and it can be written as:

where
*λ = (r ^{2}_{p} - 1) / (r^{2}_{p} + 1)*, is the fiber shape factor, and

*r*is the fiber aspect ratio, i.e., the ratio of fiber's length

_{p}= L / d*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
^{Folgar and Tucker, 1984}; ^{Doi and Edwards, 1988}; ^{Dinh and Armstrong, 1984}):

where *D _{r}* is the rotary diffusion tensor and here is considered to be isotropic and replaced by a scalar

*D*which is called the rotary diffusivity (

_{r}= C_{i}|ẏ|^{Folgar and Tucker, 1984}),

*C*is a constant describing the random interaction between fibers,

_{i}*ẏ = ▽u*is the shear rate tensor, and

^{T}+ ▽u = 2EPutting Eq. (5) into Eq. (3), the continuity equation can be written as:

where

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:

where <**A(p)**> is the macroscopic value of **A(p)**. The moment tensor of orientation is a valuable property for the fiber suspensions. When

Based on the spherical symmetry, it is easy to write *Ψ( p,t) = Ψ(p,t)*, that is, the function

*Ψ(*is an even function. For odd numbers of

**p**,t)*n*, one has

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:

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 _{1}^{2} + p_{2}^{2} + p^{2}_{3} = 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:

Taking the time partial derivative
**∮A(p)**Ψ*d***p**
*d***p**, and noting the moment equation Eq. (7), Eq. (10) becomes:

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

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*, respectively. The quadratic closure approximation is written as:

_{ijkl}

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

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*into Eq. (6), the dimensionless continuity equation can be written as:

_{r}

where the Péclet number

Finite Volume Discretization

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

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

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

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

where e_{1} (r), e_{2} (θ),and e_{3} (Φ) are the unit differential components in the spherical coordinates. For a unit spherical surface, Eq. (19) is reduced to

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

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

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

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

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

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:

where the internal diffusive flux through each face is:

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

After rearrangement, Eq. (28) becomes:

where *a _{E}, a_{W}, a_{S}, a_{N}* and

*a*are coefficients listed in the Appendix. When

_{P}*α = 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
*θΔθΔΦ* 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:

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

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:

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*. Similarly, at the cell center (

_{i,j}*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:

And rewriting Eq. (33) as:

Taking A_{i,j} - C_{i,j}E_{i,j-1} = U_{i,j}, an iterative relationship is obtained as:

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

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

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

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

Substituting Eq. (39) into Eq. (38), the expression for *Ψ _{i,N+1}* becomes:

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

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

Multiplying the first equation in Eq. (41) by *p _{k,N}* and noting p

_{k,N}Ψ

_{k,N+1}, it leads to:

Thus, the boundary value is written as

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:

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

Grids/ a_{2} |
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

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

where *Δ = E ^{-1}*, and

*E*is the deformation gradient, calculated as

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

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

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

To test the numerical code with high Pe number, the main components a_{11}, a_{12} and a_{22} are calculated with 1000×2000 grids and compared to the analytical solutions (Pe=10^{8}), as shown in Figure 4.

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.

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.