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.


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 andTucker (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), andVerWeyst (1998).The quadratic closure, the orthogonal closure, the natural closure, and the hybrid closure are most commonly Brazilian Journal of Chemical Engineering 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 threedimensional 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.

Fiber Orientation Distribution Model
Using spherical coordinates, the fiber orientation is represented by a unit vector, i.e.,

   
1 2 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   ,t  p 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 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  p , p  is the change rate of p and it can be written as: where  4) is the wellknown Jeffery equation (Jeffery, 1922).
To model the fiber-fiber interactions or Brownian motion in fiber suspensions, a random term r D    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): where r D is the rotary diffusion tensor and here is considered to be isotropic and replaced by a scalar  which is called the rotary diffusivity (Folgar and Tucker, 1984), i C is a constant describing the random interaction between fibers, 2 is the shear rate tensor, and   is the effective shear rate.
Putting Eq. ( 5) into Eq.(3), the continuity equation can be written as: where 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: where The moment tensor of orientation is a valuable prop-erty for the fiber suspensions.When    the orientation tensor is: Based on the spherical symmetry, it is easy to write 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: 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 t   over the integral term , 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 ij a and ijkl a , respectively.The quadratic closure approximation is written as: 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.into Eq.( 6), the dimensionless continuity equation can be written as:

Dimensionless
where the Péclet number Pe

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 Brazilian Journal of Chemical Engineering Vol. 34, No. 01, pp. 307 -316, January -March, 2017 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   1 e r ,   2 e  ,and   3 e  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 E a , W a , S a , N a and P a 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

Boundary Condition Implementation
On the sphere boundary, there are two pole points where singularity is always encountered because of the term 1 sin 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  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: 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 , , and the other terms in Eq. ( 29) enter into , i j D .Similarly, at the cell center (i, j-1), assuming and inserting it into Eq.

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.
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 Brazilian Journal of Chemical Engineering 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.

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 , and E is the deformation gradient, calculated as When the fiber shape factor  is unity, Eq. ( 48) becomes: , and the shear rate is t     , E can be written as: , and 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 11 a , 12 a and 22 a are calculated with 1000×2000 grids and compared to the analytical solutions (Pe=10 8 ), as shown in Figure 4.
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 fiberfilled 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 Pe=∞ 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.

Figure 1 :
Figure 1: Schematic of a unit vector p.
Form of the Fokker-Planck EquationSubstituting the nondimensional terms * / effective shear rate to the rotary diffusivity.

Figure 2 :
Figure 2: Spherical coordinates and the control volume.

Figure 3 :
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.
, such that the time step can adapt to the varied Pe number.