# Simulation of the fiber orientation through a finite element approach to solve the Fokker-Planck equation

Nazih Assaad Al Ayoubi, Hugues Digonnet, Luisa Silva, Christophe Binetruy, Thierry Renault, Sebastien Comas-Cardona

► To cite this version:

Nazih Assaad Al Ayoubi, Hugues Digonnet, Luisa Silva, Christophe Binetruy, Thierry Renault, et al. Simulation of the fiber orientation through a finite element approach to solve the Fokker-Planck equation. Journal of Non-Newtonian Fluid Mechanics, 2024, 331, pp.105284. 10.1016/j.jnnfm.2024.105284. hal-04786739

HAL Id: hal-04786739

https://hal.science/hal-04786739v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Full length article

# Simulation of the fiber orientation through a finite element approach to solve the Fokker-Planck equation

Nazih Assaad Al Ayoubi<sup>a</sup>, Hugues Digonnet<sup>a</sup>, Luisa Silva<sup>a,*</sup>, Christophe Binetruy<sup>a</sup>, Thierry Renault<sup>b</sup>, Sebastien Comas-Cardona<sup>a</sup>

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

### Keywords:

Fiber orientation

Probability distribution function

Fokker-Planck equation

Finite element method

Numerical modeling

## ABSTRACT

This work aims to introduce a groundbreaking approach by directly computing the Fokker-Planck equation, providing a mesoscopic scale orientation indicator based on the 2D-probability density function (PDF) of the fibers' orientation state. Unlike conventional methods that rely on pre-averaged quantities and closure approximations, our method offers enhanced accuracy and information preservation. The model's enhanced accuracy can be served as a foundational tool for future studies, enabling the development of comprehensive models describing the fluid-flow coupling problem with precision. Consequently, this advancement facilitates the simulation of real-case scenarios, such as the dynamic motion of fibers during the injection phase of molten thermoplastics within a mold cavity. The novelty of this work lies in its application of the Streamline-Upwind/Petrov-Galerkin (SUPG) finite element method, on both orientation and physical spaces. Our model shows the potential to improve the understanding and prediction of fiber behavior in industrial applications, offering valuable insights into process optimization and design. Implemented within a finite element framework, a comprehensive investigation is conducted into the effects of mesh refinement, time scheme, and time stepping on the computational modeling of the PDF evolution, aiming to strike an optimal balance between model precision and computational efficiency. The validation tests were conducted for the case of simple shear flow to examine the influence of the interaction coefficient $C_I$ and the fiber shape factor $\lambda$ on the resolution of the probability distribution function. The numerical results demonstrate the evolution of fiber orientation over time under Poiseuille flow conditions.

## 1. Introduction

Thermoplastics reinforced with short fibers have gained prominence in engineering, notably in automotive applications, offering a path to weight reduction and cost efficiency while optimizing mechanical attributes [1,2]. The microstructure generated during manufacturing significantly influences the material's properties [3,4]. For injection-molded parts, fiber orientation profoundly impacts mechanical characteristics, including tensile strength [5,6]. Consequently, predicting fiber orientation using numerical methods has been a critical research focus [7-9]. Prioritizing the study of fiber orientation is essential, providing a foundational understanding before delving into the complexities of flow-fiber coupling.

The bibliographic work [10,11] has highlighted the complexity of the phenomena that govern the orientation of fibers. In the pursuit of

predicting and simulating fiber orientation, various orientation evolution models have been proposed across different scales, encompassing meso, micro, and macroscales. These models aim to reveal the intricate rotational dynamics and behaviors exhibited by discontinuous fibers. At the microscale, a singular model described through the Jeffery equation, is employed to explore the individual dynamics of fiber orientation, with intersections among fibers omitted. This study [12] was one of the first to describe the orientation dynamics of a particle, while providing a simplified yet insightful representation of the orientation dynamics of ellipsoidal particles within a flow. Its model is written as:

$$ \dot{\mathbf{p}} = -\omega \cdot \mathbf{p} + \lambda (\mathbf{E} \cdot \mathbf{p} - \mathbf{E} : \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p}) \quad (1) $$

where $\mathbf{p}$ is the unit vector that describes a single fiber orientation, $\dot{\mathbf{p}}$ is the material derivative of $\mathbf{p}$, $\mathbf{E}$ and $\omega$ are the strain rate and vorticity

<footnote><sup>a</sup> Nantes Université, Ecole Centrale de Nantes, CNRS, Institut de Recherche en Génie Civil et Mécanique (GeM), UMR 6183, 1, rue de la Noë, Nantes, F-44000, France</footnote>

<footnote><sup>b</sup> FORVIA, 14 Rue des Petits Bois, Saint-Malo, 35400, France</footnote>

<footnote>* Corresponding author.
E-mail address: luisa.rocha-da-silva@ec-nantes.fr (L. Silva).</footnote>
tensors respectively, and $\lambda$ is the shape factor which depends on the geometry of the fiber:

$$ 
\begin{aligned}
\mathbf{E} &= \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T) \\
\omega &= \frac{1}{2}(\nabla \mathbf{u}^T - \nabla \mathbf{u}) \\
\lambda &= \frac{\beta^2 - 1}{\beta^2 + 1}
\end{aligned}
\qquad (2) $$

where $\mathbf{u}$ is the fluid velocity vector and $\beta$ is the aspect ratio of the fiber (length/diameter). However, its effectiveness hinges on certain assumptions like rigid ellipsoidal shapes, uniform velocity changes, and adhesive fluid-particle contact. As a result, its accuracy might be limited when dealing with intricate flow patterns, complex particle geometries, or particle interactions.

In commercial software, a model at macroscale is employed, involving the use of the second-order orientation tensor $\mathbf{a}_2$ as the principal variable for characterizing fiber orientation. This tensor is defined as the moment of the probability distribution function $\varphi$ [13], which expresses the probability of finding a fiber with an orientation between ($\mathbf{p}$ and $d\mathbf{p}$) at time t, as it must satisfy conditions of periodicity, normalization, and continuity.

$$ P = \int_{\mathbf{p}}^{\mathbf{p}+d\mathbf{p}} \varphi(\mathbf{p})d\mathbf{p} \qquad (3) $$

Furthermore, orientation tensors in (2D) possess intriguing properties, such as they are completely symmetrical ($\mathbf{a}_{xy} = \mathbf{a}_{yx}$) and their trace is unitary ($a_{xx}+a_{yy} = 1$). The diagonal terms of $\mathbf{a}_2$ signify fiber alignment relative to the coordinate directions $\mathbf{x}$ and $\mathbf{y}$, whereas off-diagonal elements quantify orientation distribution asymmetry. Therefore to model particle motion within semi-concentrated solutions on a macroscopic scale, the orientation state is described through a phenomenological approach allowing for the evolution of the orientation tensor. This model, proposed by Folgar and Tucker [14], is represented in the following form:

$$ \frac{D\mathbf{a}_2}{Dt} = (\mathbf{a}_2 \cdot \omega - \omega \cdot \mathbf{a}_2) + \lambda(E \cdot \mathbf{a}_2 + \mathbf{a}_2 \cdot E - 2E : \mathbf{a}_4) + 2C_I \gamma (\mathbf{I} - 3\mathbf{a}_2) \quad (4) $$

where $\gamma$ is the generalized shear rate, and $C_I$ is the interaction coefficient. Nevertheless, this model possesses a drawback, necessitating the application of a closure approximation. It is needed as a dependency of the lower even order fiber orientation tensor on the tensor of the next higher even-order is recurred, as the fourth order tensor in Eq. (4) is a function of a six-order tensor and so forth. Consequently, different closure approximations have been proposed in the literature [8,15–17], but often resulting in inaccuracy and loss of information.

To mitigate the limitations associated with employing models at both the micro and macroscale, the development of a mesoscopic evolution model is necessary to describe the orientation behavior of a fiber population. This model is governed by the Fokker–Planck equation, requiring the solution of the distribution function $\varphi$ which offers a thorough and precise description of the orientation state of a suspension at a given material point.

Advani and Tucker [18] addressed the Fokker–Planck equation through the application of finite difference method and implicit time integration. They employed central differences for $\mathbf{p}$ derivatives and Gauss–Seidel relaxation for computations at each time step. But notably, re-normalization of the PDF $\varphi(\mathbf{p}, t)$ was necessary due to boundary conditions imposing zero values at angles 0 and $\pi$. However, this implies negligible probabilities for fibers oriented at both angles, a simplification that may not always hold.

Chinesta and co-workers [19,20] employed the meshless particle method for discretizing the Fokker–Planck equation dominated by advection. This approach has been extended to handle a convection–diffusion equation where diffusive behavior is described by random motions [21,22]. Nonetheless, its applicability was constrained to cases involving minimal diffusion effects, as the large number of particles

required for solution stabilization rendered this technique impractical for intricate industrial applications.

Galerkin methods employing spherical harmonics as basis functions within orientation space [23–25] demonstrate efficiency and notable performance in situations characterized by a dominant diffusion regime. However, such approaches are less advisable for challenges driven by advection-dominated transport phenomena. This is due to the potential requirement of a considerable number of modes/moments to achieve positivity preservation for probability distribution functions exhibiting sharp peaks [26].

Several studies [27–29] have employed an approximation approach to incorporate Fokker–Planck resolution within a finite volume framework, yielding promising outcomes. However, in scenarios with heightened convection dominance and the presence of sharp gradients, the SUPG finite element method can offer advantages. This method is particularly proficient at managing strongly aligned orientation states, which frequently arise in practical situations involving fibers with quasi-infinite shapes ($\lambda \to 1$). Moreover, FEM-SUPG facilitates the implementation of adaptive mesh refinement, allowing the concentration of finer meshes in regions of interest. This refinement strategy, as exemplified in situations of solution concentration for specific fiber orientations, can substantially reduce computational expenses while elevating solution accuracy.

This study presents a new numerical approach using FEM-SUPG to solve the Fokker–Planck equation, focusing on fiber orientation dynamics. While currently restricted to orientation mechanisms, the framework holds potential for future coupling simulations. Section 2 outlines the methodology for computing the evolution equation of the probability distribution function in a 2D domain, detailing the reformulation and discretization via full FEM-SUPG of the Fokker–Planck equation. Section 3 presents numerical outcomes, including validation tests and assessment of parameter influences on the probability distribution function solution. The study concludes by discussing prospects for determining the fourth-order tensor and its integration within the Stokes equation. Additionally, the presented analysis showcases the evolution of fiber orientation within the context of Poiseuille flow conditions, serving as a foundational study poised to pave the way for a comprehensive examination of the intricate fluid-flow coupling phenomenon in future work.

## 2. Mathematical modeling

### 2.1. Fokker–Planck equation

Generally, in statistical mechanics, the Fokker–Planck equation is defined as a partial differential equation that describes the temporal evolution of the probability density function of the particle’s velocity under the impact of drag forces. In short fiber reinforced thermoplastic processing, fibers can orient in different directions induced by the injection flow. Prager [30] was the first to get interested in the statistics of particle orientations by defining an orientation distribution function. In kinetic theory, a complete description of the orientation state of the composite requires calling an orientation distribution function $\varphi(\mathbf{p}, t)$ which, as discussed above, expresses the probability of a fiber for having a certain orientation $\mathbf{p}$ at time t. From Bird et al. [31], if $\mathbf{p}$ ($d\mathbf{p}/dt$) denotes the material time derivative of $\mathbf{p}$, then the conservation equation of $\varphi$ is written, neglecting the Brownian motion of the fibers as follows:

$$ \frac{\partial \varphi}{\partial t} = -\nabla_{\mathbf{p}} \cdot (\mathbf{p} \varphi) \qquad (5) $$

This equation is valid in the particular case of dilute fiber suspensions. On the other hand, in the case of semi-diluted or concentrated fiber suspensions, the hydrodynamic type interactions should be taken into account by adding in Eq. (5) a term analogous to a pseudo-Brownian diffusivity. In order to take into account fiber interactions, Folgar and Tucker [14] added to the Jeffery equation a phenomenological
<figure><img src="image_3.png" /><figcaption>Fig. 1. Illustration: time evolution of PDF with $\Delta t = 0.5$ s, $\dot{\gamma} = 1 \text{s}^{-1}$, $\lambda = 0.8$, and $C_I = 0.01$, for a simple shear flow.</figcaption></figure>

diffusion term. Consequently, the time evolution of particles can be expressed generally as [13]:

$$ \dot{\mathbf{p}}_{\text{FT}} = \underbrace{-\omega \cdot \mathbf{p} + \lambda(\mathbf{E} \cdot \mathbf{p} - \mathbf{E} : \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p})}_{\mathbf{p}_{\text{Jeffery}}} - \frac{D_r}{\varphi} \nabla \varphi \quad (6) $$

By combining Eq. (5) with Eq. (6), the well known Fokker-Planck equation is introduced as:

$$ \frac{\partial \varphi}{\partial t} = -\nabla_{\mathbf{p}} \cdot (\dot{\mathbf{p}}\varphi) + D_r \nabla_{\mathbf{p}}^2 \varphi \quad (7) $$

where $\dot{\mathbf{p}}$ is the orbital velocity and $D_r$ is the rotary diffusivity coefficient. The expression of the latter is given in [13]:

$$ D_r = C_I \dot{\gamma} \quad (8) $$

where $C_I$ is empirically determined and describes the rate of interaction. After all the Fokker-Planck Eq. (7) expresses the rate of change for the PDF $\varphi$, using the Jeffery equation for a single fiber (1) and the continuity Eq. (5). Fibers are molded as independent random variables with identical distribution and zero mean. Each interaction results in a change of orientation in both fibers. So the right hand side of Eq. (7) is the sum of two contributions: the convection and the diffusion parts. An illustration of PDF time evolution is shown in Fig. 1.

## 2.2. Reformulation of Fokker-Planck equation

### 2.2.1. Equation recasting

The Fokker-Planck equation usually presents challenges in evaluating its parameters, and the higher number of dimensions involved results in often prohibitively high computational costs. To address this issue, this section describes the methodology employed to reduce the complexity of Eq. (7), thereby minimizing the computational burden of its resolution. Notably, in this study, The Fokker-Planck equation is addressed within a 2D orientation domain, where the probability distribution function is described on a circle, and it can be solved for each single spatial node in the spatial domain (fluid domain). So at first, Eq. (7) is reformulated. Through the knowledge that divergence is a linear operator and using the product rule derived from the ordinary differentiation rules of calculus [32], the convective term is expanded and expressed in the following form:

$$ \nabla_{\mathbf{p}} \cdot (\dot{\mathbf{p}}\varphi) = \dot{\mathbf{p}} \cdot (\nabla_{\mathbf{p}}\varphi) + \varphi \nabla_{\mathbf{p}} \cdot \dot{\mathbf{p}} \quad (9) $$

Acknowledging that the orientation of a vector is defined as the angle formed by its line segment with the horizontal axis, the representation of a single fiber's orientation is achieved, as mentioned earlier, through

<figure><img src="image_4.png" /><figcaption>Fig. 2. 2D polar domain where all the possibilities of a fiber orientation exist.</figcaption></figure>

the use of the normal direction vector, denoted as $\mathbf{p}$. So $\mathbf{p}$ defined in a polar coordinates is described as function of angle ($\theta$ in 2D) and radius ($r = 1$). Consider a fiber initially oriented at an angle $\theta_0$ and gradually rotated to an angle $\theta_1$. Fig. 2 illustrates the fiber orientations at $\theta = \theta_0$ and $\theta = \theta_1$ denoted by $\mathbf{p}_0$ and $\mathbf{p}_1$, respectively. By determining the slope (gradient) of $\mathbf{p}$ using trigonometric expressions related to the polar angle, the change in $\mathbf{p}$ direction representing the fiber orientation corresponds to the change in its angle of rotation. Thus, an association between $\mathbf{p}$ and $\theta$ is established as $\nabla_{\mathbf{p}} = \nabla_{\theta}$. Consequently, in the 2D case, $\varphi$ can also be expressed as a function of the orientation angle $\theta$. Combining this expression with Eq. (9), the reformulated Fokker-Planck equation is obtained:

$$ \frac{\partial \varphi}{\partial t} + \dot{\mathbf{p}} \cdot (\nabla_{\theta} \varphi) + \varphi \cdot \nabla_{\theta} \dot{\mathbf{p}} - D_r \nabla_{\theta} \cdot (\nabla_{\theta} \varphi) = 0 \quad (10) $$

### 2.2.2. Angular velocity calculation

The orbital velocity $\dot{\mathbf{p}}$ contributes to fiber motion, both in terms of rotation and deformation. Thus, the first term in Eq. (10) leads to pure rotation of the fibers, the second term generates deformation(elongation) which causes also a rotation and the third term contributes to only elongation. This can be mathematically proven by simple calculations of each term of $\mathbf{p}$. Taking for example a simple 2D shear flow in the x direction, its velocity profile can be given by:

$$ \mathbf{u} = \begin{bmatrix} u_x \\ u_y \end{bmatrix} = \dot{\gamma} \cdot \begin{bmatrix} y \\ 0 \end{bmatrix} \quad (11) $$

where $\mathbf{u}$ has one non-zero component, $u_x$, whose gradient is perpendicular to the velocity itself and equal to $\dot{\gamma}$. For that case, this implies that Eq. (10) can be solved separately for every single spatial node as the velocity gradient is constant over the domain, which means that all fibers orient themselves with the same rotation rate. Firstly, vorticity and deformation tensors are computed from expressions (2), knowing that $(\nabla u)_{ij}$ is expressed as $\partial u_i / \partial x_j$ in index notation (with $i, j \in [x, y]$), and then plugged into Jeffery's part of Eq. (6) to calculate each term of $\dot{\mathbf{p}}$, respectively. The components of the vector $\mathbf{p}$ are expressed in polar coordinates and the expressions of the first two terms in Eq. (6) can be determined:

$$ \begin{align} \mathbf{E} \cdot \mathbf{p} &= \frac{1}{2} \begin{bmatrix} 0 & \dot{\gamma} \\ \dot{\gamma} & 0 \end{bmatrix} \cdot \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} = \frac{1}{2} \dot{\gamma} \begin{bmatrix} \sin \theta \\ \cos \theta \end{bmatrix} \tag{12} \\ \omega \cdot \mathbf{p} &= \frac{1}{2} \begin{bmatrix} 0 & -\dot{\gamma} \\ \dot{\gamma} & 0 \end{bmatrix} \cdot \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} = \frac{1}{2} \dot{\gamma} \begin{bmatrix} -\sin \theta \\ \cos \theta \end{bmatrix} \end{align} $$

The third term is re-written in a simpler manner to facilitate its calculation, without necessarily computing the third order tensor resulting from the tensor product between the three unit vectors $\mathbf{p}$:

$$ \mathbf{E} : \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} = (\mathbf{p}^T \cdot \mathbf{E} \cdot \mathbf{p}) \cdot \mathbf{p} \quad (13) $$

where $\mathbf{p}^T$ is the transposed vector of $\mathbf{p}$. Then, similarly:

$$ (\mathbf{p}^T \cdot \mathbf{E} \cdot \mathbf{p}) \cdot \mathbf{p} = \dot{\gamma} (\cos \theta \sin \theta) \cdot \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} \quad (14) $$

A vector field is defined as a function that maps points to vectors, thus the fields in Eqs. (12)–(14) which describe the physical motions of a
<figure><img src="image_5.png" /><figcaption>Fig. 3. p-effect applied to fiber motion: (a) pure rotation (Eq. (12)); (b) elongation+rotation (Eq. (12)); (c) pure elongation (Eq. (13)).</figcaption></figure>

<table><caption>Table 1
Representation of vector fields.</caption><thead><tr><th>Fields</th><th>θ</th><th>0</th><th>π/4</th><th>π/2</th><th>π</th></tr></thead><tbody><tr><td>E.p</td><td></td><td>(0, β)</td><td>(β√2/2, β√2/2)</td><td>(β, 0)</td><td>(0, -β)</td></tr><tr><td>ω.p</td><td></td><td>(0, β)</td><td>(-β√2/2, β√2/2)</td><td>(-γ, 0)</td><td>(0, -β)</td></tr><tr><td>(p<sup>T</sup>.E.p).p</td><td></td><td>0</td><td>(γ√2/2, γ√2/2)</td><td>0</td><td>0</td></tr></tbody></table>

fiber, can be sketched by plugging some coordinates of different points as shown in the following table (see Table 1), where $β = γ/2$ and $γ = γ(\cos θ \sin θ)$. Starting from the configuration presented in Fig. 2, contributions of the three terms on the vector fields are drawn in Fig. 3.

Furthermore, the vector field corresponding to Eq. (14) and the unit vector **p** are shown to be co-linear, physically meaning that the fiber experiences only deformation considering the effect of Eq. (14). This implies that this term can be negligible if fibers are very stiff. Hence, equation Eq. (8) is recast as follows:

$$ \dot{\mathbf{p}} = -\omega \cdot \mathbf{p} + \lambda(\mathbf{E} \cdot \mathbf{p}) \quad (15) $$

With the help of Eqs. (2), (15) is written in the index notation where $E_{ij}$ and $ω_{ij}$ are the components of the strain rate and the vorticity tensors respectively. This results in:

$$ \dot{\mathbf{p}} = \begin{bmatrix} \cos \theta (\lambda E_{xx} - \omega_{xx}) + \sin \theta (\lambda E_{xy} - \omega_{xy}) \\ \cos \theta (\lambda E_{yx} - \omega_{yx}) + \sin \theta (\lambda E_{yy} - \omega_{yy}) \end{bmatrix} \quad (16) $$

where $(ω_{xx} = ω_{yy} = 0)$. As aforementioned, the rotational effect of **p** is only sought in order to predict the fiber orientation. Thus, the angular velocity $\dot{θ}$, which is the rate of change of angle in time, is determined in order to compute the fiber rotation. As seen in Fig. 4, **p**

<figure><img src="image_6.png" /><figcaption>Fig. 4. Projection of **p**.</figcaption></figure>

can be decomposed into radial and tangential orientation velocities by projecting it on unit vectors $V_r$ and $V_t$ respectively. So, the projection of $\dot{θ}$ is determined from the following formula:

$$ \text{proj}_{V_t} \dot{\mathbf{p}} = \left( \frac{\dot{\mathbf{p}} \cdot V_t}{\|V_t\|^2} \right) \cdot V_t \quad (17) $$

where $(\|V_t\|^2 = 1)$ in case of the unit circle. Eq. (16) is plugged into expression (17), **p** is projected and $\dot{θ}$ is computed as follows:

$$ \dot{\theta} = \frac{1}{2} \cos^2 \theta (\lambda E_{xy} - \omega_{yx}) - \frac{1}{2} \sin^2 \theta (\lambda E_{xy} - \omega_{xy}) + \lambda \cos \theta \sin \theta (E_{yy} - E_{xx}) \quad (18) $$

The final form of the reformulated Fokker-Planck equation, which will be later discretized in the $θ$-space, is:

$$ \frac{\partial \varphi}{\partial t} + \varphi \nabla_{\theta} \dot{\theta} + \dot{\theta} (\nabla_{\theta} \varphi) - D_r \nabla_{\theta} \cdot (\nabla_{\theta} \varphi) = 0 \quad (19) $$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Projecting 2D to 1D discretization domain.</figcaption></figure>

<table><caption>**Table 2**
Nodal connectivity matrix.</caption><thead><tr><th>Element</th><th>Topology</th></tr></thead><tbody><tr><td>1</td><td>1,2</td></tr><tr><td>2</td><td>2,3</td></tr><tr><td>.</td><td>.</td></tr><tr><td>.</td><td>.</td></tr><tr><td>.</td><td>.</td></tr><tr><td>N</td><td>N-1,N</td></tr><tr><td>N</td><td>N 1</td></tr></tbody></table>

## 2.3. Discretization of Fokker-Planck equation

### 2.3.1. 2D-1D configuration

Initially, the configurational domain is represented by a unit circle in the 2D plane, to reduce the complexity and the computational effort. The 2D domain is converted into a 1D domain by projecting the unit circle into a 1D segment whose length is $2\pi$, as shown in Fig. 5. As seen in Table 2, the numbers of nodes and elements are equal, and the last node is connected to the first node.

### 2.3.2. Finite element approach

As mentioned previously, there are a large number of studies in the literature concerning the numerical calculation of fiber orientation by solving the Fokker-Planck equation (FPE). On the contrary, few of them were able to get an accurate solution while maintaining its stabilization. The discretization of the geometries plays a key role for an accurate resolution of the FPE. The initial step is the definition of a proper computational domain (see Fig. 5), representing as good as possible the characteristics and shape of the simulated problem. So, the FEM was used in this work to compute the FPE which allows theoretical studies such as uniqueness, existence, consistency, and stability of solution to be performed, due to its strong mathematical developments. Also it assists in taking advantages from techniques like parallelization, a technique that use full capacity of the modern day computers for obtaining quick results, which have also been developed for enhancing in particular the FEM discretization technique; the parallelization technique will be advantageous especially for future works related to the flow-fiber coupling problem. So in order to solve Eq. (19), the SUPG approach [33] instead of the standard Galerkin method, as the latter either produces excessive numerical diffusion, or non-physical oscillations as a first order differential operator also appears for FPE. The SUPG consists in disturbing the test functions for the convection terms by adding a term which is similar to the artificial diffusion:

$$ \tilde{\phi}_{supg} = \tilde{\phi} + \gamma \dot{\theta} \cdot \nabla \tilde{\phi} \quad (20) $$

where $\gamma$ is a stabilizing coefficient that controls the artificial diffusion, expressed as:

$$ \gamma = \frac{h}{(2|\dot{\theta}| + 2 \frac{D_{\kappa}}{h})} \quad (21) $$

where $h$ is the mesh size. To begin with, the various spaces needed to express the weak formulation of FPE's physical problem, are introduced. The spaces used in ICI-tech (solver libraries from Centrale Nantes, France) are recalled:

* $C^0(\Omega_\theta)$: the space of functions continuous on $\Omega_\theta$

* $P^1(\Omega_\theta)$: the space of polynomials of degree 1 on $\Omega_\theta$

* $L^2(\Omega_\theta)$: the Lebesgue space of square-integrable functions on $\Omega_\theta$

$$ L^2(\Omega_\theta) = \{w \in L^2(\Omega_\theta), D w \in L^2(\Omega_\theta)\} \quad (22) $$

* $H^1(\Omega_\theta)$: the Sobolev space of functions whose order 1 derivatives belong to $L^2(\Omega_\theta)$

$$ H^1(\Omega_\theta) = \{w \in L^2(\Omega_\theta), D w \in L^2(\Omega_\theta)\} \quad (23) $$

Thus, the spaces for the PDF and its arbitrary function are defined as:

$$ \begin{align} V &= \{\varphi \in H^1(\Omega_\theta)\} \tag{24} \\ W &= \{\tilde{\varphi} \in H^1(\Omega_\theta), \dot{\theta} \cdot \nabla \tilde{\varphi} \in H^1(\Omega_\theta)\} \nonumber \end{align} $$

Galerkin's variational formulation of the problem is then constructed by multiplying Eq. (19) by the test function $\tilde{\varphi}(x) \in W$ and integrating over the entire domain $\Omega_\theta$, results in:

$$ \begin{align} & (\partial_t \varphi, \tilde{\varphi}) + (\gamma \dot{\theta} \cdot \partial_t \varphi, \nabla_\theta \tilde{\varphi}) + (\nabla_\theta \dot{\theta} \cdot \varphi, \tilde{\varphi}) + (\dot{\theta} \cdot \nabla_\theta \varphi, \tilde{\varphi}) \nonumber \\ & \qquad + (\dot{\theta} \cdot \nabla_\theta \varphi, \gamma \dot{\theta} \cdot \nabla_\theta \tilde{\varphi}) + D_r(\nabla_\theta \varphi, \nabla_\theta \tilde{\varphi}) = 0 \tag{25} \end{align} $$

where the notation $(\cdot, \cdot)$ defines the scalar product in the sense of the $L^2(\Omega_\theta)$:

$$ \forall (u, v) \in L^2(\Omega_\theta), (u, v) = \int_{\Omega_\theta} uv d\Omega_\theta \quad (26) $$

To solve this problem, the domain is discretized and broken down into simplexes $\mathcal{K}$ formed by 1D-elements (see Fig. 5). Thus, the preceding spaces of continuous functions are approximated on $\mathcal{K}$ by discrete subspaces of finite dimensions constituted by piece-wise continuous linear functions.

$$ \begin{align} \mathcal{V}_h &= \{\varphi_h | \varphi_h \in C^0(\Omega_\theta), \forall \mathcal{K}, \varphi_{h|\mathcal{K}} \in P^1(\mathcal{K})\} \tag{27} \\ \mathcal{W}_h &= \{\tilde{\varphi}_h | \tilde{\varphi}_h \in \mathcal{V}_h\} \nonumber \end{align} $$

The full-discrete Galerkin problem consists of finding $\varphi_h \in \mathcal{V}_h$ such as $\tilde{\varphi}_h \in \mathcal{W}_h$ for the following problem:

$$ \begin{align} & (\partial_t \varphi_h, \tilde{\varphi}_h) + (\gamma \dot{\theta} \cdot \partial_t \varphi_h, \nabla_\theta \tilde{\varphi}_h) + (\nabla_\theta \dot{\theta} \cdot \varphi_h, \tilde{\varphi}_h) + (\dot{\theta} \cdot \nabla_\theta \varphi_h, \tilde{\varphi}_h) \nonumber \\ & \qquad + (\dot{\theta} \cdot \nabla_\theta \varphi_h, \gamma \dot{\theta} \cdot \nabla_\theta \tilde{\varphi}_h) + D_r(\nabla_\theta \varphi_h, \nabla_\theta \tilde{\varphi}_h) = 0 \tag{28} \end{align} $$

The approximation of the solution on an element is determined using the values at the nodes. Then the solution is approached by:

$$ \begin{align} \varphi_{h|\mathcal{K}} &= \sum_{i=1}^{D_{\mathcal{K}}} \varphi_i(t) N^i(x) \tag{29} \\ \tilde{\varphi}_{h|\mathcal{K}} &= \sum_{j=1}^{D_{\mathcal{K}}} \tilde{\varphi}_j(t) N^j(x) \nonumber \end{align} $$
where $i$ represents the node, $N(x)$ the shape function on that node and $D_K$ the number of nodes in the $K$-simplex. The existence of a time derivative in the FPE requires a time discretization. The FPE problem is solved at different instants of the time, separated by a time step $\Delta t$. The choice of $\Delta t$ at each instant needs to be well thought to capture the time rate of change of fiber orientation required by both the variational formulation and $\theta$-space discretization. Thus, three different temporal discretization schemes (backward Euler, Crank–Nicolson and Rosenbrock) were implemented and compared with each other (see Section 3). Eventually Eq. (29) is introduced into Eq. (28), which generates a linear system of algebraic equations.

So once the stabilized weak formulation is well written, the equivalent matrix formulation is deduced. In ICI-tech, the method adopted for the resolution of the system is an iterative method of conjugate gradient type resolution. To do this, the library PETSc is used. PETSc [34] is a library of algorithms for solving linear and non-linear systems, both sequentially and in parallel, using iterative numerical methods. In the context of solving partial differential equations such as the ones in the article, PETSc offers interesting methods including preconditioners. Local arrays are created in ICI-tech as local solvers. Subsequently, PETSc supports them for preconditioning. A preconditioner, using an incomplete factorization method called the decomposition ILU($k$) of a matrix $A$ in the form $A=LU$, is used. The parameter $k$ makes it possible to control the filling ratio of the hollow matrix (sparse matrix), thus making it possible to solve even problems of a poorly conditioned nature. It should be noted that the resolution time with such an iterative process varies according to the number of iterations performed at the end of the resolution. C++11 is used as a programming language to write the FEM code of solving the probability distribution function.

# 3. Numerical results: solving the Fokker-Planck equation for fiber orientation

## 3.1. Validation tests

In this section, we conduct a validation study using a simple shear flow case characterized by a constant generalized shear rate ($\dot{\gamma} = 1s^{-1}$). The primary objectives of this validation are to confirm the accuracy of the computed angular velocity obtained from Eq. (18) and to verify the reliability of the numerical solution of the Fokker-Planck equation. This validation includes assessing the normalization condition and investigating the influence of the shape factor $\lambda$ and the interaction coefficient $C_I$ on the numerical results.

### 3.1.1. Normalization condition

The 2D orientation distribution function is introduced, as described above, to represent the probability of finding a fiber between the angles ($\theta$ and $\theta + d\theta$), given as:

$$P = \int_{0}^{\theta+d\theta} \varphi(\theta)d\theta \qquad (30)$$

Moreover, the integral over all possible orientations of the distribution function must be equal to unity satisfying the normalization condition induced by the definition of a probability:

$$P = \int_{0}^{2\pi} \varphi(\theta)d\theta = 1 \qquad (31)$$

The latter expression was verified, as illustrated in Fig. 6. When solving a pure convection Fokker-Planck equation ($D_r = 0$) across the entire domain $\Omega_\theta$, the integral of the probability distribution function exhibits remarkable constancy, remaining close to 1 throughout varying time $t$ for different shape factors $\lambda$. When examining Fig. 6, it becomes evident that the error is extremely small, fluctuating within the order of $10^{-11}$ for all different shape factors $\lambda$. Notably, a distinct pattern is observed in the error plot for the case of $\lambda = 1$, where significant fluctuations occur compared to other shape factors. This behavior

arises due to the physical phenomenon of fiber concentration or high gradient in the resolution of the Fokker-Planck equation (see Fig. 8). In this particular case, the fibers tend to align at a certain moment, leading to a concentration of the solution at a specific orientation angle. Nonetheless, despite the fluctuations, the error magnitude remains remarkably small and within the order of $10^{-11}$ across all shape factor values, as this observation also highlights the effectiveness of the SUPG finite element method in stabilizing critical cases involving convection-dominated problems, thereby ensuring the preservation of the normalization condition.

### 3.1.2. Influence of $\lambda$-parameter

Fig. 8 shows the variation of the angular velocity for different values of $\lambda$ parameter. When $\lambda = 0$, the fiber is represented by a circular shape (Fig. 7), and it rotates with an uniform circular motion as the angular velocity remains constant at any angle $\theta$. The fibers will keep subsequently rotating over themselves with a steady velocity if fibers do not intersect. As for $\lambda = 1$, the case where the fibers are assumed quasi-infinite, the angular velocity is zero at 0 and $\pi$. Therefore, the fibers that are oriented in the direction of the velocity field will not rotate. This is logical as the velocity gradient along the y-axis does not affect the rotation of fibers as their thickness tends to be zero. Therefore, we conclude that fibers initially oriented at a small angle with respect to the flow, will take longer to rotate in the case of shear flow if $\lambda$ value is close to 1. In other words, it could be verified that the angular velocity values will be the highest for $\lambda = 1$ when $\theta \in ]\pi/4,3\pi/4[ \cup ]5\pi/4,7\pi/4[$. It is also noticeable that the maximum angular velocity is always at $\pi/2$ and $3\pi/2$ where the torque exerted by the flow on a fiber is at its maximum and the latter state is perpendicular to the velocity streamlines.

The Fokker-Planck equation has been computed in order to determine the probability distribution function that evolves over time to understand the change of fiber orientation with respect to time. The Fokker-Planck equation is an advection-diffusion equation, where the advection term represents the deterministic dynamics of the fibers, and the diffusion term represents the stochastic effects due to the interaction between the fibers. To evaluate the Fokker-Planck equation, we initialize a probability distribution function at time $t = 0$ that characterizes a scenario wherein the fibers are aligned as depicted in Fig. 9. This alignment scenario is characterized by fibers being oriented around angles 0 and $\pi$, where the specific angular ranges are defined as follows: $\pi - c\Delta\theta$ to $\pi + c\Delta\theta$ for fibers oriented around $\pi$, and range that spans from 0 to $0 + c\Delta\theta$ and from $2\pi - \Delta\theta(1 + c)$ to $2\pi - \Delta\theta$ for fibers around 0. This range accounts for the periodicity in the domain at angle 0, where $2\pi - \Delta\theta$ returns to 0, signifying that the angle $2\pi$ is equivalent to the angle 0. Thus, the complete representation of the total probability for this alignment scenario is defined by the integral equation:

$$P_{\text{total\_alignment}} = P(\theta \le 0 + c\Delta\theta) + P(\theta \ge 2\pi - \Delta\theta(1 + c)) \\ + P(\pi + c\Delta\theta \ge \theta \le \pi + c\Delta\theta) \qquad (32)$$

In this expression, $P_{\text{total\_alignment}}$ gives us the representation of the total alignment within $+/-c$. Then, the Fokker-Planck equation is resolved, and from its solution of the probability distribution function, the probability of finding fibers oriented around angles 0 and $\pi$ is computed at each time step. The range is deliberately kept small with ($c = 32$ and $\Delta\theta = \frac{2\pi}{2048}$) in order to effectively capture the variation of the probability of finding aligned fibers over time. The shape factor of fibers plays a significant role in the evolution of the probability distribution function over time as it is shown in Fig. 10. Starting with an initial distribution function that depicts a unidirectional fiber orientation case, where most of the fibers (about 90%) are oriented with angles spanning both 0 and $\pi$, the probability distribution function decreases proportional to different values of shape factor of fibers over time. Physically, the fibers tend to orient faster around 0 and $\pi$ as the shape factor gets smaller ($\lambda \to 0$) representing a more circular
<figure><img src="image_8.png" /><figcaption>Fig. 6. Time evolution and error analysis of φ integral for different λ: ensuring normalization.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7. λ-effect on fibers shape.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Variation of $\dot{\theta}$ versus $\theta$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Initialization of probability distribution function at $t = 0$ with $c = 32$ and $\Delta\theta = \frac{2\pi}{2048}$.</figcaption></figure>

shape of fibers, compared to fibers with a shape factor approaching 1, which is consistent with previous results (see Fig. 8). This is due to the fact that with shape factors ($\lambda \to 1$) represent more elongated fibers, they have a greater tendency to remain aligned and oriented in the direction of the flow, in this case in both 0-direction and $\pi$-direction, for a longer period of time. Conversely, fibers with small shape factors tend to orient slower than fibers with quasi-infinite shape factors when far from angles 0 and $\pi$, which leads to a slower alignment with the initial angles 0 and $\pi$. This interpretation holds true in the absence of the diffusion term. However, when the diffusion term is included, depicted in Fig. 10, fibers with shape factors closer to zero undergo subsequent rotation until the probability distribution function becomes more uniform, indicating a more isotropic orientation. In contrast,

fibers with a shape factor ($\lambda = 1$) tend to remain aligned with angles 0 and $\pi$, serving as a counterforce to the diffusion term, even though a decrease of the peak of the probability distribution function near angles 0 and $\pi$ occurred, but still resulting in a higher probability of finding fibers oriented in both angles 0 and $\pi$ compared to the case of shape factors less than 1. These results demonstrate the important role of the shape factor in the evolution of fiber orientation over time.

### 3.1.3. Influence of $C_1$-parameter

An analysis is also done on the effect of the diffusion term on the fiber orientation to understand how it affects the probability distribution function. Thus, Fig. 11 represents the evolution of the probability
<figure><img src="image_12.png" /><figcaption>Fig. 10. Time evolution of the total probability of finding fibers oriented around angles 0 and π for different λ.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Time evolution of the total probability of finding fibers oriented around angles 0 and π for different C$_{l}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Variation of the probability distribution function φ at time t = 10 s for different C$_{l}$.</figcaption></figure>

distribution centered around 0 and π with some degree of spread cΔθ with (c = 32 and Δθ = 2π/2048) by solving the full Fokker-Planck equation (advection-diffusion equation), with different interaction coefficients

<figure><img src="image_15.png" /><figcaption>Fig. 13. Fibers orientation for the Poiseuille flow case.</figcaption></figure>

C$_{l}$ to investigate their effects on the solution of the Fokker-Planck equation. The initial probability distribution function represents a scenario where a large majority of fibers (approximately 90%) are oriented around angles of 0 and π. The diffusion term arises from the Brownian motion of fibers and causes the distribution to spread out over time, resulting in a decrease in the peak of the distribution. This behavior can be observed by comparing the solutions of the Fokker-Planck equation with and without the diffusion term. Fig. 11 demonstrates the effect of the interaction coefficient on the evolution of the probability of finding fibers oriented around angles 0 and π. When the C$_{l}$ is set to zero, the probability exhibits a periodic behavior where it returns to its initial value after a fixed period of time. However, when the C$_{l}$ is increased, the probability evolves in a way that results in a decrease in the peak value before flattening out as it is illustrated in Fig. 12. This is because diffusion acts to smooth out fluctuations in the fiber orientation and promotes a more uniform distribution of orientations indicates an isotropic distribution. As a result, the probability distribution becomes less peaked and more scattered, indicating that the fibers will be oriented in a wider range of angles meaning that there is no preferred orientation and the orientation of the fibers is considered to be random. This effect is particularly evident in the long-term behavior of the system, where the probability distribution tends to approach a steady state that is less sensitive to the initial conditions. Thus, knowing that the diffusion term in the Fokker-Planck equation captures the stochastic character of the system and describes how the probability density evolves over time due to random fluctuations, the diffusion term plays a role in the stabilization phenomenon, by allowing the system to explore different orientations of the fiber and eventually converge to a stationary distribution, as the diffusion term represents the stochastic effects of the fiber-fiber interactions.

## 3.2. Numerical convergence analysis and computational efficiency

In this section, we showcase numerical results concerning the resolution of the Fokker-Planck equation within a Poiseuille flow scenario, depicted in Fig. 13. This scenario is characterized by a position-dependent generalized shear rate, denoted as $\dot{\gamma}$, which varies across the position y. Its specific form is expressed as:

$$ \dot{\gamma} = \frac{8y}{H^2} V_{max} \qquad (33) $$

where H is the distance between the two horizontal parallel plates, $V_{max}$ is the maximum axial velocity. Throughout this section, the parameter values remain consistent, with H set to 1 m, $V_{max}$ set to 10 m/s, while the shape factor λ is assigned a value of 0.98, and the interaction coefficient C$_{l}$ is held at 0. The scope of these numerical results encompasses a thorough convergence analysis, focusing on both mesh refinement and time discretization schemes. Furthermore, a comprehensive comparison is drawn to discern the distinct impact of these schemes on computational efficiency.
### 3.2.1. Convergence study: mesh and time discretization schemes

To conduct the numerical convergence study, the computation of the second-order orientation tensor $\mathbf{a}_2$ was undertaken. This tensor is derived through the resolution of the Fokker-Planck equation, where the probability distribution function $\varphi$ satisfies the requirements of periodicity, normalization, and continuity. The computation of $\mathbf{a}_2$ involves integrating dyadic products with the probability distribution function $\varphi$, capturing the spatial average of the double tensor product of $\mathbf{p}$ [35], expressed as follows:

$$ \mathbf{a}_2 = \int_{\Omega_0} \varphi(\theta) \mathbf{p} \otimes \mathbf{p} \, d\mathbf{p} \qquad (34) $$

Hence, through the calculation of $\varphi(\theta)$ using Eq. (19), the planar orientation of the second-order tensor can be deduced using Eq. (34). This approach unveils the fiber orientation state within the primary flow direction, eliminating the necessity for employing closure approximations that can undermine accuracy. Consequently, the components of the $\mathbf{a}_2$ tensor can be expressed as:

$$ \begin{aligned} a_{xx} &= \sum_{i=1}^{N} \varphi_i \cos^2(\theta_i) \\ a_{yy} &= \sum_{i=1}^{N} \varphi_i \sin^2(\theta_i) \end{aligned} \qquad (35) $$

In pursuit of a comprehensive numerical convergence assessment, depicted in Fig. 14, the time evolution of the orientation tensor $\mathbf{a}_2$ is computed at $y = 0.5$, corresponding to a generalized shear rate $\dot{\gamma} = 5 \text{ s}^{-1}$ based on Eq. (33). The initial condition is set to represent an isotropic state, with $a_{xx} = a_{yy} = 0.5$, and the period of its evolution is $T = 4\pi$s. To validate the results, a comparative analysis is performed, involving a comparison between the obtained numerical results and the analytical solution derived in [36].

A comparative study, depicted in Fig. 14(a), is initially conducted using various temporal discretization schemes, seeking to assess the sensitivity of the obtained numerical solution to various time discretization approaches. This study reveals that the employment of the backward Euler scheme “beuler”, recognized as a first-order technique, results in significant loss of information and introduces inaccuracies that amplify over time, compared to the analytical solution. In contrast, the Crank-Nicolson scheme “cn” and the Rosenbrock-W scheme “rosw” exhibit improved performance. The observed discrepancies can be attributed to the inherent limitations of the backward Euler method in accurately capturing swift variations or oscillatory patterns. The Crank-Nicolson and Rosenbrock-W methods, distinguished by their respective local truncation errors of second and third order [37,38], exhibit similar solution trends, apparent from the overlapping curves. Despite this, they yield better solutions with less damping, with their curves approaching the analytical solution and displaying smaller errors.

Fig. 14(b) examines the influence of altering $\Delta t$ on the time evolution solution of $a_{xx}$. Employing the Crank-Nicolson time discretization scheme, the time step order is adjusted across the range of $(2 \times 10^{-2} \text{ to } 2 \times 10^{-4} \text{ s})$. Interestingly, decreasing the time step does not contribute to the convergence or enhancement of solution accuracy. This observation may be attributed to the circumstance that in situations where dominant effects are associated to convection or rapid changes, reducing the time step might not result in significant improvements, particularly if the underlying numerical scheme already provides adequate accuracy.

Lastly, an investigation into mesh convergence is conducted, as illustrated in Fig. 14(c), involving a range of mesh sizes from $(\frac{2\pi}{512} \text{ to } \frac{2\pi}{2048})$ and using a Crank-Nicolson time scheme. The objective is to address damping effects in the results. The plot demonstrates that as mesh refinement is undertaken, the numerical results progressively align with the analytical solution. Notably, for the case of a mesh size of $h = \frac{2\pi}{2048}$, the curves of the analytical and numerical results coincide, signifying the importance of mesh refinement to mitigate damping. This necessity arises due to the solution of the Fokker-Planck equation, which encompasses a probability distribution function

<table><caption>**Table 3**
Computational cost analysis.</caption><thead><tr><th>Scheme</th><th>Beuler</th><th>cn</th><th>Rosw</th></tr></thead><tbody><tr><td>h</td><td></td><td></td><td></td></tr><tr><td>512</td><td>0.16</td><td>0.172</td><td>0.365</td></tr><tr><td>1024</td><td>0.24</td><td>0.265</td><td>0.524</td></tr><tr><td>2048</td><td>0.423</td><td>0.469</td><td>1</td></tr></tbody></table>

with sharp gradients occurring at different time points around angles $\theta = 0$ and $\theta = \pi$. Thus, neglecting mesh refinement can lead to compromised solution accuracy. So to demonstrate the origins of the damping phenomenon, Fig. 15 captures the evolution of $\varphi$ at a time equals to half the period ($t = T/2$), corresponding to the state of $a_{xx}$ at its minimum value ($a_{xx} = 0.5$), representing randomly oriented fibers. Notably, as the mesh size is reduced, the oscillations in the probability distribution function become more stabilized, approaching a constant value of $1/(2\pi)$ with refined meshes. This trend aligns with the isotropic orientation distribution state, thereby converging toward the analytical solution.

### 3.2.2. Computational efficiency: impact of refinement and higher order temporal discretization schemes

A comprehensive computational cost analysis has been carried out, covering a comparison of various mesh sizes $h$ and distinct time discretization schemes. Notably, the examination of different time steps has been excluded, as evidenced by the findings in Fig. 14(b), where it was established that altering the time step does not significantly affect the accuracy of Fokker-Planck equation resolution. As a result, the focus of the time computation comparative study will remain on a fixed time step of $\Delta t = 2 \times 10^{-2}$, ensuring a thorough exploration of the impact of mesh refinement and time discretization schemes on computational efficiency. The comparative study will be based on the resolution of the Fokker-Planck equation for half period $t_{end} = T/2$.

It is worth noting that the computation time for solving the Fokker-Planck equation on one spatial node (fluid domain) is relatively small. To better understand the effects of mesh refinement and time schemes, the comparative study will be conducted relative to the last case, where the time discretization scheme is considered Rosenbrock-W and the mesh size $h$ is $\frac{2\pi}{2048}$, representing the largest computational time ($t_{computed} = 0.915 \text{ s}$). This approach provides a comprehensive overview of the computational time impact, serving as a baseline for future work related to coupling problems, where the Fokker-Planck equation will be computed for all spatial nodes in the fluid domain.

Table 3 presents first the relative computational time variations resulting from transitioning from the backward Euler scheme (beuler) to the Crank-Nicolson scheme (cn) and the Rosenbrock-W scheme (rosw). Interestingly, shifting to the Crank-Nicolson scheme leads to a much smaller increase in computational time than shifting to the Rosenbrock-W scheme. This is observed across all three mesh sizes. On the other hand, with the refinement of the mesh via a decrease in mesh size $h$, a substantial increase in relative computational time is evident, scaling by a factor of approximately 1.7 for all cases of temporal schemes. The combined insights from Fig. 14 and the data within Table 3 indicate that opting for higher accuracy in solution while maintaining computational efficiency is best achieved by employing a mesh size of $h = \frac{2\pi}{2048}$, coupled with the Crank-Nicolson time discretization scheme. Notably, the Crank-Nicolson scheme demonstrates enhanced accuracy results with nearly equivalent computational time compared to the backward Euler scheme, and it significantly outperforms the Rosenbrock-W scheme, by consuming almost half the computational time for providing the same accuracy. Furthermore, for $h = \frac{2\pi}{2048}$, despite its relatively higher computation time, it offers notably better results compared to the other two mesh sizes. It is noteworthy that the computational time discrepancy can be managed through the implementation of mesh adaptive methods, which strategically refine the mesh in specific regions of interest, thereby effectively mitigating the increased computational time associated with finer mesh sizes.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Numerical convergence analysis: second order orientation evolution with varied temporal schemes (backward Euler, beuler, Rosenbrock-W, rosuv, and Crank–Nicolson, cn), and time steps and mesh sizes (for Crank–Nicolson time scheme).</figcaption></figure>

### 3.3. Fourth order tensor computation

The anticipation of fiber orientation in industrial contexts, as previously elucidated, relies upon a constructed model employing a tensor framework. This entails the resolution of the Folgar–Tucker equation to calculate the temporal evolution of $a_2$ [8,39–42]. However, the computational advantage of directly addressing the Fokker–Planck equation lies in its capacity to compute the higher-order orientation tensor $a_4$ without necessitating recourse to closure approximations. Thus, $a_4$ is considered to be a better choice than a second-order tensor $a_2$, as it provides more detailed information on the distribution of fiber orientations within the material, which can be used to well understand the changes in fiber orientation over time and its influence on the flow during the injection process. As a result, using the fourth-order tensor is likely to be more effective in order to resolve the complex coupling problem in future work.

In a manner analogous to the computation of the second order orientation tensor $a_2$ in (Section 3.2.1), we proceeded to calculate the fourth order orientation tensor $a_4$, defined as follows:

$$a_4 = \int_{\Omega_0} \varphi(\mathbf{p}, t) \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \otimes \mathbf{p} \, d\mathbf{p} \qquad (36)$$

To validate the numerical results, an analytical solution from [36] was employed for comparison. The comparative analysis was carried out under the conditions of a simple shear flow with a generalized shear rate $\dot{\gamma} = 1 \text{ s}^{-1}$, an interaction coefficient $C_I = 0$, and two different shape factor values: $\lambda = 1$ and $\lambda = 0.98$. The initial orientation was set to an isotropic state. For the numerical Fokker–Planck resolution, following the insights obtained in (Section 3.2.2), the Crank–Nicolson

scheme was used for time discretization, with a time step $\Delta t$ of $2 \times 10^{-2}$ s and a mesh size $h$ of $\frac{2\pi}{1024}$.

The results, depicted in Fig. 16 for the $\lambda = 1$ case demonstrate a high level of accuracy, evidenced by the nearly seamless alignment between the analytical and numerical solutions for all plotted components of $a_4$. As for the results presented in Fig. 17, they highlight the behavior of the numerical solutions for the $\lambda = 0.98$ case. A notably small error is observed when comparing it with the analytical solution, particularly within the periodicity zone of the time evolution of $a_4$, where each period of evolution is completed. As previously discussed in Section 3.2.1, this error is attributed to the damping phenomenon, which can be mitigated through mesh refinement. Contrary, this damping phenomenon is absent in the $\lambda = 1$ case, where the time evolution of $a_4$ exhibits non-periodic behavior, and the fibers tend to align over time.

Thus, in both cases of $\lambda = 0.98$ and $\lambda = 1$, the computed results exhibited significantly improved accuracy with minimal error when compared to the results obtained from solving the Folgar and Tucker equation using closure approximations as demonstrated in [36]. This validates the effectiveness of the direct solver approach using the finite element method FEM-SUPG to solve the Fokker–Planck equation for determining the evolution of the fiber fourth-order orientation tensor. The results also provide confidence in the ability of the approach to accurately capture the complex coupling between fiber orientation and flow in the injection process.

### 4. Flow-induced orientation

To further improve and move on to simulating and predicting fiber orientation induced by the flow, the Fokker–Planck equation resolution
<figure><img src="image_17.png" /><figcaption>Fig. 15. Evolution of $\phi$ at isotropic orientation state ($\frac{1}{T} = 0.5$) with varying mesh refinement.</figcaption></figure>

has been coupled to fluid dynamics. The preliminary results obtained presented in this section will be a precursor for further advanced and detailed investigations in simulating fiber behavior during the injection process.

In the context of fiber-flow coupling, accurate prediction of system behavior requires the consideration of multiple length scales. On the macroscopic scale, fluid flow of is described using conservation of momentum, under the hypothesis of incompressibility and absence of inertia effects:

$$ 
\begin{aligned}
\nabla \cdot \mathbf{u} &= 0 \quad \text{in } \Omega \\
\nabla \cdot \sigma &= 0 \quad \text{in } \Omega
\end{aligned} 
\qquad (37) $$

However, for suspensions, definition of the stress tensor involves accounting for the orientation and concentration of fibers, obtained by solving the Fokker-Planck equation. This information is then incorporated into the stress constitutive equation to assess fiber influence on macroscale simulations. Based on the slender body theory [43], the latter can be expressed in the following form:

$$ \sigma = -PI + 2\eta(E + N_p E : a_4) \qquad (38) $$

where $N_p$ is the rheological coupling coefficient. Conservation equations are then solved using stable finite elements using a method based on former works for multiphase flows and for suspension flows [44,45], by modifying the second-hand side to include fiber contribution.

Hence, in order to finalize the system of equations describing the full model, it is also pertinent to reiterate that the Fokker-Planck

equation is employed to calculate the fiber orientation tensor at discrete time steps within each spatial domain. Subsequently, the advection of fibers by the fluid flow completes orientation determination and is implemented adding an advection equation in the physical space, solved using a SUPG-stabilized finite element method, and expressed as follows:

$$ \frac{\partial \varphi}{\partial t} + \mathbf{u} \nabla \varphi = 0 \qquad (39) $$

To test the method, a square domain with a side length of 1 m was taken, discretized into 2092 space nodes, as depicted in Fig. 18. The simulations were conducted under the conditions to obtain a Poiseuille flow, with specified boundary conditions: the pressure at the inlet $P_0$ was set to 10 Pa, the pressure at the outlet $P_1$ was set to 0 Pa, and a zero-slip boundary condition was imposed on the top and bottom wall. The dynamic viscosity of the fluid $\mu$ was defined as 0.15 Pa s. Concerning the advection equation, a completely random fiber orientation distribution was enforced at the inlet as a Dirichlet boundary condition ($\varphi = \frac{1}{2\pi}$). For the Fokker-Planck equation, the angular domain was discretized with a mesh size $h$ of $\frac{2\pi}{4096}$. Temporal discretization employed the Crank-Nicolson scheme, with a time step $\Delta t$ of $5 \times 10^{-3}$ s. The shape factor $\lambda$ was set to 0.98, and the interaction coefficient $C_I$ was set to 0.1. Finally, the simulations were conducted for different rheological coupling coefficients ($N_p = 0, N_p = 2, N_p = 5$ and $N_p = 10$).

The prediction of fiber orientation is initially investigated within the context of an uncoupled flow field scenario ($N_p = 0$). This scenario
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Analytical vs numerical: a comparative study of fourth order fiber Orientation tensor solutions for $\lambda = 1$ and $C_f = 0$.</figcaption></figure>

assumes a steady flow under Poiseuille flow conditions, where the streamlines remain parallel both to each other and to the walls. In this case, depicted in Fig. 19, the time evolution of fiber orientation within the flow is firstly illustrated. Notably, the absence of coupling effects on the conservation equations results in an unchanged velocity profile over time, maintaining the parabolic distribution across the y-direction in the fluid domain.

The ellipses, representing average fiber orientations, are derived from the second-order orientation tensor $\mathbf{a}_2$. The eigenvalues and eigenvectors of $\mathbf{a}_2$ dictate the major axes of the ellipse in the xy-plane, indicating the degree of orientation along these directions. The ellipse color reflects the eigenvalue magnitude in the θ-direction. White signifies isotropic conditions ($a_{xx} = a_{yy} = 0.5$), while black represents full fiber alignment. Close to the centerline axis, fiber orientation appears quasi-random, whereas near the wall, the ellipses flatten and align. Additionally, partial alignment arises due to the nonzero value of $C_f$, preventing perfect fiber alignment along the streamlines, consistent with predictions from the Folgar–Tucker model [14]. Moreover, in the absence of coupling ($N_p = 0$), a fiber suspension flowing between plates can be considered spatially uniform, enabling calculation of dynamic orientation history along a streamline when moving with the fluid.

When accounting for the coupling effect ($N_p \neq 0$), the fiber suspension transforms into a non-Newtonian fluid due to the introduction of particle stress contributions within the momentum Eq. (38). This effect leads to an increase in the total viscosity within the core region, altering the flow behavior. The transition to non-Newtonian behavior is evident in Figs. 20, 21, and 22, where an elevation in $N_p$ corresponds

to a reduction in the maximum velocity. Specifically, as $N_p$ increases from 2 to 10, the velocity profile becomes flatter at its peak, positioned between the core region characterized by unaligned fiber orientation and the wall where fibers align more closely. This flattening of the velocity profile signifies the emergence of non-Newtonian behavior within the flow.

In terms of fiber orientation, as fibers are advected into the domain, those with smaller $N_p$ values reach the outlet more swiftly. This aligns with expectations, as increased fiber concentration results in slower advection. This trend is especially notable when comparing Figs. 19 and 20 with Fig. 22. Additionally, the flow streamlines become progressively disturbed proportional to the rise in $N_p$, gradually returning to parallel alignment as all fibers approach the outlet. Moreover, the coupling between flow and orientation demonstrates minimal impact, primarily due to the dominance of shear flow within the chosen simple geometry (square domain).

**5. Summary and concluding remarks**

A novel 2D numerical model has been constructed, embedded within a finite element framework employing the SUPG method as an alternative to the standard Galerkin method. This model is tailored to explore the intricate motion of fibers, with the potential for future integration into the broader context of flow-fiber coupling, particularly within the scope of molten thermoplastics during the injection phase within mold cavities. The devised Direct Numerical Simulation (DNS) approach is specifically aimed at resolving the Fokker–Planck equation,
<figure><img src="image_19.png" /><figcaption>Fig. 17. Analytical vs numerical: a comparative study of fourth order fiber Orientation tensor solutions for $\lambda = 0.98$ and $C_f = 0$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18. Square domain: FE mesh.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. 19. Evolution of fiber orientation within flow: time-dependent changes in fiber alignment and distribution with $N_p = 0$.</figcaption></figure>

providing a deeper insight into fiber orientation at the mesoscopic scale. Unlike conventional methods reliant on pre-averaged quantities like the second and fourth order orientation tensors ($a_2, a_4$), which often introduce inaccuracies and information loss due to closure approximations, this DNS approach offers a more accurate and detailed representation of fiber orientation dynamics.

The numerical approach accuracy is validated, as a series of rigorous validation tests were conducted. The normalization condition, ensuring the integral of the probability distribution function equals 1, was consistently met. Exploration of parameter dependencies, such as $\lambda$ and $C_I$, yielded results consistent with established expectations. A comprehensive convergence study identified optimal mesh sizes and time steps, striking a balance between accuracy and computational efficiency. Furthermore, comparisons with existing literature, including analytical solutions, underscored the proposed solver's better precision and remarkable agreement with theoretical predictions, affirming its robustness and effectiveness. To extend the proposed approach, the computed orientation tensors are integrated into the fluid problem, creating a comprehensive fiber-flow coupling model. Simulations in a simple geometry under Poiseuille flow conditions are executed to test the coupling effect while accurately accounting for the dynamic fiber motion within the flow.

### CRediT authorship contribution statement

**Nazih Assaad Al Ayoubi:** Methodology, Software, Validation, Writing – original draft, Writing – review & editing. **Hugues Digonnet:** Software, Supervision, Validation, Writing – original draft. **Luisa Silva:**

Investigation, Methodology, Supervision, Validation, Writing – original draft, Writing – review & editing. **Christophe Binetruy:** Funding acquisition, Methodology, Validation, Writing – original draft, Writing – review & editing. **Thierry Renault:** Funding acquisition, Methodology. **Sebastien Comas-Cardona:** Methodology, Supervision, Writing – original draft, Writing – review & editing.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Data availability

Data will be made available on request.

### Acknowledgments

The authors gratefully acknowledge GLiCID Scientific Computing Infrastructure for the access to LIGER supercomputer, GENCI for the access to Joliot Curie supercomputer, and FORVIA for the financial support of this study.
<figure><img src="image_22.png" /><figcaption>Fig. 20. Evolution of fiber orientation within flow: time-dependent changes in fiber alignment and distribution with $N_p = 2$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. 21. Evolution of fiber orientation within flow: time-dependent changes in fiber alignment and distribution with $N_p = 5$.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Fig. 22. Evolution of fiber Orientation within flow: time-dependent changes in fiber alignment and distribution with $N_p$ = 10.</figcaption></figure>
References

[1] A. Andriyana, N. Billon, L. Silva, Mechanical response of a short fiber-reinforced thermoplastic: Experimental investigation and continuum mechanical modeling, Eur. J. Mech. A Solids 29 (6) (2010) 1065-1077.

[2] C. Unterweger, O. Brüggemann, C. Fürst, Synthetic fibers and thermoplastic short-fiber-reinforced polymers: Properties and characterization, Polym. Compos. 35 (2) (2014) 227-236.

[3] C.-T. Huang, X.-W. Chen, W.-W. Fu, Investigation on the fiber orientation distributions and their influence on the mechanical property of the co-injection molding products, Polymers 12 (1) (2019) 24.

[4] U. Gandhi, S. Goris, T.A. Osswald, Y.-Y. Song, Discontinuous Fiber-Reinforced Composites: Fundamentals and Applications, Carl Hanser Verlag GmbH Co KG, 2020.

[5] C.L. Tucker III, Fundamentals of Fiber Orientation: Description, Measurement and Prediction, Carl Hanser Verlag GmbH Co KG, 2022.

[6] G. Ramorino, S. Cecchel, G. Cornacchia, Effect of fiber orientation and residual stresses on the structural performance of injection molded short-fiber-reinforced components, Adv. Syst. Sci. Appl. 20 (2) (2020) 1-19.

[7] H.-C. Tseng, R.-Y. Chang, C.-H. Hsu, Numerical prediction of fiber orientation and mechanical performance for short/long glass and carbon fiber-reinforced composites, Compos. Sci. Technol. 144 (2017) 51-56.

[8] C.L. Tucker III, Planar fiber orientation: Jeffery, non-orthotropic closures, and reconstructing distribution functions, J. Non-Newton. Fluid Mech. 310 (2022) 104939.

[9] A. Moosaie, M. Manhart, A direct numerical simulation method for flow of Brownian fiber suspensions in complex geometries, J. Dispers. Sci. Technol. 34 (3) (2013) 427-440.

[10] S.K. Kugler, A. Kech, C. Cruz, T. Osswald, Fiber orientation predictions—a review of existing models, J. Compos. Sci. 4 (2) (2020) 69.

[11] P.J. Krochak, J. Olson, D. Martinez, The orientation of semidilute rigid fiber suspensions in a linearly contracting channel, Phys. Fluids 20 (7) (2008).

[12] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proc. Royal Soc. Lond. Ser. A Contain. Pap. Math. Phys. Character 102 (715) (1922) 161-179.

[13] S.G. Advani, C.L. Tucker III, The use of tensors to describe and predict fiber orientation in short fiber composites, J. Rheol. 31 (8) (1987) 751-784.

[14] F. Folgar, C.L. Tucker III, Orientation behavior of fibers in concentrated suspensions, J. Reinf. Plast. Compos. 3 (2) (1984) 98-119.

[15] T. Karl, M. Schneider, T. Böhlke, On fully symmetric implicit closure approximations for fiber orientation tensors, J. Non-Newton. Fluid Mech. 318 (2023) 105049.

[16] J.H. Phelps, C.L. Tucker III, An anisotropic rotary diffusion model for fiber orientation in short-and long-fiber thermoplastics, J. Non-Newton. Fluid Mech. 156 (3) (2009) 165-176.

[17] H.-C. Tseng, R.-Y. Chang, C.-H. Hsu, The use of principal spatial tensor to predict anisotropic fiber orientation in concentrated fiber suspensions, J. Rheol. 62 (1) (2018) 313-320.

[18] S.G. Advani, C.L. Tucker III, Closure approximations for three-dimensional structure tensors, J. Rheol. 34 (3) (1990) 367-386.

[19] F. Chinesta, G. Chaidron, A. Poitou, On the solution of Fokker-Planck equations in steady recirculating flows involving short fiber suspensions, J. Non-Newton. Fluid Mech. 113 (2) (2003) 97-125.

[20] K. Chiba, A. Ammar, F. Chinesta, On the fiber orientation in steady recirculating flows involving short fibers suspensions, Rheol. Acta 44 (2005) 406-417.

[21] J.D. Eldredge, Numerical simulation of the fluid dynamics of 2D rigid body motion with the vortex particle method, J. Comput. Phys. 221 (2) (2007) 626-648.

[22] B. Mokdad, E. Pruliere, A. Ammar, F. Chinesta, On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach, Appl. Rheol. 17 (2) (2007) 26494-1-26494-14.

[23] J. Gillissen, B. Boersma, P. Mortensen, H. Andersson, On the performance of the moment approximation for the numerical computation of fiber stress in turbulent channel flow, Phys. Fluids 19 (3) (2007) 035102.

[24] S. Montgomery-Smith, D.A. Jack, D.E. Smith, A systematic approach to obtaining numerical solutions of Jeffery's type equations using spherical harmonics, Composites A 41 (7) (2010) 827-835.

[25] C. Lohmann, Efficient algorithms for constraining orientation tensors in Galerkin methods for the Fokker-Planck equation, Comput. Math. Appl. 71 (5) (2016) 1059-1073.

[26] S. Dahm, C. Helzel, Hyperbolic systems of moment equations describing sedimentation in suspensions of rod-like particles, Multiscale Model. Simul. 20 (3) (2022) 1002-1039.

[27] J. Ferec, D. Mezi, S.G. Advani, G. Ausias, Axisymmetric flow simulations of fiber suspensions as described by 3D probability distribution function, J. Non-Newton. Fluid Mech. 284 (2020) 104367.

[28] J. Férec, M. Heniche, M. Heuzey, G. Ausias, P. Carreau, Numerical solution of the Fokker-Planck equation for fiber suspensions: application to the Folgar-Tucker-Lipscomb model, J. Non-Newton. Fluid Mech. 155 (1-2) (2008) 20-29.

[29] T. Johnson, P. Röyttä, A. Mark, F. Edelvik, Simulation of the spherical orientation probability distribution of paper fibers in an entire suspension using immersed boundary methods, J. Non-Newton. Fluid Mech. 229 (2016) 1-7.

[30] S. Prager, Stress-strain relations in a suspension of dumbbells, Trans. Soc. Rheol. 1 (1) (1957) 53-62.

[31] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, Volume 2: Kinetic Theory, Wiley, 1987.

[32] R. Weinstock, Calculus of Variations: With Applications to Physics and Engineering, Courier Corporation, 1974.

[33] T.J. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, Comput. Methods Appl. Mech. Engrg. 73 (2) (1989) 173-189.

[34] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, et al., PETSc Users Manual, Argonne National Laboratory, 2019.

[35] G.L. Hand, A theory of anisotropic fluids, J. Fluid Mech. 13 (1) (1962) 33-46.

[36] M.C. Altan, L. Tang, Orientation tensors in simple flows of dilute suspensions of non-Brownian rigid ellipsoids, comparison of analytical and approximate solutions, Rheol. Acta 32 (1993) 227-244.

[37] P. Perrochet, D. Bérod, Stability of the standard Crank-Nicolson-Galerkin scheme applied to the diffusion-convection equation: Some new insights, Water Resour. Res. 29 (9) (1993) 3291-3297.

[38] J. Lang, J. Verwer, ROS3P—an accurate third-order Rosenbrock solver designed for parabolic problems, BIT Numer. Math. 41 (4) (2001) 731-738.

[39] J.M. Park, S.J. Park, Modeling and simulation of fiber orientation in injection molding of polymer composites, Math. Probl. Eng. 2011 (2011).

[40] P. Laure, L. Silva, M. Vincent, 20 - Modelling short fibre polymer reinforcements for composites, in: P. Boisse (Ed.), Composite Reinforcements for Optimum Performance, in: Woodhead Publishing Series in Composites Science and Engineering, Woodhead Publishing, 2011, pp. 616-650.

[41] S.K. Kugler, G.M. Lambert, C. Cruz, A. Kech, T.A. Osswald, D.G. Baird, Macroscopic fiber orientation model evaluation for concentrated short fiber reinforced polymers in comparison to experimental data, Polym. Compos. 41 (7) (2020) 2542-2556.

[42] J. Kattinger, S. Joas, F. Willems, M. Kreutzbruck, C. Bonten, Application of the Folgar-Tucker model to predict the orientation of particles of different aspect ratios in polymer suspensions, J. Polym. Eng. 41 (7) (2021) 528-536.

[43] G. Batchelor, Slender-body theory for particles of arbitrary cross-section in Stokes flow, J. Fluid Mech. 44 (3) (1970) 419-440.

[44] L. Silva, T. Coupez, H. Digonnet, Massively parallel mesh adaptation and linear system solution for multiphase flows, Int. J. Comput. Fluid Dyn. 30 (6) (2016) 431-436.

[45] J. Férec, P. Laure, L. Silva, M. Vincent, 20 - Short fiber composite reinforcements, in: P. Boisse (Ed.), Composite Reinforcements for Optimum Performance (Second Edition), second ed., in: Woodhead Publishing Series in Composites Science and Engineering, Woodhead Publishing, 2021, pp. 627-669.