A Particle Strategy for Solving the Fokker-Planck Equation Governing the Fiber Orientation Distribution in Steady Recirculating Flows Involving Short Fiber Suspensions

In this work we have analyzed the application of a deterministic approximation of the diﬀusion term in the Fokker-Planck equation using smooth particles for computing its steady solution in a steady recirculating ﬂow. The main idea of this approach lies in the introduction of the Fokker-Planck diﬀusion term into the ad-vection one, which allows to proceed in a Lagrangian deterministic manner without a mesh support requirement


Introduction
As indicated by Chaubal et al. [CHA97], "complex fluid" is the term commonly used to describe a wide class of liquid-like materials, in which the relaxation time towards the equilibrium state occurs sufficiently slowly that significant changes in the microstructural configuration, and thus in the macroscopic properties, can be induced by the flow.Viscoelastic fluids or short fiber suspensions may be considered as examples of complex fluids.The Fokker-Planck formalism is a commonly used description of kinetic theory problems, for describing the evolution of the configuration distribution function.This function represents the probability of finding the microstructural element in a particular configuration.
In the case of a short fiber suspension, the configuration distribution function (also known as orientation distribution function) gives the probability of finding the fiber in a given direction.Obviously, this function depends on the physical coordinates (space and time) as well as on the configuration coordinates, that taking into account the rigid character of the fibers, are defined on the surface of the unit sphere.Thus, we can write Ψ (x, t, p), where x defines the position of the fiber center of mass, t the time and p the unit vector defining the fiber orientation.The evolution of the distribution function is given by the Fokker-Planck equation where d/dt represents the material derivative, D r is a diffusion coefficient and ṗ 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's equation where Ω and D are the vorticity and the strain rate tensors respectively, associated with the fluid flow undisturbed by the presence of the fiber, and k is a scalar which depends on the fiber aspect ratio λ (ratio between the fiber length and the fiber diameter) Although this work focuses in the flow of short fiber suspensions, the numerical procedures here developed could be applied successfully to other fluids whose microstructure is described by similar kinetic theory models.This is the case for example of some viscoelastic fluids.In the FENE model [BIR87], the probability distribution function depends on the physical coordinates and on the conformation coordinates that in this case are defined by the end-to-end molecule vector q.The equation governing its evolution results: where the evolution of q can be evaluated from with the spring force F given by being b the maximum molecule length.The FENE-P model is obtained by averaging the square norm of the molecule length in the expression of the spring force [BIR80].
In these cases the Fokker-Planck equations remain linear, even if the advection fields depend in a non-linear manner on the configuration coordinates (see for instance Eqs.(1.2) and (1.5)).
Another more complex example consists in the model describing the behavior of Liquid Crystals Polymers (LCP) due to Doi.This model has been treated using a particle technique in [CHA97].In this case the microstructural element is a rigid axis-symmetric rod of infinite aspect ratio.Under some simplifying assumptions the orientation distribution (also defined on the surface of the unit sphere) evolves according to where the advection field ṗ is given by As the potential U depends on the second moment of the distribution function p ⊗ p defined by and the Fokker-Planck equation results in this case non-linear.Many of the experimental and industrial flows show recirculating areas or recirculate themselves.For example, many rheometric devices involve this type of flows, whose steady and recirculating character introduce some additional difficulties in their numerical simulation.Actually, the Fokker-Planck equation which defines an advection problem in physical coordinates, 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 [CHI03] the discretisation of the advection dominated Fokker-Planck equation, governing the fiber orientation in short fiber suspension flows, was carried out using a particle technique, where the diffusion term was modelled from random motions.It was pointed out that the number of fibers required in this stochastic simulation to describe the fiber distribution increases significantly with the diffusion coefficient D r .Thus, it was argued that for practical applications the use of the particle method in the framework of a stochastic simulation, is restricted to very slight diffusion effects.When the diffusion becomes dominant, continuous approximations using a fixed mesh seem to be suitable, but in this case accurate stabilizations are required for dealing with small diffusion effects, and a lack of accuracy is noticed in the treatment of the advection dominated case.
Chaubal et al. [CHA97] propose the use of the SPH (Smooth Particle Hydrodynamics) to solve the dynamics of a liquid crystalline polymer (LCP) using the Doi's model.Thus, the diffusion term is treated in a determinist manner, and high accurate results were obtained due to the meshless and Lagrangian character of the SPH technique considered by the authors.
In this work we examine the application of this kind of techniques for solving the steady Fokker-Planck equation in steady recirculating flows, but we will limit our attention to the 2D linear Fokker-Planck equations.

A Particle Discretisation
In this section we will consider the simplest linear form of the Fokker-Planck equation (1.1), that in the 2D case results where Ψ (x, t, ϕ), φ = φ(x, t, ϕ) and D r is assumed constant.Due to the steady character of the flow kinematics we can remove the temporal variable in the expression of φ.A general approximation of the steady probability distribution function at point x 0 can be written as Due to the linearity and homogeneity of Eq. (2.1) we can compute its solution for each function F i (ϕ) along the pathline related to the point x 0 .These solutions are denoted by Ψ i (x(x 0 , t 0 ; t), ϕ), where the notation x(x 0 , t 0 ; t) refers to the position at time t of a particle located at point x 0 at time t 0 .The trajectory is defined by the equation Thus, the general solution of Eq. (2.1) results The particle method lies in taking the Dirac masses in Eq. (2.2).Thus, if we consider the N directions defined by where α i (x 0 ) represents the "mass" of the particle aligned in the ϕ i direction.

The Advection Equation
In the case of neglecting all the diffusion effects, the solution Ψ i (x(x 0 , t 0 ; t), ϕ) results where ϕ(x 0 , t 0 , ϕ i ; t) denotes the orientation at time t of a fiber located at time t 0 in x 0 and whose orientation was defined by ϕ i .Obviously, the spatial location of a fiber refers to the position of its center of mass, that at time t is given by Eq. (2.3).
In the pure advection case ϕ(x 0 , t 0 , ϕ i ; t) becomes Thus, Eqs.(2.3) and (2.7) describe the position and orientation of each fiber along the flow trajectory.Several discretisations of these equations exist, being the simplest one the backward Newton method, that consider the time interval [t 0 , t] divided into M intervals [t m , t m+1 ] of length Δt such that M Δt = t − t 0 , and the fibers updating given by (2.8) where (2.9)

Steady Recirculating Flows
In the case of a steady recirculating flow the particles come back to the departure position after a time T , that corresponds with the period of the considered trajectory.Thus, we can write x(x 0 , t 0 ; T ) = x(x 0 , t 0 ; t 0 ) = x 0 (2.10)However, the final orientation will be, in general, different to the initial one, i.e. ϕ(x 0 , t 0 , ϕ i ; T ) = ϕ i , ∀i (2.11) On the other hand, the steady solution of the probability distribution at point x 0 , Ψ (x 0 , ϕ), requires its periodicity along the closed streamline (2.12) but this expression, (Eq.(2.12)), cannot be used in its present form because the Dirac masses are concentrated in different angles at t 0 and T .
To use this expression we need to transfer the masses concentrated in the directions ϕ(x 0 , t 0 , ϕ i ; T ) to the initial ones ϕ i .In [CHI03] this transfer is performed from each angle ϕ(x 0 , t 0 , ϕ i ; T ) towards the neighbor directions, being the mass transferred to each neighbor direction proportional to the distance between them.Thus, we can finally write (see [CHI03] for more details) where β ij depends in the considered probability transfer from ϕ(x 0 , t 0 , ϕ i ; T ) to ϕ j .If a linear transfer is used, then it results The periodicity condition becomes in this case where δ ij is the unit tensor.The previous equation implies j=N j=1 that with the normality condition allows to compute the N coefficients α i (x 0 ), and then, the steady solution of the probability distribution at point x 0

Introducing Diffusion Effects Using a Deterministic Approximation
When diffusion effects are retained in the kinetic model the application of a particle technique becomes delicate as we will describe in this section.The first possibility for introducing of the diffusion effects lies in the use of a stochastic technique [CHI03].Thus, a great amount of particles are considered at point x 0 aligned in each direction ϕ i .Then, each fiber is subjected to three actions during its movement along its closed trajectory: (i) the advection of its center of mass, (ii) the rotation induced by the term ṗ and (iii) the rotation associated with the diffusion that is modelled from a random motion.In spite of the simplicity of its computational implementation and the very accurate results obtained for slight diffusion coefficients, the great number of particles required when the diffusion coefficient increases makes this strategy useless for treating the problems encountered in practical applications.
To circumvent these computational drawbacks, we transform the linear Fokker-Planck advection-diffusion equation in a pure advection problem.

Ψ
(2.20) Now, Eq. (2.21) could be used for computing the different solutions associated with the Dirac's distributions δ(ϕ − ϕ i ), ∀i, but two difficulties appear: (i) the derivative and the ratio of linear combinations of Dirac masses cannot be defined in a proper way, and (ii) due to the distribution function derivative required to compute the modified advection field φ * in Eq. (2.20), a coupling between the different fibers takes place, and consequently, in this case, the evolution of each fiber cannot be computed independently of the other ones.
In order to circumvent the first difficulty we introduce a smoothed approximation of the Dirac mass given by To adapt its sharpness we modify this function in the following manner which verifies also the normality condition.Fig. 2.1 depicts this function for three values of the parameter.

Representation of the cut-off function for different values of the parameter
In this form, the gradients and ratios can be computed in a proper way.On the other hand to avoid the second limitation two simple possibilities exist: (i) a direct procedure that compute N times, the evolution history of N particles along the closed trajectory, and (ii) an iteration algorithm that computes the evolution of each particle along the closed trajectory independently of the other ones due to the fact that the coupling term is evaluated from the solution at the previous iteration.In the following sections we describe both procedures.
A Direct Procedure.We write the searched solution at point x 0 in a matrix form where the unknown coefficients α i (x 0 ), contained in the vector α(x 0 ), must be determined in order to satisfy the periodicity condition imposed by both the steady character of the searched solution and the steady and recirculating character of the flow.We define the following N alpha-vectors: where ν and μ can be chosen arbitrarily with the only restriction of verifying the normality condition, i.e. (N −1)ν+μ = 1.We can notice that the uncoupled procedure, applied when diffusion effects are neglected, uses ν = 0 and μ = 1.Introducing this notation, Eq. (2.25) can be rewritten as Now, we can compute the evolution along a closed trajectory of each solution which in fact requires the tracking of N particles.For this purpose, we need to integrate the evolution of the orientation for each fiber (represented by the subscript i) in each problem (noted by the superscript k).If we consider, for a sake of simplicity, the backward Newton method again, it results After a complete turn of period T we obtain Now the weights associated to each orientation ϕ(x 0 , t 0 , ϕ k i ; T ), β k α k i (x 0 ), must be transferred to the angles ϕ i to enforce the periodicity.This operation can be accurately carried out by using the algorithm proposed in [CHI03] as previously described.
This algorithm offers the steady solution of the problem from the resolution of the N problems, requiring each one of them, the integration of the orientation evolution of N fibers along the closed trajectory.In this form the algorithm results of order N 2 .
A Fixed Point Iteration Algorithm.In this case we start from the solution obtained assuming D r = 0.This solution can be written as Now, the iteration k results where in this case the superscript k ≥ 1 refers to the iteration.
This algorithm results of order N but the convergence is not ensured.Some improvements can be introduced, as for example the introduction of the diffusion coefficient step by step.

Numerical Examples
In this section two numerical examples involving steady recirculating flows will be considered.The first one is defined by the following kinematics: The steady solution of the probability distribution is searched at point (0, 1).For this problem the exact solution can be computed due to the solution symmetry as described in [CHI03].Thus, for a diffusion coefficient D r = 0.2 and fibers with an aspect ratio related to k = 0.6, the fixed point algorithm proposed in the section 2.3 with = 0.5 and N = 72 converges in around 10 iterations to the reference solution as depicted in Fig. 3.2.Now, the direct method described in the section 2.3 is applied.Obviously, the numerical solution accuracy will depends on the μ coefficient (for μ ∈ the extact and the computed solutions are in very good agreement).Figs.3.3 and 3.4 depicts the fiber distribution related to μ = 1.2/N and μ = 10/N respectively as well as the exact solution (dashed line).The case of μ = 0 (where N − 1 fibers with equal weights are considered) is depicted in Fig. 3.5.The solution accuracy can be improved by increasing the number of particles N involved in the simulation.A first order convergence is noticed when the number of particles N is increased.In these examples we don't take μ in the interval [3,5] because in that case the numerical and reference solutions are completely superposed.Fig. 3.6 depicts the fiber distribution in some points on the closed streamline related to point (0, 1) computed using the following model parameters: k = 0.8, = 0.5, D r = 0.2 and N = 72.The fiber orientation is then defined on the unit circle, where the orientation of each fiber is depicted by a small For this purpose we define the new advection field ṗ * ṗ * = ṗ − D r standard linear advection-diffusion Fokker-Planck equation (1.1).In the 2D case developed in the previous section, both equations result in the following purely-advection problem φ * = φ − D r ∂Ψ ∂ϕ