# On the solution of Flokker-Planck equations in steady recirculating flows involving short fibers suspensions

Francisco Chinesta, Guillaume Chaidron, Arnaud Poitou

► To cite this version:

Francisco Chinesta, Guillaume Chaidron, Arnaud Poitou. On the solution of Flokker-Planck equations in steady recirculating flows involving short fibers suspensions. Journal of Non-Newtonian Fluid Mechanics, 2003, 113 (2-3), pp.97-125. 10.1016/S0377-0257(03)00100-9. hal-01007155

HAL Id: hal-01007155

https://hal.science/hal-01007155

Submitted on 11 Dec 2016

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
# On the solution of Fokker-Planck equations in steady recirculating flows involving short fiber suspensions

F. Chinesta<sup>a</sup>, G. Chaidron<sup>b</sup>, A. Poitou<sup>c</sup>

<sup>a</sup> Laboratoire de Mécanique des Systèmes et des Procédés, UMR 8106 CNRS-ENSAM-ESEM, 151 Boulevard de l'Hôpital, F-75013 Paris, France

<sup>b</sup> Institute of Non-Newtonian Fluid Mechanics, The University of Wales, Aberystwyth, Ceredigion, SY23 3BZ Wales, UK

<sup>c</sup> LMM—Ecole Centrale de Nantes, 1 rue de la Noe, BP 92101, 44321 Nantes Cedex 3, France

Numerical modeling of non-Newtonian flows typically involves the coupling between equations of motion characterized by an elliptic character, and the fluid constitutive equation, which is an advection equation linked to the fluid history. Thus, the numerical modeling of short fiber suspensions flows requires a description of the microstructural evolution (fiber orientation) which affects the flow kinematics and that is itself governed by this kinematics. There are different ways to describe the microstructural state. The use of orientation tensors, whose evolution is governed by advection equations, is widely considered. Other possibility is to describe the fiber orientation state from its probability density whose evolution is described by the Fokker-Planck equation (which does not involve any closure relation). The numerical treatment of advection problems is not a simple matter. Moreover, many industrial flows show often steady recirculating areas which introduce some additional difficulties, because in these flows neither boundary conditions nor initial conditions are known. In some of our former papers, we have proposed accurate techniques to solve the linear and non-linear advection equations governing the evolution of the second-order orientation tensor, in steady recirculating flows. These techniques combine a standard treatment of the non-linearity with a more original treatment of the associated linear problems imposing the periodicity condition. In this paper, we generalize those numerical strategies to compute the steady solution of the Fokker-Planck equation (which governs the evolution of the fiber orientation probability distribution) in steady recirculating flows.

**Keywords:** Multidimensional advection equation; Steady recirculating flows; Short fiber suspensions; Fokker-Planck equation; Particle method; Random walks; Fourier spectral method
# 1. Introduction

Numerical modeling of non-Newtonian fluid flows usually involves the coupling between motion equations, which leads to an elliptic problem, and the fluid constitutive equation, which introduces an advection problem related to the fluid history.

## 1.1. Short fiber suspension flow: constitutive equation

In short fiber suspension flow models, the extra-stress tensor $\underline{\tau}$ depends on the fiber orientation [1–4]

$$ \underline{\tau} = 2\mu\{\underline{D} + N_p\underline{a}: \underline{D}\}, \qquad (1) $$

where $\mu$ is the equivalent suspension viscosity, $N_p$ a scalar parameter depending on both the fiber concentration and the fiber aspect ratio, $\underline{D}$ the strain rate tensor, “:” the product defined by $(\underline{a}: \underline{D})_{ij} = a_{ijkl}D_{kl}$, and $\underline{a}$ is the fourth-order orientation tensor defined by

$$ \underline{a} = \oint \underline{p} \otimes \underline{p} \otimes \underline{p} \otimes \underline{p} \Psi(\underline{p}) d\underline{p}, \qquad (2) $$

where $\underline{p}$ is the unit vector aligned in the fiber axis direction, $\otimes$ denotes the tensorial product $(\underline{p} \otimes \underline{p})_{ij} = p_i p_j$ and $\Psi(\underline{p})$ is the orientation distribution function (also known as configuration distribution function or orientation probability density), which satisfies the normality condition

$$ \oint \Psi(\underline{p}) d\underline{p} = 1. \qquad (3) $$

The evolution of the orientation distribution is given by the Fokker-Planck equation. The simplest linear expression of this equation is given by

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial \underline{p}}(\Psi \dot{\underline{p}}) + \frac{\partial}{\partial \underline{p}} \left( D_r \frac{\partial \Psi}{\partial \underline{p}} \right), \qquad (4) $$

where $d/dt$ represents the material derivative, $D_r$ is a diffusion coefficient and $\dot{\underline{p}}$ is the fiber rotation velocity. When the fibers are assumed with an ellipsoidal shape and the suspension is dilute enough, the rotation velocity can be obtained from the Jeffery equation

$$ \dot{\underline{p}} = \underline{\Omega}\underline{p} + k\underline{D}\underline{p} - k(\underline{p}^T\underline{D}\underline{p})\underline{p}, \qquad (5) $$

where $\underline{\Omega}$ is the vorticity tensor associated with the fluid flow undisturbed by the presence of the fiber, and $k$ is a scalar which depends on the fiber aspect ratio $\lambda$ (ratio between the fiber length and the fiber diameter)

$$ k = \frac{\lambda^2 - 1}{\lambda^2 + 1}. \qquad (6) $$

If $\Psi(\underline{p}) = \delta(\underline{p} - \hat{\underline{p}})$, with $\delta()$ the Dirac distribution, all the orientation probability is concentrated in the direction defined by $\hat{\underline{p}}$. In this case we can verify easily that the fourth order orientation tensor can be expressed from the tensorial product of the second-order orientation tensor $\underline{a}$

$$ \underline{a} = \underline{a} \otimes \underline{a}, \qquad (7) $$
where the second-order orientation tensor is defined by

$$ \underline{a} = \oint \underline{p} \otimes \underline{p} \Psi(\underline{p}) d\underline{p}. \qquad (8) $$

From a physical point of view, we can consider that the eigenvalues of the orientation tensor $\underline{a}$ represent the probability of finding the fiber in the direction of the corresponding eigenvectors.

In the case of a concentrated orientation distribution, the extra-stress tensor can be expressed as

$$ \underline{\tau} = 2\mu\{\underline{D} + N_p(\underline{a} : \underline{D})\underline{a}\}. \qquad (9) $$

Now, taking the time derivative of Eq. (8) and introducing the Jeffery's equation (Eq. (5)), the Fokker-Planck equation (Eq. (4)) and the closure relation (Eq. (7)) we obtain a non-linear advection equation governing the evolution of the second-order orientation tensor

$$ \frac{d\underline{a}}{dt} = \underline{\Omega}\underline{a} - \underline{a}\underline{\Omega} + k(\underline{D}\underline{a} + \underline{a}\underline{D} - 2(\underline{a} : \underline{D})\underline{a}) - 4D_r \left(\underline{a} - \frac{\underline{I}}{2}\right). \qquad (10) $$

When the orientation probability is not concentrated in the same direction at each point of the flow domain, the expression (7) is not exact, and in this case it constitutes a quadratic closure relation.

## 1.2. Numerical treatment of models involving advection equations

Usually the solution of a coupled model is found by means of an uncoupled fixed-point strategy. The numerical discretization of the equations of motion does not introduce great difficulties. However, the discretization of the advection equation is more complex. It can be carried out either integrating its Lagrangian description by means of the method of characteristics [5–10], or using the Eulerian discretization of its variational formulation: Streamline Upwind (SU) or Streamline Upwind Petrov-Galerkin (SUPG) finite elements [11–15], discontinuous finite elements [16] or discontinuous finite volumes [15]. However, for the orientation equations as encountered in short fiber suspensions models, the interpolation of the orientation tensors (which is required if standard Eulerian discretization techniques are used) introduces non-physical orientation effects [17].

Such numerical models are made in order to simulate real industrial processes and to predict resulting mechanical properties of conformed parts afterwards. Many of the experimental and industrial flows show recirculating areas or recirculate themselves. For example, in a lid-driven cavity problem, the flow recirculates under the influence of a moving plate. In a contraction or in an expansion flow (as encountered, for example, in extrusion processes), various recirculating areas are observed [18–20]. Many rheometric devices involve this type of flows. Most of these phenomena are associated with a steady state of the flow, which introduces some additional difficulties in the numerical simulation. Actually, the advection equation is supposed to have a steady solution in these steady recirculating flows but neither boundary conditions nor initial conditions are known in such flows.

In a former paper [10], we have proved that, in general two-dimensional steady recirculating flows and for fibers with a quasi-infinite aspect ratio, the local orientation of the fibers with the flow is the only possible solution. In [21] the case of a finite fiber aspect ratio was treated, and we have proved that the existence of a solution in the form $\Psi(x, p) = \delta(p - \hat{p}(x))$ (fiber orientation locally concentrated) depends on the considered flow. However, even if this kind of solution does not exist, a periodic solution along the closed streamlines of the orientation tensor $\underline{a}$ could exist.
Later, in [22] we have proved that some linear steady advection problems defined in steady recirculating flows have only one solution when the kinematics differs from a rigid motion. We also give in [22] a numerical procedure to determine this steady solution, based on the fact that the steady solution must be periodic along the closed trajectories of the flow.

However, even if this result concerns only linear advection equations, the numerical procedure proposed can be applied to solve non-linear advection equations. In this last case the solution is searched by an iterative algorithm, which at each iteration solves a linear advection problem. Thus, accurate and efficient solvers for non-linear advection equations defined in steady recirculating flows are presented in [23].

In these works we have computed different solutions of the fiber orientation in steady recirculating flows using Eq. (10) and the numerical strategies proposed in [22,23], for both quasi-infinite and finite fiber aspect ratios, with and without diffusion effects. However, sometimes the solution found for the second-order orientation tensor is not in agreement with the quadratic closure relation (7), and in these cases we can not conclude about the accuracy of the orientation solution computed.

In this work we propose a different approach which avoids the necessity of proceeding with a closure relation. Thus, the Fokker-Planck equation, Eq. (4), will be solved in steady recirculating flows using some of the ideas reported in [22]. Szeri and Leal [24,25] proposed and applied successfully a similar technique to solve 3D transient and non-recirculating flows involving microstructured fluids. This technique, based on a Lagrangian conservation statement for the distribution function, allows the computation of the Brownian contribution so that, an explicit knowledge of the distribution function is not required.

In this paper, only 2D orientation states will be considered. The discretization of the advection dominated Fokker-Planck equation using a particle technique is treated in the next section, which is followed of a section where some numerical results are presented and discussed. Some numerical strategies based on continuous approximations of the orientation distribution are analyzed in the fourth section, whose associated numerical results are presented in the last section, just before the general conclusions of this paper.

## 2. A particle discretization of the advection dominated Fokker-Planck equation

If we consider the 2D case, the orientation of a fiber is defined by the unit vector $\underline{p}$

$$ \underline{p} = \begin{pmatrix} \cos \varphi \\ \sin \varphi \end{pmatrix}, \qquad (11) $$

where $\varphi$ is the angle between the fiber axis and the x-coordinate axis.

In this case, if we denote by $\Omega_{ij}$ and $D_{ij}$ (with $i, j \in [1, 2]$) the components of the vorticity and the strain rate tensors respectively, the Jeffery equation (Eq. (5)) results

$$ \dot{\varphi} = -\Omega_{12} - kD_{12} - 2kD_{11}\sin\varphi\cos\varphi + 2kD_{12}\cos^2\varphi. \qquad (12) $$

For fibers with infinite aspect ratio, $k = 1$, in simple shear flows ($\underline{v}^T = (\dot{\gamma}y, 0)$, $D_{12} = \Omega_{12} = \dot{\gamma}/2$ and $D_{11} = 0$) the rotation velocity $\dot{\varphi}$ results

$$ \dot{\varphi} = -\dot{\gamma}(1 - \cos^2\varphi), \qquad (13) $$
where $\dot{\gamma}$ is the shear rate. In this case, we can observe that the rotation velocity is only zero for $\varphi = 0$ or $\pi$ which implies the local alignment of the fibers with the flow. In the case of fibers with a finite aspect ratio, $k \neq 1$, the fibers rotate continuously.

The Fokker-Planck equation can be written in the 2D case as

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial\varphi}(\Psi\dot{\varphi}) + \frac{\partial}{\partial\varphi}\left(D_r\frac{\partial\Psi}{\partial\varphi}\right), \qquad (14) $$

where the rotation velocity $\dot{\varphi}$ is given by Eq. (12).

It is important to notice that the orientation distribution $\Psi$ depends on the fiber position $\underline{x}$, on the time $t$ and on the angle $\varphi$ between the fiber axis and the $x$-axis at the considered point. Thus, we will write $\Psi = \Psi(\underline{x}, t, \varphi)$.

In Eq. (14) we observe that the first-term of the right member introduces advection effects, whereas the second one is a diffusion term. Thus, the election of an accurate discretization technique depends on the relation between both terms. The advection dominated case is the most unfavorable scenario from the point of view of its numerical treatment. A numerical strategy appealing in this situation lies on the use of the particle method to integrate accurately the Fokker-Planck equation along the closed streamlines.

With this aim in view, we approximate the steady solution at point $\underline{x}_0$, $\Psi(\underline{x}_0, \varphi) = \Psi^0(\varphi)$, by a discrete distribution of concentrated weights

$$ \Psi^0(\varphi) = \sum_{i=0}^{2N-1} \delta(\varphi - \varphi_i)\Psi^0(\varphi_i), \qquad (15) $$

where $\varphi_i = \pi i / N$ and $\delta(\varphi - \varphi_i)$ denotes the Dirac's distribution, verifying the normality condition

$$ \int_{0}^{2\pi} \Psi^{0}(\varphi) d\varphi = \sum_{i=0}^{2N-1} \Psi^{0}(\varphi_{i}) = 1. \qquad (16) $$

We can consider that $\Psi^0(\varphi_i)$ represents the orientation probability associated with the interval $[\varphi_i - h_\varphi/2, \varphi_i + h_\varphi/2]$, where $h_\varphi = \pi/N$.

Due to the linearity and homogeneity of the Fokker-Planck equation, the solution after a complete turn of period $T$, along the closed streamline related to the point $\underline{x}_0$, results

$$ \Psi_T^0(\varphi) = \sum_{i=0}^{2N-1} \delta(\varphi - \varphi_i^T)\Psi^0(\varphi_i), \qquad (17) $$

where $\varphi_i^T$ denotes the final orientation of the fiber initially orientated in the $\varphi_i$ direction.

To determine $\varphi_i^T$, we take a single fiber located at point $\underline{x}_0$ and aligned in the $\varphi_i$ direction at time $t = 0$, and we compute its angular rotation induced by the flow by using the Jeffery equation, during a complete turn of period $T$ along the closed streamline related to the point $\underline{x}_0$.

Thus, the movement of each individual fiber is subjected to a double advection, the first one along the closed streamline and the second one due to its rotation around its center of mass. If one proceeds by using a fully explicit and first order technique (other possibilities exist, as for example the use of high order explicit or implicit techniques) the fiber position and orientation updating results:
1. Flow advection: The flow advection only changes the position of the center of mass of each fiber, i.e.

$$ (\underline{x}_0^t, \phi_i^t) \rightarrow (\underline{x}_0^{t+\Delta t} = \underline{x}_0^t + \underline{v}(\underline{x}_0^t)\Delta t, \phi_i^t), \qquad (18) $$

where the superscript indicates the time from the beginning of the movement, the subscript “0” denotes the streamline associated with the point $\underline{x}_0$ and the subscript “i” is used to identify the fiber initially aligned in the $\phi_i$ direction.

2. Angular rotation: Due to the flow kinematics the fibers rotate with an angular velocity that can be obtained using the Jeffery Eq. (12)

$$ (\underline{x}_0^{t+\Delta t}, \phi_i^t) \rightarrow (\underline{x}_0^{t+\Delta t}, \phi_i^{t+\Delta t} = \phi_i^t + \dot{\phi}(\phi_i^t)\Delta t). \qquad (19) $$

Obviously, the distribution $\Psi_T^0(\varphi)$ given by Eq. (17) verifies the normality condition, i.e.

$$ \int_{0}^{2\pi} \Psi_{T}^{0}(\varphi) \, d\varphi = \sum_{i=0}^{2N-1} \Psi^{0}(\varphi_{i}) = 1. \qquad (20) $$

Now, we are looking for a discrete expression of the solution $\Psi_T^0(\varphi)$ using the same angles that the initial orientation distribution

$$ \Psi_T^0(\varphi) = \sum_{i=0}^{2N-1} \delta(\varphi - \varphi_i)\alpha^0(\varphi_i). \qquad (21) $$

To transfer the orientation probability concentrated at the angles $\varphi_i^T$, $i \in [0, 2N - 1]$ towards the angles $\varphi_i$, $i \in [0, 2N - 1]$ we proceed as follows: if $\varphi_i^T \in [\varphi_j, \varphi_{j+1}]$ we can transfer the orientation probability related to the angle $\varphi_i^T$ to their neighboring angles $\varphi_j$ and $\varphi_{j+1}$, in the following manner:

$$ \Psi^0(\varphi_i) \frac{\varphi_{j+1} - \varphi_i^T}{h_\varphi} \rightarrow \varphi_j, \qquad (22) $$

$$ \Psi^0(\varphi_i) \frac{\varphi_i^T - \varphi_j}{h_\varphi} \rightarrow \varphi_{j+1}. \qquad (22) $$

Thus, the coefficients $\alpha^0(\varphi_j)$, $j \in [0, 2N - 1]$ result

$$ \alpha^0(\varphi_j) = \sum_{\substack{i=0 \\ \varphi_i^T \in [\varphi_{j-1}, \varphi_{j+1}]} }^{2N-1} \frac{h_\varphi - |\varphi_i^T - \varphi_j|}{h_\varphi} \Psi^0(\varphi_i) = \sum_i F_{ji} \Psi^0(\varphi_i), \qquad (23) $$

where $\varphi_{-1} \equiv \varphi_{2N-1}$ and $\varphi_{2N} \equiv \varphi_0$.

Finally, from Eqs. (21) and (23) we can write

$$ \Psi_T^0(\varphi) = \sum_i \delta(\varphi - \varphi_i)\alpha^0(\varphi_i) = \sum_i \delta(\varphi - \varphi_i) \sum_j F_{ij} \Psi^0(\varphi_j). \qquad (24) $$

It is easy to verify that the new expression of $\Psi_T^0(\varphi)$ verifies the normality condition.

Now, in order to determine the steady solution at point $\underline{x}_0$ of the fiber orientation distribution, we impose the periodicity condition, i.e.

$$ \Psi^0(\varphi_i) = \alpha^0(\varphi_i), \quad \forall i \in [0, 2N - 1], \qquad (25) $$
which can be written in the matrix form

$$ \Psi = F\Psi, \qquad (26) $$

where the *i*-row of $\Psi$ contains the value $\Psi^0(\varphi_i)$. Thus, the solution of $\Psi^0(\varphi_i)$, $\forall i$ is obtained from the system

$$ (F - I)\Psi = 0, \qquad (27) $$

and the normality condition (16).

Before analyzing some numerical results, we introduce the following remarks to the just proposed strategy:

*   **Taking into account slight diffusion effects:** When the diffusion coefficient $D_r$ is not neglected, the Fokker-Planck equation defines an advection-diffusion problem. In this case we can apply again the particle method where the diffusive term is modeled by random motions [26,27] or by using deterministic Lagrangian particle method [24,25,28,29] (see also their associated references). When the diffusion coefficient becomes dominant, discretization strategies using continuous approximations of the orientation distribution function can be accurately applied, as described later in Section 4.

In the following lines we summarize a general particle strategy using random walks in the treatment of slight diffusion effects.

We consider $M$ ($M \ge 1$) individual fibers located at $x_0$ and aligned in each angle $\varphi_i$, $i \in [0, \dots, 2N - 1]$. We are going to describe the numerical procedure for the $M$ fibers located at the initial time in $x_0$ and aligned in the $\varphi_i$ direction. The movement of each one of these fibers is subjected to three actions, that using a fully explicit and first-order strategy can be summarized as follows:

1.  The flow advection only changes the position of the center of mass of each fiber, i.e.

$$ (x_0^t, \varphi_{i,m}^t) \rightarrow (x_0^{t+\Delta t} = x_0^t + \nu(x_0^t)\Delta t, \varphi_{i,m}^t), \quad \forall m \in [1, \dots, M] \qquad (28) $$

where the subscript “m” is used to identify each fiber.

2.  **Angular rotation:** Due to the flow kinematics the fibers rotate with an angular velocity that can be obtained using the Jeffery equation

$$ (x_0^{t+\Delta t}, \varphi_{i,m}^t) \rightarrow (x_0^{t+\Delta t}, \varphi_{i,m}^* = \varphi_{i,m}^t + \dot{\varphi}(\varphi_{i,m}^t)\Delta t), \quad \forall m \in [1, \dots, M] \qquad (29) $$

where $\varphi_{i,m}^*$ is the updated angle.

3.  **Diffusion effects:** The diffusion effects induce an angular rotation which results from a random process

$$ (x_0^{t+\Delta t}, \varphi_{i,m}^*) \rightarrow (x_0^{t+\Delta t}, \varphi_{i,m}^{t+\Delta t} = \varphi_{i,m}^* + \beta(0, \sqrt{2D_r\Delta t})), \quad \forall m \in [1, \dots, M] \qquad (30) $$

where $\beta(0, \sqrt{2D_r\Delta t})$ is a Gaussiany distributed random variable with zero mean and variance $2D_r\Delta t$.

As previously, finally it results

$$ \alpha^0(\varphi_j) = \sum_{i=0}^{2N-1} \sum_{m=1}^{M} \frac{h_{\varphi} - |\varphi_{i,m}^T - \varphi_j|}{Mh_{\varphi}} \Psi^0(\varphi_i). \qquad (31) $$
* *The non-linear Fokker-Planck equation:* In the non-linear case we could proceed firstly by linearizing the Fokker-Planck equation, following some of the ideas described in [23], in order to apply the numerical procedures just proposed to each linear problem.

* *Steady versus transient approach:* Effectively, in process simulation an appealing approach to the recirculation problem would be to use a numerical method that treats transient orientation states, from any initial condition, until reaching the steady state. However, it can be proved that from some initial orientation states the steady solution cannot be reached. For example, if we consider the case of a suspension with $k \neq 1$ and $D_r = 0$, and as initial state the orientation distribution given by

$$ \Psi(\underline{x}, t = 0, \varphi) = \delta(\varphi - \varphi_0), \qquad (32) $$

the fiber orientation distribution at time $t$, resulting from the integration of the Fokker-Planck equation, is

$$ \Psi(\underline{x}, t, \varphi) = \delta(\varphi - \varphi_t), \qquad (33) $$

where

$$ \varphi_t = \varphi_0 + \int_0^t \dot{\varphi}(\varphi(t')) dt'. \qquad (34) $$

We can notice that the fiber orientation does not reach a steady state. However, the Fokker-Planck equation has a steady solution, which results from the integration of

$$ \frac{\partial}{\partial \varphi}(\dot{\varphi}(\varphi)\Psi(\varphi)) = 0 \rightarrow \Psi(\varphi) = \frac{A}{\dot{\varphi}(\varphi)}, \qquad (35) $$

where A must be taken to verify the normality condition

$$ \int_{0}^{2\pi} \Psi(\varphi) d\varphi = 1. \qquad (36) $$

In our opinion, the low enough incidence of this fact in real process simulation, justifies the application of a transient approach.

* *Extension to 3D orientations:* In fact, a 2D fiber orientation state in a steady recirculating flow represents a very idealized situation, with a quite narrow practical interest, but this paper is a first tentative for solving this kind of problems. The proposed technique can be applied to 3D problems, assuming that the fiber orientation is now defined by two angles, $\varphi$ and $\theta$, whose evolution is governed by the Jeffery's equation. In this case the discrete distribution of the orientation probability (used in the particle method) results

$$ \Psi^0(\varphi, \theta) = \sum_{i=0}^{2N-1} \sum_{j=1}^{M-1} \delta(\varphi - \varphi_i) \delta(\theta - \theta_j) \Psi^0(\varphi_i, \theta_j), \qquad (37) $$

where $\varphi_i = \pi i/N$ and $\theta_j = \pi j/M$. However, in the 3D case some specific difficulties are encountered:

* Leal and Hinch state in [30] (Section 3), that when a 3D orientation state is considered, the complete neglect of a Brownian diffusion leaves an indeterminate problem for the steady probability orientation distribution. In the 2D case this problem does not appear because a steady solution given by Eq. (35) exists.

* From the numerical point of view, the use of spherical polar coordinates leads to major difficulties associated with the singularity of this coordinate system at the poles $\theta = \pm \pi/2$ [25].
* *Numerical treatment of general flows:* The proposed numerical strategy can be applied in general flows, containing or not recirculating areas. Thus, if we look for the fiber orientation distribution at point $\underline{x}_0$, we proceed, as previously described, by integrating along its associated streamline in the upstream direction. Then, two possibilities exist:

1. The inflow boundary is reached at point $\underline{X}$. The numerical solution at that point results

$$ \Psi(\underline{X}, \varphi) = \sum_{i=0}^{2N-1} \delta(\varphi - \varphi_i^X) \Psi^0(\varphi_i), \qquad (38) $$

where $\varphi_i^X$ denotes the orientation at point $\underline{X}$ of the fiber initially orientated in the $\varphi_i$ direction at point $\underline{x}_0$. Now, if we identify the solution given by Eq. (38) with the orientation distribution known on the inflow boundary, expressed by

$$ \Psi(\underline{X}, \varphi) = \sum_{i=0}^{2N-1} \delta(\varphi - \varphi_i) \Psi^X(\varphi_i), \qquad (39) $$

we can determine the coefficients $\Psi^0(\varphi_i)$ and then the searched solution at point $\underline{x}_0$.

2. The departing point is reached. In this case, $\underline{x}_0$ is located in a recirculating area, and then, the procedure described in this section can be applied directly.

### 3. The particle method: numerical results and discussion

In this section, some numerical examples are considered using the particle technique just described.

#### 3.1. Simple shear non-recirculating flow

We consider the limit case where the diffusion vanishes ($D_r = 0$) and we are looking for the steady solution for different fiber aspect ratios $0 \le k \le 1$ in the flow defined by

$$ \underline{v} = \begin{pmatrix} y \\ 0 \end{pmatrix}. \qquad (40) $$

The numerical strategy used consists of looking for an initial orientation distribution, such that for any time, the orientation distribution coincides with the initial one. If this orientation distribution does not depend on the time considered, the solution obtained can be considered as the steady solution of the Fokker-Planck equation.

Fig. 1(a) and (b) depict the fiber orientation distribution for fibers with $k = 0.6$ and $0.8$, respectively. In these resolutions we take $N = 1500$ and $D_r = 0$.

We can notice in the previous figures that the orientation probability is maximum around $\varphi = 0$ and $\varphi = \pi$. However, in order to smooth these distributions we take a new angular mesh with $N_s$ intervals $S_i$ of size $h_s$, and we redefine the probability of each one as

$$ \tilde{\Psi}_i = \int_{S_i} \Psi(\varphi) \, d\varphi, \quad \forall i \in [1, \dots, N_s]. \qquad (41) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Steady fiber orientation in a simple and non-recirculating shear flow without diffusion effects ($D_t = 0$) computed using the particle method: from top to bottom (a) $\Psi(\varphi)$, $k = 0.6$; (b) $\Psi(\varphi)$, $k = 0.8$; (c) smoothed $\Psi(\varphi)$, $k = 0.6$; (d) smoothed $\Psi(\varphi)$, $k = 0.8$; (e) reference solution $\Psi(\varphi)$, $k = 0.6$ and (f) reference solution $\Psi(\varphi)$, $k = 0.8$.</figcaption></figure>

Thus, the smoothing distributions corresponding to $k = 0.6$ and $0.8$ are represented in Fig. 1(c) and (d), respectively.

In order to evaluate the accuracy of these solutions we can compute a good approximation of the exact solution, which will be used as reference solution, using the ergodic principle. Following this principle, the steady orientation distribution, $\Psi(\varphi)$, coincides in a simple shear flow with the spatial average of a single fiber orientation evolution, whose discrete form results

$$ \Psi_{\text{ref}}(\varphi_i) = \frac{t_i}{T}, \quad \forall i \in [0, \dots, 2N-1], \qquad (42) $$

where $T$ is the rotation period and $t_i$ is the time spent by the fiber during its rotation into the interval $[\varphi_i - h_\varphi/2, \varphi_i + h_\varphi/2]$. We can notice that the previous definition verifies the discrete normality condition.

The reference solutions for $k = 0.6$ and $0.8$ are shown in Fig. 1(e) and (f), respectively. In the case of fibers with quasi-infinite aspect ratio ($k = 1$), all the orientation probability is concentrated, as expected, at the angles $\varphi = 0$ and $\pi$.

To compare the reference solutions and the ones computed with the particle strategy just described, we proceed to define and evaluate the following discrete error which depends on the angular mesh, the time
step and the fiber aspect ratio

$$E(N, \Delta t, k) = \sum_{i=0}^{2N-1} |\Psi(\varphi_i) - \Psi_{\text{ref}}(\varphi_i)|. \qquad (43)$$

We have verified that the error decreases as the time step decreases and the angular degrees of freedom increase. In the analyzed cases, we obtain

$$\begin{cases} E(1500, 0.001, 0.6) = 0.01 \\ E(1500, 0.001, 0.8) = 0.03, \end{cases} \qquad (44)$$

which proves the high accuracy and the negligible numerical diffusion characteristic of the particle method applied to advection problems.

**Remark.** It can be noticed that a steady fiber orientation distribution $\Psi(\varphi)$ exists for $k \neq 1$ and $D_r = 0$, verifying

$$\frac{\partial}{\partial \varphi}(\dot{\varphi}(\varphi)\Psi(\varphi)) = 0, \qquad (45)$$

whose solution has been introduced in Eqs. (35) and (36), in spite of the fact that a single fiber immersed in such flow never reaches a steady orientation, according to Eq. (12).

## 3.2. Shear recirculating flow

We consider the flow defined by the following velocity field

$$\underline{v} = \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} -y\sqrt{x^2 + y^2} \\ x\sqrt{x^2 + y^2} \end{pmatrix}, \qquad (46)$$

and the limit case where all the diffusion effects are neglected ($D_r = 0$).

We are looking for the steady solution of the orientation distribution at point $x = 0, y = 1$; imposing the periodicity condition and the particle technique.

When the fibers have an infinite aspect ratio $k = 1$, we have proved in a former work, [10], that the orientation distribution is concentrated at the angles $\varphi = 0$ and $\pi$ which corresponds to the local alignment of the fibers with the flow, i.e. the orientation distribution can be written as

$$\Psi(\varphi) = \delta(\varphi - \varphi(t)), \qquad (47)$$

where $\varphi(t)$ represents the angle associated with the unit vector $\underline{t}$, tangent to the streamline at the considered point

$$\underline{t} = \frac{\underline{v}}{\|\underline{v}\|}. \qquad (48)$$

For other aspect ratios $k < 1$, we have also proved [21] that in general flows the orientation distribution solution at any point cannot be expressed by means of a Dirac's function, i.e. the orientation probability is not concentrated in a certain angle.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Smoothed fiber distribution in a recirculating shear flow at point (0, 1) without diffusion effects $D_r = 0$ computed using the particle method: from top to bottom (a) fibers with $k = 0.6$ and (b) fibers with $k = 0.8$.</figcaption></figure>

Thus, for $k = 0.6$ and $0.8$, Fig. 2(a) and (b) show the smoothed steady orientation distribution obtained using the particle method with $N = 1500$, $N_s = 100$, $\Delta t = 0.001$ and $D_r = 0$ at point $x = 0$, $y = 1$.

With the steady orientation distribution computed at any point of the steady recirculating flow, we can evaluate easily the second and fourth-order orientation tensors, from the discrete form of the expressions (8) and (2)

$$ \left\{ \begin{aligned} a_{11} &= \sum_{i=0}^{2N-1} \cos^2(\varphi_i) \Psi(\varphi_i) \\ a_{12} &= \sum_{i=0}^{2N-1} \cos(\varphi_i) \sin(\varphi_i) \Psi(\varphi_i) \\ a_{22} &= \sum_{i=0}^{2N-1} \sin^2(\varphi_i) \Psi(\varphi_i), \end{aligned} \right. \qquad (49) $$

with $a_{12} = a_{21}$ and

$$ \left\{ \begin{aligned} a_{1111} &= \sum_{i=0}^{2N-1} \cos^4(\varphi_i) \Psi(\varphi_i) \\ a_{1112} &= \sum_{i=0}^{2N-1} \cos^3(\varphi_i) \sin(\varphi_i) \Psi(\varphi_i) \\ a_{1122} &= \sum_{i=0}^{2N-1} \cos^2(\varphi_i) \sin^2(\varphi_i) \Psi(\varphi_i) \\ a_{1222} &= \sum_{i=0}^{2N-1} \cos(\varphi_i) \sin^3(\varphi_i) \Psi(\varphi_i) \\ a_{2222} &= \sum_{i=0}^{2N-1} \sin^4(\varphi_i) \Psi(\varphi_i), \end{aligned} \right. \qquad (50) $$
with

$$ \left\{ \begin{array}{l} a_{1112} = a_{1121} = a_{1211} = a_{2111} \\ a_{1122} = a_{2211} = a_{1221} = a_{2112} = a_{1212} = a_{2121} \\ a_{1222} = a_{2122} = a_{2212} = a_{2221}. \end{array} \right. \qquad (51) $$

Thus, for $k = 0.6$ we obtain the following expression for the second-order orientation tensor

$$ \underline{q}(0, 1) = \begin{pmatrix} 0.66 & 0.00 \\ 0.00 & 0.33 \end{pmatrix}. \qquad (52) $$

To evaluate the accuracy of this solution, we can compute the exact solution, which can be obtained due to the very simple kinematics considered. Due to the flow symmetry, the orientation distribution in cylindrical coordinates must be the same along the closed streamlines. The Lagrangian derivative results in this case

$$ \frac{d\Psi}{dt} = \frac{\partial\Psi}{\partial\varphi}\omega, \qquad (53) $$

where $\omega$ represents the angular velocity related to the considered streamline, that in our case results $\omega = -1$. Thus, the Fokker-Planck equation without diffusion effects results in this case

$$ \frac{\partial}{\partial \varphi}((\varphi - 1.0)\Psi) = 0, \qquad (54) $$

where for fibers with $k = 0.6$, the angular velocity of the fibers at point (0, 1) results from Eq. (12)

$$ \dot{\varphi} = 1.8 - 0.6 \cos^2 \varphi. \qquad (55) $$

The integration of Eq. (54) results in

$$ \Psi(\varphi) = \frac{A}{\dot{\varphi} - 1.0} = \frac{A}{0.8 - 0.6 \cos^2 \varphi}, \qquad (56) $$

where A does not depend on the angular coordinate $\varphi$. Now, the value of A can be calculated to verify the normality condition

$$ 1 = \int_{0}^{2\pi} \Psi(\varphi) \, d\varphi = A \int_{0}^{2\pi} \frac{1}{0.8 - 0.6 \cos^2 \varphi} \, d\varphi = 15.69A. \qquad (57) $$

With A known, the orientation distribution becomes perfectly defined, and we can proceed to obtain the value of the different orientation tensors. In particular, the second-order orientation tensor at point (0, 1) results

$$ \left\{ \begin{array}{l} a_{11} = \int_{0}^{2\pi} \cos^2(\varphi) \Psi d\varphi = \int_{0}^{2\pi} \cos^2(\varphi) \frac{A}{\dot{\varphi}} d\varphi = 0.66 \\ a_{12} = \int_{0}^{2\pi} \cos(\varphi) \sin(\varphi) \Psi d\varphi = \int_{0}^{2\pi} \cos(\varphi) \sin(\varphi) \frac{A}{\dot{\varphi}} d\varphi \approx 0 \\ a_{22} = \int_{0}^{2\pi} \sin^2(\varphi) \Psi d\varphi = \int_{0}^{2\pi} \sin^2(\varphi) \frac{A}{\dot{\varphi}} d\varphi = 0.33, \end{array} \right. \qquad (58) $$

which proves the accuracy of the particle strategy described in this paper.
As commented in the first section, sometimes the fourth-order orientation tensor is assumed to derive from the second-order one (quadratic closure relation), i.e.

$$ \underline{a} = \underline{a} \otimes \underline{a} \quad (a_{ijkl} = a_{ij}a_{kl}, \quad i, j, k, l \in [1, 2]). \qquad (59) $$

In order to evaluate the accuracy of the quadratic closure relation we define the ratio $E_{qc}(\underline{x}, k)$ at point $\underline{x}$ and for fibers with aspect ratio $k$

$$ E_{qc} = \frac{\| \underline{a} - \underline{a} \otimes \underline{a} \|_{L^2}}{\| \underline{a} \|_{L^2}}, \qquad (60) $$

where $\| \underline{b} \|_{L^2}^2 = \sum_{i,j,k,l} b_{ijkl}^2$ for any fourth-order tensor $\underline{b}$.

We obtain $E_{qc}((0, 1), 0.6) = 0.48$ which proves that the quadratic closure relation is not very accurate for fibers with $k = 0.6$. As commented in the introduction section, when the fibers have an infinite aspect ratio, $k = 1$, the quadratic closure relation is exact. Thus, for $k = 1$ we obtain at point $(0, 1)$

$$ \underline{a}(0, 1) = \begin{pmatrix} 1.0 & 0.0 \\ 0.0 & 0.0 \end{pmatrix}, \qquad (61) $$

and $E_{qc}((0, 1), 1.0) = 0.0$.

For fibers with $k = 0.8$, the second orientation tensor associated with the orientation distribution depicted in Fig. 2(b) results

$$ \underline{a}(0, 1) = \begin{pmatrix} 0.75 & 0.00 \\ 0.00 & 0.25 \end{pmatrix}, \qquad (62) $$

and $E_{qc}((0, 1), 1.0) = 0.37$. Effectively, more the orientation distribution is concentrated, better is the quadratic closure relation.

If one looks for the steady solution of the second-order orientation tensor imposing the periodicity in its evolution Eq. (10), which involves the quadratic closure relation, following the procedure described in [22], the solution obtained is far from the reference solution found using the particle procedure. Thus, the steady solution of Eq. (10) with $D_r = 0$ for the shear recirculating flow defined in Eq. (46), obtained by imposing its periodicity along the closed streamline associated with the point $x = 0, y = 1$; for fibers with $k = 0.6$, results [22]

$$ \underline{a}(0, 1) = \begin{pmatrix} 0.80 & 0.00 \\ 0.00 & 0.20 \end{pmatrix}, \qquad (63) $$

with a significant deviation from the expression (52) obtained using the particle method to impose the solution periodicity on the Fokker-Planck equation.

From these results, we can conclude that the steady solution of the fiber orientation, for fibers with finite aspect ratio, in steady recirculating flows must be searched preferably without introducing the quadratic closure relation.
In this work only the quadratic closure relation is considered for illustrating the utility of accurate solutions of the fiber orientation distribution in steady recirculating flows. In [31] we have computed the error introduced by using different closure relations (quadratic, linear, hybrid [32] and natural [33]) in some simple flows. In those preliminary results the most accurate solutions are obtained when the natural closure relation is considered.

### 3.3. Introducing slight diffusion

In this example we consider the flow kinematics defined by Eq. (46) and a suspension characterized by $k = 0.6$ and $D_r = 0.18$. The particle method has been applied for solving the Fokker-Planck equation, according to the algorithm described in Section 2, and the computed solutions at point (0, 1) for different values of $M$ are shown in Fig. 3.

In order to evaluate the accuracy of these numerical solutions we compute a reference solution. For this, we consider that the fiber orientation distribution does not change along the streamline if a cylindrical coordinates system is considered. In this case, at point (0, 1) the Fokker-Planck

<figure><img src="image_4.png" /><figcaption>Fig. 3. Numerical solution of the fiber orientation distribution ($k = 0.6$ and $D_r = 1.8$) in a shear recirculating flow at point (0, 1) computed using the particle method coupled with a random walk strategy using $M = 1, 10, 100$ and 400.</figcaption></figure>
equation becomes

$$D_r \frac{d^2 \Psi}{d\varphi^2} - (\dot{\varphi}|_{(0,1)} - 1.0) \frac{d\Psi}{d\varphi} - \left( \left. \frac{d\dot{\varphi}}{d\varphi} \right|_{(0,1)} \right) \Psi = 0, \qquad (64)$$

where the solution $\Psi(\varphi)$ at point (0, 1) must verify the normality condition, i.e.

$$\int_{0}^{2\pi} \Psi((0, 1), \varphi) d\varphi = 1. \qquad (65)$$

Eqs. (64) and (65) can be solved accurately by finite differences taking into account the solution periodicity, leading to the reference solution which will be used to evaluate the numerical solutions depicted in Fig. 3. We have verified that the error norm decreases as the value of $M$ increases, which proves the convergence of the particle method coupled with a random walk strategy to take into account slight diffusion effects. However, the number of particles $M$ required for a given error threshold increases significantly with the diffusion coefficient. Thus, in the numerical example described in this section, more than 400 fibers for the 1600 angles considered (640,000 fibers) are required to describe accurately the orientation distribution at point (0, 1). For this reason, for practical applications, the use of the particle method is restricted to slight diffusion effects.

## 4. Numerical discretization strategies based on continuous approximations of the orientation distribution

We consider again the discretization of the 2D Fokker-Planck equation

$$\frac{d\Psi}{dt} = -\frac{\partial}{\partial\varphi}(\Psi\dot{\varphi}) + \frac{\partial}{\partial\varphi}\left(D_r\frac{\partial\Psi}{\partial\varphi}\right), \qquad (66)$$

using different numerical schemes based on continuous approximations of the orientation distribution function. These strategies are appealing in the case of considering diffusion coefficients lying in a large interval.

### 4.1. Centered schemes: a Fourier spectral technique

When the order of magnitude of the convection term is lower than the diffusion one, we can use centered schemes in the discretization of the convection term. As the orientation distribution must be periodic, we can use for the angular discretization a Fourier spectral technique [34]. Thus, the fiber orientation distribution for any angle $\varphi$ can be interpolated from their values in the $2N$ angles $\varphi_i$

$$\varphi_i = \frac{\pi i}{N}, \quad i \in [0, \dots, 2N-1], \qquad (67)$$

from the expression

$$\Psi(\underline{x}, t, \varphi) = \sum_{i=0}^{2N-1} \Psi_i(\underline{x}, t) N_i(\varphi), \qquad (68)$$
where $\Psi_i(\underline{x}, t) = \Psi(\underline{x}, t, \varphi_i)$ and the Fourier's interpolation functions $N_i(\varphi)$ are given by

$$N_i(\varphi) = \frac{1}{2N} \sin(N(\varphi - \varphi_i)) \cot\left(\frac{\varphi - \varphi_i}{2}\right), \qquad (69)$$

whose different derivatives will be denoted by $(Dk)_{ji}$

$$(Dk)_{ji} = \left. \frac{d^k N_i(\varphi)}{d\varphi^k} \right|_{\varphi=\varphi_j}, \qquad (70)$$

where

$$(D1)_{ji} = \begin{cases} \frac{1}{2}(-1)^{j+i} \cot\left(\frac{\varphi_j - \varphi_i}{2}\right), & j \neq i \\ 0, & j = i, \end{cases} \qquad (71)$$

and

$$(D2)_{ji} = \begin{cases} \frac{1}{2}(-1)^{j+i+1} \frac{1}{\sin^2(\varphi_j - \varphi_i/2)}, & j \neq i \\ -\frac{2N^2 + 1}{6}, & j = i. \end{cases} \qquad (72)$$

Using a collocation technique at each angle $\varphi_i$, Eq. (66) results in the following linear system of ordinary differential equations (ODE):

$$\begin{aligned} \frac{d\Psi_i(\underline{x}, t)}{dt} = & -\dot{\varphi}|_{\varphi=\varphi_i} \sum_{j=0}^{2N-1} (D1)_{ij} \Psi_j(\underline{x}, t) - \left.\frac{d\dot{\varphi}}{d\varphi}\right|_{\varphi=\varphi_i} \Psi_i(\underline{x}, t) \\ & + D_r \sum_{j=0}^{2N-1} (D2)_{ij} \Psi_j(\underline{x}, t), \quad i \in [0, \dots, 2N-1], \end{aligned} \qquad (73)$$

which can be written in the matrix form

$$\frac{d}{dt} \Psi = A^{\text{cen}} \Psi, \qquad (74)$$

where the superscript “cen” indicates that the discretization has been carried out by using a centered scheme. In the previous expression $\Psi$ is given by

$$\Psi = \begin{pmatrix} \Psi_0(\underline{x}, t) \\ \Psi_1(\underline{x}, t) \\ \vdots \\ \Psi_{2N-1}(\underline{x}, t) \end{pmatrix}, \qquad (75)$$

and the i, j-component of $A^{\text{cen}}$, $A^{\text{cen}}_{ij}$, results

$$A^{\text{cen}}_{ij} = -\dot{\varphi}|_{\varphi=\varphi_i} (D1)_{ij} - \left.\frac{d\dot{\varphi}}{d\varphi}\right|_{\varphi=\varphi_i} \delta_{ij} + D_r (D2)_{ij}, \qquad (76)$$
where $\delta_{ij}$ is the Kroenecker's delta. The components of $A^{\text{cen}}$ depend on the point $\bar{x}$ due to the dependence of the fiber rotation velocity $\dot{\varphi}$ and its gradient $d\dot{\varphi}/d\varphi$ on both the vorticity and the strain rate tensors, which obviously differ in general flows from one point to another.

## 4.2. Stabilisation by upwinding

When lower diffusion effects are considered, it is well known that centered schemes are inaccurate. In this case, upwinding is widely applied in order to stabilize the numerical discretization. An alternative formula for the convective derivative is

$$ \left. \frac{\partial(\Psi \dot{\varphi})}{\partial \varphi} \right|_{\varphi=\varphi_i} = \frac{\Psi_i \dot{\varphi}|_{\varphi=\varphi_i} - \Psi_{i-1} \dot{\varphi}|_{\varphi=\varphi_{i-1}}}{\varphi_i - \varphi_{i-1}}, \quad \text{if } \dot{\varphi}|_{\varphi=\varphi_i} > 0 \qquad (77) $$

or

$$ \left. \frac{\partial(\Psi \dot{\varphi})}{\partial \varphi} \right|_{\varphi=\varphi_i} = \frac{\Psi_{i+1} \dot{\varphi}|_{\varphi=\varphi_{i+1}} - \Psi_i \dot{\varphi}|_{\varphi=\varphi_i}}{\varphi_{i+1} - \varphi_i}, \quad \text{if } \dot{\varphi}|_{\varphi=\varphi_i} < 0 \qquad (78) $$

where in these equations $\Psi_{2N} \equiv \Psi_0$ and $\Psi_{-1} \equiv \Psi_{2N-1}$.

The second derivative can be discretized using the finite difference

$$ \left. \frac{\partial^2 \Psi}{\partial \varphi^2} \right|_{\varphi=\varphi_i} = \frac{\Psi_{i+1} - 2\Psi_i + \Psi_{i-1}}{h_{\varphi}^2}, \qquad (79) $$

where $h_{\varphi} = 2\pi/2N = \pi/N$.

As previously, the final linear system of ODE can be written in the form

$$ \frac{d}{dt} \Psi = A^{\text{up}} \Psi \qquad (80) $$

where the superscript “up” indicates that the discretization has been carried out by using an upwinding technique.

**Remark.** Other higher order stabilizations exist, however, in this work, we will consider the simplest first-order schema just described, in spite of its higher numerical diffusion.

## 4.3. Imposing the periodicity along the closed streamlines

Now, the resulting linear system of ODE can be integrated using the method of characteristics along the closed trajectories (which correspond with the streamlines in steady flows). The only difficulty is that, as the trajectories are closed, we do not have an initial condition to start this integration. This difficulty was treated in our previous work [22].

Let

$$ \frac{d}{dt} \bar{\Psi} = A \Psi, \qquad (81) $$

be the system of equations to solve, where $A$ is $A^{\text{cen}}$ or $A^{\text{up}}$, depending on the considered discretization technique, and $\Psi^n$ is the solution of Eq. (81) when the initial condition
$$ \Psi^n(\underline{x}_0) = \begin{pmatrix} \Psi_0^n(\underline{x}_0) \\ \Psi_1^n(\underline{x}_0) \\ \vdots \\ \Psi_{n-1}^n(\underline{x}_0) \\ \Psi_n^n(\underline{x}_0) \\ \Psi_{n+1}^n(\underline{x}_0) \\ \vdots \\ \Psi_{2N-1}^n(\underline{x}_0) \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}, \qquad (82) $$

is considered at point $\underline{x}_0$.

Then, the solution $\Psi^n(\underline{x}(\underline{x}_0, t))$ can be computed using the method of characteristics at any point of the closed trajectory $\underline{x} = \underline{x}(\underline{x}_0, t)$. The streamline passing through $\underline{x}_0$ can be reconstructed by integrating

$$ \underline{v}(\underline{x}) = \frac{d\underline{x}}{dt}, \quad \underline{x}(t=0) = \underline{x}_0, \qquad (83) $$

where $\underline{v}$ is the velocity field.

In order to integrate Eq. (81) from the initial condition imposed at point $\underline{x}_0$ the simplest strategy is the use of a first-order fully explicit technique

$$ \left\{ \begin{aligned} \frac{d}{dt} \Psi(\underline{x}(\underline{x}_0, m \Delta t)) &= \underline{A} \Psi(\underline{x}(\underline{x}_0, m \Delta t)) \\ \Psi(\underline{x}(\underline{x}_0, (m+1) \Delta t)) &= \Psi(\underline{x}(\underline{x}_0, m \Delta t)) + \frac{d}{dt} \Psi(\underline{x}(\underline{x}_0, m \Delta t)) \Delta t \\ \underline{x}(\underline{x}_0, (m+1) \Delta t) &= \underline{x}(\underline{x}_0, m \Delta t) + \underline{v}(\underline{x}(\underline{x}_0, m \Delta t)) \Delta t. \end{aligned} \right. \qquad (84) $$

Higher order fully explicit techniques (Runge-Kutta, etc.) are also available, as well as different types of semi-implicit or fully-implicit techniques. In all cases the time step must be adjusted to guarantee stability, and convergence must be checked.

Now, with the $2N$ solutions $\Psi^n$, $n \in [0, \dots, 2N-1]$ computed, the general solution $\Psi$ of Eq. (81)), taking into account its linearity and homogeneity, results

$$ \Psi = \sum_{n=0}^{2N-1} \alpha_n \Psi^n, \quad \alpha_n \in \mathbb{R}. \qquad (85) $$

In order to determine the coefficients $\alpha_n$, we enforce the periodicity of $\Psi$ on the closed trajectory, i.e. $\Psi(\underline{x}_0) = \Psi(\underline{x}(\underline{x}_0, T))$, where $T$ is the period of the closed trajectory considered. Thus, this condition implies

$$ \sum_{n=0}^{2N-1} \alpha_n \Psi^n(t=T) = \sum_{n=0}^{2N-1} \alpha_n \Psi^n(t=0) = \sum_{n=0}^{2N-1} \underline{\alpha}, \qquad (86) $$
where $\alpha$ is defined by

$$ \alpha = \begin{pmatrix} \alpha_0 \\ \alpha_1 \\ \vdots \\ \alpha_{2N-1} \end{pmatrix}. \qquad (87) $$

The periodicity condition expressed by Eq. (86) can be written in the following matrix form

$$ (\underline{B} - \underline{I})\underline{\alpha} = \underline{0}, \qquad (88) $$

where the n-column of $\underline{B}$ contains the vector $\Psi^n(t = T)$ and $\underline{I}$ is the unit matrix. If we define the matrix $\underline{C}$ as $\underline{C} = \underline{B} - \underline{I}$, then the previous linear system results

$$ \underline{C}\underline{\alpha} = \underline{0}. \qquad (89) $$

Due to the homogeneous character of this system, two possibilities exist:

1. First possibility

$$ \det \underline{C} \neq 0, \qquad (90) $$

which implies that the only solution is given by $\underline{\alpha} = \underline{0}$. The orientation distribution associated with $\underline{\alpha} = \underline{0}$ is $\underline{\Psi} = \underline{0}$ that is in contradiction with the normality property of the fiber orientation distribution given by Eq. (3).

2. Second possibility

$$ \det \underline{C} = 0. \qquad (91) $$

In this case one of the equations in the system (89) can be written as a linear combination of the others. So, we must eliminate one of the equations, imposing at its place the discrete form of the normality condition (3). To this purpose we can write, taking into account the regular distribution of the collocation points

$$ \int_{0}^{2\pi} \Psi(x_0, \varphi) \, d\varphi = \sum_{i=0}^{2N-1} \alpha_i h_{\varphi} = 1, \qquad (92) $$

where $h_{\varphi} = \pi/N$.

Thus, the modified linear system results

$$ \underline{K}\underline{\alpha} = \underline{F}, \qquad (93) $$

where $\underline{K}$ results from $\underline{C}$ replacing for example in its last row each component by $h_{\varphi}$, and $\underline{F}$ contains zeros except its last component which is the unit value.

From the solution of the linear system (93) we obtain the vector $\underline{\alpha}$, from which the orientation distribution at the point $x_0$ results finally perfectly known.
### 4.4. An efficient and accurate implicit strategy

Standard implicit formulations require the computation of a $2N \times 2N$ inverse matrix at each time step. In this section, we propose an implicit numerical scheme which does not require any matrix inversion, with a significant reduction in the CPU time.

#### 4.4.1. An implicit procedure

The 2D Fokker-Planck Eq. (66) can be written as

$$ \frac{d\Psi}{dt} = \mathbb{L}_{\varphi}(\Psi), \qquad (94) $$

where $\mathbb{L}_{\varphi}$ represents the linear differential operator

$$ \mathbb{L}_{\varphi}(\Psi(x, \varphi)) = -\frac{\partial}{\partial \varphi}(\Psi \dot{\varphi}) + \frac{\partial}{\partial \varphi} \left( D_r \frac{\partial \Psi}{\partial \varphi} \right). \qquad (95) $$

Now, the implicit time discretization of Eq. (94) results

$$ \frac{\Psi^{n+1}(\varphi) - \Psi^n(\varphi)}{\Delta t} = \mathbb{L}_{\varphi}(\Psi^{n+1}(\varphi)), \qquad (96) $$

where $\Psi^n(\varphi)$ is known from the initial condition or from the solution computed at the previous iteration at point $x = x_0, n\Delta t$, and $\Psi^{n+1}(\varphi)$ represents the unknown solution at point $x = x_0, (n + 1)\Delta t$.

Eq. (96) can be rewritten as

$$ -\Delta t \mathbb{L}_{\varphi}(\Psi^{n+1}(\varphi)) + \Psi^{n+1}(\varphi) = \Psi^n(\varphi). \qquad (97) $$

Now, we define the new linear differential operator $\mathbb{G}_{\varphi}$

$$ \mathbb{G}_{\varphi}(\Psi(\varphi)) = -\Delta t \mathbb{L}_{\varphi}(\Psi(\varphi)) + \Psi(\varphi), \qquad (98) $$

such that, Eq. (97) results

$$ \mathbb{G}_{\varphi}(\Psi^{n+1}(\varphi)) = \Psi^n(\varphi), \qquad (99) $$

with $\Psi^{n+1}(\varphi) = 2\pi$-periodic.

The periodicity condition can be treated using the strategy previously described, and that we summarize in the following lines.

We define the functions $\Psi^{1,h}$ and $\Psi^{2,h}$ solution of the problems:

$$ \mathbb{G}_{\varphi}(\Psi^{1,h}(\varphi)) = 0, \quad \Psi^{1,h}(\varphi = 0) = 1 \quad \text{and} \quad \left. \frac{d\Psi^{1,h}}{d\varphi} \right|_{\varphi=0} = 0, \qquad (100) $$

and

$$ \mathbb{G}_{\varphi}(\Psi^{2,h}(\varphi)) = 0, \quad \Psi^{2,h}(\varphi = 0) = 0 \quad \text{and} \quad \left. \frac{d\Psi^{2,h}}{d\varphi} \right|_{\varphi=0} = 1, \qquad (101) $$

respectively. In this case, the general solution of the Eq. (99) without second member, $\Psi^h(\varphi)$, is given by

$$ \Psi^h(\varphi) = \alpha_1 \Psi^{1,h}(\varphi) + \alpha_2 \Psi^{2,h}(\varphi), \quad \alpha_1, \alpha_2 \in \mathbb{R}. \qquad (102) $$
A particular solution of Eq. (99), $\Psi^c(\varphi)$, can be obtained from

$$ \mathbb{G}_{\varphi}(\Psi^c(\varphi)) = \Psi^n(\varphi), \quad \Psi^c(\varphi = 0) = 0 \quad \text{and} \quad \left. \frac{d\Psi^c}{d\varphi} \right|_{\varphi=0} = 0. \qquad (103) $$

All these solutions $\Psi^{1,h}(\varphi)$, $\Psi^{2,h}(\varphi)$ and $\Psi^c(\varphi)$ can be obtained easily using finite differences. When the integration is carried out in the downstream direction, using upstream finite difference formulas, stabilization by upwinding is introduced.

In this form, the general solution of Eq. (99), $\Psi^{n+1}(\varphi)$, results

$$ \Psi^{n+1}(\varphi) = \Psi^c(\varphi) + \Psi^h(\varphi) = \Psi^c(\varphi) + \alpha_1 \Psi^{1,h}(\varphi) + \alpha_2 \Psi^{2,h}(\varphi). \qquad (104) $$

Taking into account the definition of the previous problems, $\Psi^{n+1}(\varphi = 0) = \alpha_1$ and $((d\Psi^{n+1})/(d\varphi))|_{\varphi=0} = \alpha_2$. The periodicity condition requires the verification of the relations

$$ \left\{ \begin{array}{l} \Psi^{n+1}(\varphi = 0) = \Psi^{n+1}(\varphi = 2\pi) \\ \left. \frac{d\Psi^{n+1}}{d\varphi} \right|_{\varphi=0} = \left. \frac{d\Psi^{n+1}}{d\varphi} \right|_{\varphi=2\pi} \end{array} \right. \qquad (105) $$

These conditions can be expressed by the following linear system

$$ \begin{pmatrix} 1 - \Psi^{1,h}(2\pi) & -\Psi^{2,h}(2\pi) \\ -\Psi_{,\varphi}^{1,h}(2\pi) & 1 - \Psi_{,\varphi}^{2,h}(2\pi) \end{pmatrix} \begin{pmatrix} \alpha_1 \\ \alpha_2 \end{pmatrix} = \begin{pmatrix} \Psi^c(2\pi) \\ \Psi_{,\varphi}^c(2\pi) \end{pmatrix}, \qquad (106) $$

where $\Psi_{,\varphi} \equiv d\Psi/d\varphi$. The scalars $\alpha_1$ and $\alpha_2$, which define completely the solution $\Psi^{n+1}(\varphi)$, are obtained solving the previous linear system.

### 4.4.2. Integration along the closed streamlines

As previously described, we can approximate the steady solution at $\underline{x}_0$, $\Psi(\underline{x}_0, \varphi) = \Psi^0(\varphi)$, by

$$ \Psi^0(\varphi) = \sum_{i=0}^{2N-1} N_i(\varphi) \Psi^0(\varphi_i), \qquad (107) $$

where $N_i(\varphi)$ are the interpolation functions. Now, we denote by $\tilde{N}_i$ the solution of the Fokker-Planck problem (66) associated with the initial condition $\tilde{N}_i(\underline{x}_0, \varphi) = N_i(\varphi)$. This solution is obtained using the implicit numerical strategy just proposed. Taking into account the linear and homogeneous character of the Fokker-Planck equation as well as the periodicity condition

$$ \Psi(\underline{x}_0, \varphi) = \Psi(\underline{x}(\underline{x}_0, t = T), \varphi), \qquad (108) $$

(T is the period related to the closed trajectory considered), we obtain

$$ \sum_{i=0}^{2N-1} N_i(\varphi) \Psi^0(\varphi_i) = \sum_{i=0}^{2N-1} \tilde{N}_i(\underline{x}(\underline{x}_0, t = T), \varphi) \Psi^0(\varphi_i), \qquad (109) $$

where $\tilde{N}_i(\underline{x}(\underline{x}_0, t = T), \varphi)$ is the solution after a turn along the closed streamline, related to the initial condition $N_i(\varphi)$, i.e. $\tilde{N}_i(\underline{x}_0, \varphi) = N_i(\varphi)$.
If we impose this equality at each angle $\varphi_j$, and taking into account that $N_i(\varphi_j) = \delta_{ij}$ ($\delta_{ij}$ is the Kroenecker's delta), we obtain

$$ \Psi^0(\varphi_j) = \sum_{i=0}^{2N-1} \tilde{N}_i(x(x_0, t=T), \varphi_j) \Psi^0(\varphi_i), \qquad (110) $$

which can be written in the matrix form

$$ \Psi = \underline{D}\Psi, \qquad (111) $$

where the i-row of $\Psi$ contains the value $\Psi^0(\varphi_i)$ and the j-component of $\underline{D}$ results $\tilde{N}_j(x(x_0, t=T), \varphi_i)$. The solution $\Psi^0(\varphi_i)$, $\forall i$, is obtained solving the linear system

$$ (\underline{D} - \underline{I})\Psi = 0, \qquad (112) $$

with the normality condition

$$ \int_{0}^{2\pi} \Psi^{0}(\varphi) \, d\varphi = 1, \qquad (113) $$

whose discrete form results

$$ \sum_{i=0}^{2N-1} \Psi^0(\varphi_i) h_{\varphi} = 1. \qquad (114) $$

# 5. Continuous approximation of the orientation distribution: numerical results and discussion

In this section, we present different results concerning simple flows in order to evaluate the capabilities and limitations of the numerical strategies just introduced.

## 5.1. Simple shear non-recirculating flow

We consider a simple shear and non recirculating flow in the limit case when diffusion vanishes, i.e. $D_r = 0$ in Eq. (66). The velocity field is defined again by

$$ \underline{v} = \begin{pmatrix} y \\ 0 \end{pmatrix}. \qquad (115) $$

In this case, we know that when the fibers have an infinite aspect ratio, the steady solution corresponds to the local alignment of the fibers with the flow. Thus, an orientation distribution concentrated around the angles $\varphi = 0$ and $\pi$, which corresponds with the discrete values $\varphi_0$ and $\varphi_N$, is expected. In fact, the steady solution obtained numerically using an explicit first-order upwind method, is
$$ \Psi^{\text{steady}} = \begin{pmatrix} \Psi(\varphi_0) \\ \Psi(\varphi_1) \\ \vdots \\ \Psi(\varphi_{N-1}) \\ \Psi(\varphi_N) \\ \Psi(\varphi_{N+1}) \\ \vdots \\ \Psi(\varphi_{2N-1}) \end{pmatrix} = \begin{pmatrix} \approx \frac{N}{2\pi} \\ \approx 0 \\ \vdots \\ \approx 0 \\ \approx \frac{N}{2\pi} \\ \approx 0 \\ \vdots \\ \approx 0 \end{pmatrix}, \qquad (116) $$

where “≈ 0” denotes values lower than $10^{-12}$.

These results are consistent because

$$ \int_{0}^{2\pi} \Psi(\varphi) \, d\varphi = \sum_{i=0}^{2N-1} \Psi(\varphi_i) h_{\varphi} = 1, \qquad (117) $$

and

$$ \Psi(\varphi_0) \Theta \left( \frac{1}{h_{\varphi}} \right), \qquad (118) $$

reproducing the Dirac’s delta behavior.

Now, fibers with a finite aspect ratio ($k < 1$) are considered. The maximum orientation probability remains around the angles $\varphi = 0$ and $\pi$, but its value decreases with respect to the case $k = 1$. In Table 1, we compare the result for different fiber aspect ratios, with $D_r = 0$ and $N = 10$. As expected, oscillations appear when centered schemes are used in the discretization of the advection term for diffusion coefficients approaching to zero, $D_r \to 0$.

<caption>Table 1</caption>

Numerical fiber orientation distribution computed using an explicit first-order method for different fiber aspect ratios ($D_r = 0, N = 10$)

<table><thead><tr><th>φ</th><th>Ψ(φ, k = 1)</th><th>Ψ(φ, k = 0.9)</th><th>Ψ(φ, k = 0.8)</th><th>Ψ(φ, k = 0)</th></tr></thead><tbody><tr><td>0</td><td>1.59</td><td>0.68</td><td>0.47</td><td>0.159</td></tr><tr><td>π/10</td><td>0.00</td><td>0.25</td><td>0.27</td><td>0.159</td></tr><tr><td>2π/10</td><td>0.00</td><td>0.09</td><td>0.12</td><td>0.159</td></tr><tr><td>3π/10</td><td>0.00</td><td>0.05</td><td>0.07</td><td>0.159</td></tr><tr><td>4π/10</td><td>0.00</td><td>0.04</td><td>0.05</td><td>0.159</td></tr><tr><td>5π/10</td><td>0.00</td><td>0.03</td><td>0.05</td><td>0.159</td></tr><tr><td>6π/10</td><td>0.00</td><td>0.04</td><td>0.05</td><td>0.159</td></tr><tr><td>7π/10</td><td>0.00</td><td>0.05</td><td>0.07</td><td>0.159</td></tr><tr><td>8π/10</td><td>0.00</td><td>0.09</td><td>0.12</td><td>0.159</td></tr><tr><td>9π/10</td><td>0.00</td><td>0.25</td><td>0.27</td><td>0.159</td></tr><tr><td>π</td><td>1.59</td><td>0.68</td><td>0.47</td><td>0.159</td></tr></tbody></table>
The last column of Table 1 corresponds with an isotropic fluid (fibers with an aspect ratio equal to one, $k = 0$). This isotropic state can be obtained also considering fibers with any aspect ratio and very high diffusion coefficients. The solutions for $k = 1$ and $k = 0$ are in perfect agreement with the corresponding exact solution, whereas the solutions for intermediate values of $k$ exhibit significant numerical diffusion. In order to evaluate the accuracy of the numerical solutions for $0 < k < 1$ we can compute a good approximation of the exact solution, which will be used as reference solution, using the ergodic principle, as described in Section 3.1.

Now, if we compare the solution for fibers with $k = 0.8$ obtained using the upwinding technique (described in Section 4) with $N = 10$ and $D_r = 0$, whose results have been grouped in Table 1, with the reference one, then the error results $||\Psi^{upw}(\varphi_i) - \Psi^{ref}(\varphi_i)||_{L^2} \approx 0.065$ with $||\Psi^{ref}(\varphi_i)||_{L^2} \approx 0.9$, which represents an error of 7%. Thus, we can conclude that the upwinding technique can introduce a significant numerical diffusion, unsurprisingly, in the limit case when diffusion vanishes (in spite of the good results obtained for $k = 0$ and 1) due to the first-order accuracy stabilization here considered.

## 5.2. Shear recirculating flow

In this case we consider the flow defined by

$$ \underline{v} = \begin{pmatrix} u \\ v \end{pmatrix} = \begin{pmatrix} -y\sqrt{x^2 + y^2} \\ x\sqrt{x^2 + y^2} \end{pmatrix}, \qquad (119) $$

and the limit case where the diffusion effects are neglected ($D_r = 0$) and the fibers have a quasi-infinite aspect ratio ($k = 1$).

We are looking for the steady solution at point $x = 0, y = 1$; that as we have proved in our former work [10] corresponds to the local alignment of the fiber with the flow, i.e. $\Psi((0, 1), \varphi) = \delta(\varphi - 0)$. However, using the numerical strategies proposed in Section 4 (explicit first-order upwind schema), the solution obtained for $N = 10$ (Table 2) does not correspond to the exact solution. A very slow convergence is obtained by changing the time step $\Delta t$ and the value of $N$. Using a finer discretization in the angular direction ($N = 100$), the computed solution $\Psi(\varphi)$ is depicted in Fig. 4.

<table><caption>Table 2
Numerical versus exact orientation distribution solution ($k = 1, D_r = 0, N = 10$)</caption><thead><tr><th>φ</th><th>Ψ<sub>num</sub>(φ)</th><th>Ψ<sub>exact</sub>(φ)</th></tr></thead><tbody><tr><td>0</td><td>0.26</td><td>1.59</td></tr><tr><td>π/10</td><td>0.16</td><td>0.00</td></tr><tr><td>2π/10</td><td>0.10</td><td>0.00</td></tr><tr><td>3π/10</td><td>0.07</td><td>0.00</td></tr><tr><td>4π/10</td><td>0.06</td><td>0.00</td></tr><tr><td>5π/10</td><td>0.06</td><td>0.00</td></tr><tr><td>6π/10</td><td>0.09</td><td>0.00</td></tr><tr><td>7π/10</td><td>0.17</td><td>0.00</td></tr><tr><td>8π/10</td><td>0.28</td><td>0.00</td></tr><tr><td>9π/10</td><td>0.33</td><td>0.00</td></tr><tr><td>π</td><td>0.26</td><td>1.59</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Numerical solution of the fiber orientation distribution in a shear recirculating flow at point (0, 1) computed using an explicit first-order upwind method ($k = 1, D_r = 0$).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Numerical solution of the fiber orientation distribution in a shear recirculating flow at point (0, 1) computed using an explicit first-order upwind method, for fibers with $k = 0.6$ and from top to bottom $D_r = 0.0, 0.18, 0.9$ and $1.8$.</figcaption></figure>
We can notice that the numerical discretization introduces a significant numerical diffusion as well as a slight displacement of the angles with a maximum orientation probability in the upstream angular convection direction. Moreover, the numerical solution converges slowly, which difficulties an accurate resolution. The only difference between both flows (the non recirculating and the recirculating one) is that in the first one the angular velocity vanishes at the angle related to the flow direction, which corresponds with the steady fiber orientation; whereas in the second one, the rotation velocity of the fibers is not zero anywhere. The use of the alternative implicit strategy described in Section 4.4 does not improve significantly these numerical results for the limit case of $D_r = 0$.

To evaluate the proposed strategy when diffusion effects become higher, we compute the steady solution at point (0, 1) related to fibers with $k = 0.6$ and different values of the parameter $D_r$. In this case, the rotation velocity of the fibers located at the indicated point results from the Jeffery's equation

$$ \dot{\varphi} = 1.8 - 0.6 \cos^2(\varphi). \qquad (120) $$

The maximum value of fiber rotation velocity is obtained for $\varphi = \pi/2$ and $\varphi = 3\pi/2$, whose value results $\dot{\varphi} = 1.8$. Thus, we propose to check our numerical strategy for $D_r = 0, 0.18, 0.9$ and $1.8$. Fig. 5 depicts the steady solutions computed at point (0, 1) for these four values of $D_r$ just indicated, using $N = 40$. We can notice that the introduction of a diffusive term and fibers with finite aspect ratio tends to smooth the fiber orientation distribution as well as to move the angles with a maximum orientation probability in the upstream angular convection direction. This fact justifies that the numerical diffusion induced by the discretization scheme in the limit case when diffusion vanishes, $D_r \to 0$, moves the fiber orientation distribution precisely in the upstream direction.

The error (using the $L^2$-norm) between the numerical solution and the reference one (computed like in Section 3.3) will depend on the mesh size $N$ and on the values of $D_r$ and $k$ ($E(N, D_r, k) = \|\Psi^{upw}(\varphi_i) - \Psi^{ref}(\varphi_i)\|_{L^2}$). Thus, we obtain the following values for the error $E(N = 40, D_r = 0.18, k = 0.6) = 0.1$, $E(N = 40, D_r = 0.9, k = 0.6) = 0.015$ and $E(N = 40, D_r = 1.8, k = 0.6) = 0.0068$, from which we can conclude that, as expected, the accuracy of the numerical solution increases as the ratio between the advection and diffusion terms decreases.

## 6. Conclusions

In this work, we have proposed firstly an accurate numerical strategy to compute steady solutions of the advection dominated Fokker-Planck equation in steady recirculating flows. This strategy combines a particle method with a more original treatment of the solution periodicity along the closed streamlines imposed by the steady regime. Moreover, this technique seems well adapted to treat also slight diffusion effects. The treatment of 3D, non-linear and slight diffusion models is a work in progress.

In the second part of this work, we have presented some discretization techniques making use of continuous approximations of the orientation distribution function, able to compute steady solutions of Fokker-Planck equations in steady recirculating flows where the diffusion coefficient lies in a large interval. The accuracy of those numerical strategies increases, as expected, as diffusion effects become dominant.
References

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

[2] E.J. Hinch, L.G. Leal, Constitutive equations in suspension mechanics. Part I, J. Fluid Mech. 71 (1975) 481-495.

[3] E.J. Hinch, L.G. Leal, Constitutive equations in suspension mechanics. Part II, J. Fluid Mech. 76 (1976) 187-208.

[4] F. Meslin, Propriétés Rhéologiques des Composites Fibres Courtes à l'Etat Fondu, Ph.D. thesis, Ecole Normale Superieure de Cachan, 1997.

[5] C. Givler, M.J. Crochet, R.B. Pipes, Numerically Predicted Fiber Orientation in Dilute Suspensions, NUMIFORM, Pineridge Press, 1982, pp. 559-575.

[6] J. Rosenberg, M. Denn, R. Keunings, Simulation of non-recirculating flows of dilute fiber suspensions, J. Non-Newtonian Fluid Mech. 37 (1990) 317-345.

[7] G. Ausias, Etude de l'Extrusion de Tubes en Polymères Thermoplastiques Chargés de Fibres Courtes, Ph.D. thesis, E.N.S.M.P, 1991.

[8] M.C. Altan, S.I. Guçeri, R.B. Pipes, Anisotropic channel flow of fiber suspensions, J. Non-Newtonian Fluid Mech. 42 (1992) 65-83.

[9] K. Chiba, K. Nakamura, Numerical solution of fiber suspensions flow through a complex channel, J. Non-Newtonian Fluid Mech. 78 (1998) 167-185.

[10] A. Poitou, F. Chinesta, R. Torres, Numerical simulation of the steady recirculating flows of fibers suspensions, J. Non-Newtonian Fluid Mech. 90 (2000) 65-80.

[11] H. Henry de Frahan, V. Verleye, F. Dupret, M.J. Crochet, Numerical prediction of fiber orientation in injection molding, Polym. Eng. Sci. 32 (4) (1992) 254-260.

[12] R.S. Bay, C.L. Tucker III, Fiber orientation in simple injection moldings. Part I. Theory and numerical methods, Polym. Composites 13 (4) (1992) 317-330.

[13] S.T. Chung, T.H. Kwon, Numerical simulation of fiber orientation in injection moldings of short-fiber-reinforced thermoplastics, Polym. Eng. Sci. 35 (7) (1995) 604-618.

[14] J. Azaiez, R. Guénette, A. Ait-Kadi, Investigation of the abrupt contraction flow of fiber suspensions in polymeric fluids, J. Non-Newtonian Fluid Mech. 73 (1997) 289-316.

[15] F. Chinesta, A. Poitou, R. Torres, Numerical prediction of the fiber orientation in steady flows, Revue Européenne des Eléments Finis 8 (4) (1999) 355-374.

[16] B. Souloumiac, Etude Rheologique, Modélisation et Simulation Numérique de l'Ecoulement des Thermoplastiques Chargés de Fibres Courtes, Ph.D. thesis, Ecole National Supérieure des Mines de Paris, 1996.

[17] F. Chinesta, A. Poitou, R. Torres, A semi-Lagrangian strategy to predict the fiber orientation in steady flows of reinforced thermoplastics, Comput. Meth. Appl. Mech. Eng. 189 (2000) 233-247.

[18] K. Walters, Overview of macroscopic viscoelastic flows, in: S.A. Lodge, M. Renardy, J.A. Nohel (Eds.), Viscoelasticity and Rheology, Academic Press, 1985, pp. 47-79.

[19] R.E. Evans, K. Walters, Further remarks on the lip-vortex mechanism of vortex enhancement in planar-contraction flows, J. Non-Newtonian Fluid Mech. 32 (1989) 95-105.

[20] P. Townsend, K. Walters, Expansion flows of non-Newtonians liquids, Chem. Eng. Sci. 49 (5) (1993) 749-763.

[21] F. Chinesta, Modélisation Numérique en Mise en Forme de Polymères et Ceramiques: Differents Probl Isquo;emesde Transport, Habilitation à Diriger des Recherches, Université Paris 6, 1999.

[22] F. Chinesta, G. Chaidron, On the steady solution of linear advection problems in steady recirculating flows, J. Non-Newtonian Fluid Mech. 98 (2001) 65-80.

[23] G. Chaidron, F. Chinesta, On the steady solution of non-linear advection equations in steady recirculating flows, Comput. Meth. Appl. Mech. Eng. 191 (2002) 1159-1172.

[24] A.J. Szeri, L.G. Leal, A new computational method for the solution of flow problems of microstructured fluids. Part 1. Theory, J. Fluid Mech. 242 (1992) 549-576.

[25] A.J. Szeri, L.G. Leal, A new computational method for the solution of flow problems of microstructured fluids. Part 2. Inhomogeneous shear flow of a suspension, J. Fluid Mechanics 262 (1994) 171-204.

[26] A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid Mech. 57 (4) (1973) 785-796.

[27] H.C. Ottinger, Stochastic Processes in Polymeric Fluids, Springer, Berlin, 1996.

[28] P. Degond, F.J. Mustieles, A deterministic approximation of diffusion equations using particles, Siam J. Sci. Stat. Comput. 11 (2) (1990) 293-310.
[29] C.V. Chaubal, A. Srinivasan, O. Egecioglu, L.G. Leal, Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part 1. Method, J. Non-Newtonian Fluid Mech. 70 (1997) 125–154.

[30] L.G. Leal, E.J. Hinch, The effect of weak Brownian rotations on particles in shear flow, J. Fluid Mech. 46 (4) (1971) 685–703.

[31] F. Chinesta, K. Chiba, New numerical strategies for solving kinetic theory problems: rheological applications, in: Proceedings of the 6th International ESAFORM Conference on Material Forming, Salerno, Italy, 2003.

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

[33] V. Verleye, F. Dupret, Prediction of fibre orientation in complex injection moulded parts, in: Developments in Non-Newtonian Flows, AMD-vol. 175, ASME, New York, 1993, pp. 139–163.

[34] J.P. Boyd, Chebyshev and Fourier Spectral Methods, Springer-Verlag, Berlin, 1989.