Etienne Prulière

Amine Ammar
email: amine.ammar@ujf-grenoble.fr

Nadia El Kissi

Etienne Pruliere

Francisco Chinesta
email: francisco.chinesta@paris.ensam.fr

Recirculating Flows Involving Short Fiber Suspensions: Numerical Difficulties and Efficient Advanced Micro-Macro Solvers

whether they are published or not. The documents may come   L'archive ouverte pluridisciplinaire

1 Introduction

Recirculating flows of complex fluids

Understanding recirculating flows of complex fluids is of importance in many areas of polymer processing such as extrusion, injection or mixing. However, in such complex flows it is difficult to measure the detailed kinematics and stresses that develop during flow. It is thus necessary to predict them, by making good measurements of basic fluid properties and using them in the appropriate constitutive equations. This means taking into account the specific properties of complex fluids. Indeed, viscosity and shear thinning, elasticity and normal stresses, elongation and strain hardening as well as relaxation phenomena induce several unusual phenomena, and recirculating flows are one of them.

They have been observed in many flow configurations. In the case of Newtonian and inelastic fluids, their occurrence and evolution are mainly related to the inertial properties of the flow. Thus, the Rayleigh-Taylor instability occurs during flow between two coaxial cylinders [START_REF] Taylor | Stability of a viscous liquid contained between two rotating cylinders[END_REF], the external one being fixed while the internal one rotates. Due to centrifugal forces, the fluid is pushed from the inner cylinder towards the external one. Beyond a certain threshold, these forces cause recirculating flows in the form of rollers developing in the gap between the two cylinders. This is the well known inertial instability of Rayleigh-Taylor. Notice that the term instability in this case indicates that a secondary effect is superimposed on the main flow in controlling the kinematics. In fact, the recirculating areas are stable in time and space. For a number of years now, many authors have studied the Rayleigh-Taylor inertial instability. Theoretical studies have analyzed the role of elasticity and showed that this property could stabilize flow, delaying, if not eliminating, the occurrence of the inertial instability [START_REF] Muller | A purely elastic transition in Taylor-Couette flow[END_REF].

It is also possible to observe recirculating flows independently of inertia. In the case of complex fluids for example, they result from properties such as viscosity, but also elasticity, elongation or relaxation phenomena. Such recirculating areas have been observed for example in flows through non-circular capillaries, typically of square or ellipsoidal section [START_REF] Oldroyd | Some steady flows of the general elastic-viscous liquid[END_REF] [45] [START_REF] Giesekus | Zur Stabilitt von Strmungen viskoelastischer Flssigkeiten[END_REF] [START_REF] Dodson | Non-Newtonian flow in pipes of non-circular crosssection[END_REF]. They have also been observed when complex fluids flow through abrupt or shaped contractions, the fluid being pushed from a large-diameter reservoir, through a contraction and then into a capillary of smaller diameter. The high ratio between the capillary and reservoir diameters is defined as the contraction ratio. This flow configuration was the subject of thorough numerical as well as experimental studies. The reasons for this are numerous: From an academic point of view, this type of flow is of interest as an appropriate test problem. They also allow the elongational properties of complex fluids, which are generally difficult to measure, to be evaluated [START_REF] Cogswell | Converging flow of polymer melts in extrusion dies[END_REF] [START_REF] Binding | An approximate analysis for contraction and converging flows[END_REF]. From an industrial point of view, contraction flows are of great importance in modelling polymer processing involving many contraction and expansion sections, such as extrusion. The occurrence, evolution, and geometrical characteristics of the recirculating areas developing upstream of the contraction have a great impact in controlling the end-use properties of the material. Indeed, the recirculating fluid undergoes a thermo-mechanical history which is different from that of the flowing fluid. This can induce heterogeneities in the final properties of the material, thus affecting its mechanical characteristics and degradation properties. It is thus essential to understand and control the origin and development of such recirculating flows, the aim being to predict and optimize the kinematics with a view to eliminating secondary flows and regions of high stress.

If Newtonian or shear-thinning inelastic fluids are considered, it has been shown experimentally and by numerical simulation [START_REF] Giesekus | Verschiedene Phnomene in Strmungen viskoelastischer Flssigkeiten durch Dsen[END_REF] [8] [START_REF] Kim-E | The roles of inertia and shear-thinning in flow of an inelastic liquid through an axisymmetric sudden contractionJ[END_REF] [START_REF] Boger | Viscoelastic flows through contractions[END_REF] that inertia is responsible for the appearance of a recirculating area localized in the corner of the contraction.

The scenario is quite different when viscoelastic fluids are considered. The occurrence, evolution and shape of the recirculating zone are controlled completely by the viscoelastic properties. At low flow rates, the fluid behaves as a Newtonian fluid, with small vortices developing due to inertia. As the flow rate increases, the shape and size of the recirculating area are completely different from what is observed with Newtonian liquids. Two main phenomena are usually observed and are a good demonstration of viscoelastic behavior vortex enhancement the occurrence of a second recirculating zone, developing at the entrance of the contraction, called the lip vortex. Compared to Newtonian behavior where only small recirculating areas are observed, the two vortex areas coexist at low flow rates. As flow rates are increased, the lip vortex grows significantly, and progressively invades the corner vortex, leading to a single recirculating area in the upstream flow. The height of this recirculating area may reach several reservoir diameters, thus invading a large proportion of the upstream flow. Thorough quantitative experimental studies, relating basic rheological measurements to the occurrence and evolution of the recirculating area patterns [START_REF] Dennison | Flow instability in polymer melts: a review[END_REF] [START_REF] Petrie | Instabilities in polymer processing[END_REF], suggested that elasticity is responsible for vortex enhancement. Represented in terms of a Weissenberg number, this was used successfully to explain the variation in vortex size and in particular to correlate the growth of the recirculating area.

This situation is also encountered when flows involving short fiber suspensions are addressed. In this case the size of the vortex appearing in abrupt contractions depends on the suspension concentration and can be accounted without introducing inertia effects [START_REF] Lipscomb | The flow of fiber suspensions in complex geometries[END_REF].

Understanding recirculating flows thus appears to be of great importance in order to optimize polymer processes as popular as extrusion or injection. Actually many rheometric devices involve this type of flows. The ability to predict the characteristics of the observed recirculating areas is highly dependent on: accurate experimental observation made with adequately characterized material, appropriate constitutive equations reflecting all the complexity of the material and in particular the extensional properties of the fluid being studied. Good numerical methods, capable of simulating the flow patterns and accurately describing the kinematics as well as the level of stresses, are needed for this purpose.

Mechanical modelling of short fiber suspensions

Numerical modelling 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.

Thus, the flow model related to a short fiber suspension is given by the the following equations [START_REF] Azaiez | State-of-the-Art on numerical simulation of fiberreinforced thermoplastic forming processes[END_REF]:

• Equilibrium equation (neglecting mass and inertia terms):

Divσ = 0 (1)
where σ represents the stress tensor.

• Fluid incompressibility:

Divv = 0 ( 2
)
where v is the velocity field.

• The constitutive equation that, for dilute suspensions of high aspect ratio fibers, can be assumed as:

σ = -pI + 2µD + τ (3)
where p denotes the pressure, I the unit tensor, D the strain rate tensor (symmetric component of the gradient of velocity tensor) and τ the anisotropic viscous component of the stress tensor due to the presence of fibers.

The fluid microstructure can be defined from the fiber orientation. From now on the fibers will be considered ellipsoidal rigid particles. If ρ denotes the unit vector aligned in the fiber axis direction, its evolution, when the fiber is immersed in a Newtonian matrix whose kinematics is defined by Gradv, is given by the Jeffery's equation:

dρ dt = Ω ρ + k(D ρ -(D : (ρ ⊗ ρ))ρ) (4)
where k = λ 2 -1 λ 2 +1 , being λ the fiber aspect ratio (fiber length to fiber diameter ratio), and Ω denotes the vorticity tensor. The tensor product ⊗ of vectors a and b is defined by

(a ⊗ b) ij = a i b j .
Now, we can also define the orientation distribution function Ψ(ρ) given the probability of finding fibers aligned in the direction ρ. This distribution function, that verifies the normality condition Ψ(ρ)dρ = 1, allows the definition of different moments:

• Second order orientation tensor:

a = ρ ⊗ ρ Ψ(ρ) dρ (5)
• Fourth order orientation tensor:

A = ρ ⊗ ρ ⊗ ρ ⊗ ρ Ψ(ρ) dρ (6)
both integrals defined on the unit sphere.

Remark. In fact the distribution function Ψ depends on the physical coordinates x, the time t and the conformation coordinates ρ, i.e. Ψ(x, t, ρ). Some times, for the sake of simplicity in the notation, we omit its dependence on x and t.

The evolution of the orientation distribution function is governed by the Fokker-Planck equation:

dΨ dt + ∂ ∂ρ Ψ dρ dt = ∂ ∂ρ D r ∂Ψ ∂ρ (7)
where the advection field dρ dt is given by the Jeffery's equation ( 4), and dΨ dt refers to the material derivative, i.e.

dΨ dt = ∂Ψ ∂t + v • Grad x Ψ (8)
where

Grad x Ψ = ∂Ψ ∂x (9)
Using a spatial homogenization and a statistical averaging as well as other simplifying hypotheses, tensor τ in Eq. ( 3) results

τ = 2µN p (A : D) (10)
Instead of the use of the microscopic model defined by Eqs. ( 10), ( 6) and ( 7), macroscopic descriptions have been widely used because its lower computational requirements. Thus, taking the time derivative of Eq. ( 5) and introducing the Jeffery's equation ( 4) as well as the Fokker-Planck equation [START_REF] Azaiez | State-of-the-Art on numerical simulation of fiberreinforced thermoplastic forming processes[END_REF] it results:

da dt = Ω a -a Ω + k(D a + a D -2(A : D)) -2dD r a - I d (11)
where d = 2 and d = 3 in 2D and 3D respectively. However, the equation governing the evolution of the second order orientation tensor uses the fourth order one, the one governing the volution of the fourth order uses the sixth order one and so on. In conclusion, a closure relation expressing the fourth order orientation tensor as a function of the lower order orientation tensors must be introduced, being the most common relations the quadratic, linear, hybrid and natural [START_REF] Advani | Closure approximations for three-dimensional structure tensors[END_REF] [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF]. It is today widely recognized the significant impact that those closure relations can introduce in the computed solutions. Thus, a purely microscopic approach based on the solution of the Fokker-Planck equation seems to be an appealing choice for treating short fiber suspension models.

This paper focuses on the resolution of the steady recirculating flows involving a short fiber suspension whose description is accomplished in the kinetic theory framework from the introduction of the Fokker-Planck (FP) equation governing the fiber orientation distribution. In our knowledge, models coupling the flow equations with a description of the microstructure evolution based on the use of deterministic kinetic theory models, are extremely rare. In the context of polymer suspensions or entangled polymers, micro-macro techniques (coupling the flow with the microstructure evolution) are performed by considering stochastic techniques. These techniques are based in the equivalence between the deterministic Fokker-Planck equation and the associate Ito stochastic equation. Thus, instead the resolution of the FP equation, one can determine the trajectories related to a large number of realizations of the stochastic process, number that depends on the model and flow considered. The main reason for preferring stochastic simulations instead the resolution of the deterministic FP partial differential equation is that this last equation involves the distribution function which depends on the physical coordinates, the time, and the conformation coordinates related to the orientation. As this orientation is defined by a unit vector, the conformation space involves only two coordinates. Thus, for each point in the physical space and time, the orientation distribution is defined on the unit sphere. Standard finite element discretizations fail in this situation due to the highly dimensional character of the problem.

In this scenario, the first works concerning micro-macro simulations are related to the CONNFESSIT stochastic method [START_REF] Öttinger | Smart polymers in finite element calculation[END_REF], [START_REF] Keunings | Micro-macro methods for the multiscale simulation viscoelastic flow using molecular models of kinetic theory[END_REF] (see also the references therein). This approach was considered in [START_REF] Somasi | Brownian dynamics simulations of bead-rod and bead-spring chains: numerical algorithms and coarse-graining issues[END_REF] for treating MBS (multi bead spring) models. A similar technique was proposed in [START_REF] Venkiteswaran | A QMC approach for high dimensional Fokker-Planck equations modelling polymeric liquids[END_REF] in the context of MBS kinetic models, which introduces a change of variable and uses a Monte-Carlo technique for accounting the diffusion term. In these techniques a high number of particles are introduced in the simulation and a stochastic technique is used for accounting the Brownian effects. A multi-scale approach using deterministic particles for treating the advection and a different set of particles for accounting diffusion effects, which leads in fact to a multi-scale approach, was considered in [START_REF] Jendrejack | A method for multiscale simulation of flowing complex fluids[END_REF]. The same idea was used in the case of short fibre suspensions flows in [START_REF] Chinesta | On the solution of the Fokker-Planck equation in steady recirculating flows involving short fiber suspensions[END_REF]. In that work, the discretization of the advection dominated Fokker-Planck equation governing the fibre orientation, was carried out using a particle technique, where the diffusion term was modelled from random walks. It was pointed out that the number of fibres required in this stochastic simulation to describe the fibre distribution increases significantly with the diffusion coefficient. 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. Other deterministic particle approach, very close to that proposed in [START_REF] Chaubal | Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems[END_REF], was analyzed in [START_REF] Ammar | A particle strategy for solving the Fokker-Planck equation governing the fiber orientation distribution in steady recirculating flows involving short fiber suspensions[END_REF] using smooth particles, but it was noticed that the incidence of smoothing on the solution can be significant.

Thus, when the diffusion effects can vary in a large interval, continuous approximations using fixed or moving meshes seem to be a suitable alternative. In this case accurate stabilizations are required for dealing with small diffusion effects. However, as just pointed out, the Fokker-Planck equation is defined in a multidimensional space. Some attempts for solving the Fokker-Planck equation using a fixed mesh discretization exist [START_REF] Lozinski | A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model[END_REF] [START_REF] Chauviere | Simulation of dilute polymer solutions using a Fokker-Planck equation[END_REF]. Due to the multidimensional character of the problem, the linear systems obtained after conventional implicit or semi-implicit space-time discretizations result to be extremely large for a practical inversion. On the other hand, explicit discretizations, which do not require matrix inversions, have the constraint of too small time steps.

In this paper we consider a splitting technique that allows to decouple the advection problem defined in the physical space that is assumed corresponding to a steady recirculating flow, and the advection-diffusion one involving the conformation coordinates. The last equation is discretized using a stabilized finite element scheme, whereas two different schemes are analyzed for solving the advection equation in the physical space, one based on the discontinuous Galerkin method and the other one based on searching directly the steady solution by imposing the solution periodicity along the closed streamline. In addition, to alleviate the computing efforts, model reduction techniques (similar to the ones proposed in [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] and [START_REF] Ammar | A new family of solvers for some clases of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF]) will be applied and evaluated. The first reduction strategy is based on the use of a proper orthogonal decomposition within an adaptive scheme, whereas the second one reduces the the computing efforts by performing separated representations of the unknown fields.

Kinematic solver

As this paper focuses on the driven-cavity flow problem, we assume the flow model defined in the cavity domain Ω, being the velocity prescribed on its whole boundary Γ. Thus, the variational formulation related to the flow kinematics (Eqs. (1) and ( 2)) is given by: Ω σ : D * dΩ = 0 [START_REF] Bungartz | Sparse grids[END_REF] where ":" denotes the tensor product twice contracted, and D * represents the strain rate tensor related to a virtual velocity field v * which vanishes on Γ.

Introducing now the constitutive equation (3) it results:

Ω (-pI + 2µ(D + N p (a : D)) : D * dΩ = 0 (13)
The variational formulation related to the incompressibility condition results:

Ω Divv p * dΩ = 0 ( 14
)
where p * is a scalar weighting function. Eqs. ( 13) and ( 14) are then discretized by using the standard Galerkin finite element method. A P 2 -P 1 interpolation is considered for approximating the velocity and pressure fields, which verifies the well known LBB stability condition.

The driven-cavity problem considered in this paper defines a fully recirculating flow. We will prove later that in these recirculating flows some particular numerical strategies can be more accurate than standard finite elements. Those strategies are based on the integration along the closed trajectories, which can be accurately described in 2D using the stream function that we denote by Ξ (there exist other techniques able to compute eventual closed streamlines in the 3D case). Thus, in the 2D case the velocity field can be written, accounting the incompressibility, as

u = ∂Ξ ∂y v = -∂Ξ ∂x (15)
From the above expression the following elliptic problem results:

∆Ξ = ∂u ∂y - ∂v ∂x (16)
whose weighted residual form results

Ξ * ∆Ξ dΩ = Ξ * ∂u ∂y - ∂v ∂x dΩ (17)
According to the prescribed velocity on the domain boundary, one must search for the appropriate value of the stream function, or its normal derivative, to be imposed on the domain boundary.

Usually, in driven-cavity flows problems the component of the velocity vector normal to the boundaries vanishes, which corresponds to a null tangential stream function derivative on the boundary. One possibility to verify this last condition consists of taking a null value of the stream function on the whole domain boundary, which allows to write the following variation formulation:

-gradΞ * gradΞ dΩ = Ξ * ( ∂u ∂y - ∂v ∂x ) dΩ (18)
where Ξ, Ξ * ∈ H 1 0 , being H 1 0 the usual Sobolev's functional space.

3 Efficient solvers for advection problems defined in a steady recirculating flow Models governing the microstructure evolution using a macroscopic description (based on the use of the orientation tensors) or the microscopic one (based on the use of the orientation distribution function whose evolution is governed by the Fokker-Planck equation) are all of them defined by an advection equation in the physical coordinates.

In the present paper we use and compare different advanced solvers. The first one based on the use of a first order discontinuous Galerkin technique in the finite element framework, and the second one based on a particular application of the method of characteristics.

In this section we will summarize both strategies on a generic advection problem of the generic vector field a ∈ R M , given by:

da dt = A(x) a + B(x) (19)
whose steady solution is searched. In the above equation the material derivative is given by

da dt = ∂a ∂t + (v • grad)a (20)

A first order discontinuous Galerkin technique

Taking into account the velocity field incompressibility, Eq. ( 19) can be written as:

∂a ∂t + Div(v ⊗ a) = A(x) a + B(x) (21)
Now, we consider a partition of Ω in a set of non-overlapping cells Ω e such that Ω = ∪ e=Ne e=1 Ω e . Integrating Eq. ( 21) in each cell, and using the divergence theorem, it results:

Ωe ∂a ∂t dΩ + Γe a (v • n) dΓ = Ωe (Aa + B) dΩ (22)
where n represents the unit outward vector defined on the cell boundary. Eq. ( 22) is equivalent to the weighted residual Galerkin formulation when one considers a constant approximation of both trial and test functions into each cell (being that approximation discontinuous across the cell boundaries). In this case, the previous equation can be simplified to:

∂a e ∂t |Ω e | + Γ + e a (v • n) dΓ + Γ - e a (v • n) dΓ = Ωe (Aa e + B) dΩ (23)
where Γ + e and Γ - e are the outflow and inflow boundaries of cell Ω e respectively. As the approximation of a is not defined on the cells boundaries, we assume that a e (x ∈ Γ + e ) = a e and a e (x ∈ Γ - e ) = a e -, where a e -denotes the value of a in the upstream neighbor cell with respect to Γ - e . Moreover, the flow incompressibility results in

Γ + e v • n dΓ + Γ - e v • n dΓ = 0 (24)
Thus, we can finally write:

Γ + e a (v • n) dΓ = a e Γ + e v • n dΓ = -a e Γ - e v • n dΓ (25)
Now, if we assume the cells boundary composed of n ed edges, then it can be stated

Γ - e a (v • n) dΓ = i=n ed i=1 a e i Γ i e γ i e (v • n) dΓ (26)
where a e i denotes the value of a in the neighbor cell whose common edge with Ω e is Γ i e . In order to restrict the sum to the inflow contributions we have introduced the parameter γ i e which is defined by:

γ i e = 1 if v • n < 0 0 if v • n > 0 ( 27
)
From all the above expressions, it finally results ∀Ω e :

∂a e ∂t + 1 |Ω e | i=n ed i=1 (a e i -a e ) Γ i e γ i e (v • n) dΓ = 1 |Ω e | Ωe
(Aa e + B) dΩ (28)

A technique based on the method of characteristics

Now, we are considering another kind of technique well adapted to compute steady solutions of advection problems in steady recirculating flows as encountered in the driven-cavity flow. We consider, as before, the advection problem for a generic vectorial unknown field (tensor equations can be rewritten in a vector form):

da dt = A(x(t))a + B(x(t)) (29)
In [START_REF] Chinesta | On the steady solution of linear advection problems in steady recirculating flows[END_REF] we have proposed an original technique for computing the steady state of that equation in a steady recirculating flow, that we summarize in the this section. This numerical procedure is based on solving the homogeneous equation and looking for a particular solution of the complete equation. Then, the searched solution results of imposing its periodicity along its associated closed streamline. Thus, the following three steps can be identified.

Computing the general solution of the homogeneous system.

We consider that vector a has M components. We denote by a i h the solution of the homogeneous system

(i = 1, • • • , M ) da i h dt = A(x(t))a i h (30)
related to the initial condition (a i h (t = 0)) j = δ ij , where (a i h ) j denotes the j th component of vector a i h and δ ij the Kroenecker's delta. For computing one of these solutions we make use of the method of characteristics, according to the following first order backward algorithm:

• Initialization: t = 0, x(t = 0) = X, (a i h (t = 0)) j = δ ij ∀i • First step: t ← t + ∆t, x ← X + v(X) ∆t, a i h ← a i h + ∆t(A(X)a i h ) • While x = X -t ← t + ∆t -x ← x + v(x) ∆t -a i h ← a i h + ∆t(A(x)a i h )
Now the general solution of Eq. ( 30) can be written as:

a h (t) = i=M i=1 α i a i h (t) (31)
3.2.2 Computing a particular solution of the complete system.

Now we integrate Eq. ( 29) from an arbitrary initial condition, as for example a p (t = 0) = 0, using the previous algorithm, from which we obtain a particular solution of that equation a p (t) E. Pruliere, A. Ammar, N. El Kissi, F. Chinesta

Imposing the solution periodicity along the closed streamlines.

As the streamlines are closed, and a steady state is assumed, after a time T (that is the period associated with the streamline related to point X that we denote by C(X)) an imaginary fluid particle comes back to the departing point, and now the solution periodicity implies:

a(t = 0) = a h (t = 0) + a p (t = 0) =0 =     α 1 α 2 . . . α M     = a(t = T ) = = i=M i=1 α i a i h (t = T ) + a p (t = T ) (32)
from which the coefficients α i can be obtained and then the steady solution at point X established. The interested reader can refers to [START_REF] Chinesta | On the steady solution of linear advection problems in steady recirculating flows[END_REF] for more details.

Concluding remarks.

Remark 1. The discontinuous Galerkin technique can be applied in both the steady and the transient cases. This discretization technique introduces some amount of numerical diffusion (certainly inevitable) that has a significant impact in the algorithm convergence in transient resolutions. In fact, it is easy to prove that in absence of diffusion the steady solution cannot be reached from an arbitrary initial condition [START_REF] Manhart | Rheology of suspensions of rigid-rod like particles in turbulent channel flow[END_REF]. However, due to the numerical diffusion introduced by the discretization scheme the discontinuous Galerkin scheme leads to a steady solution (the long time solution) even in absence of real diffusion effects.

Remark 2. The integration by the method of characteristics has several advantages: (i) the order of the method can be improved by using high order Runge-Kutta schemes for example; (ii) the time step considered (which induces a characteristic length l with l = ∆t×max x∈C(X) (v(x)) does not depend on the size of the mesh considered for the resolution of the flow kinematics. The integration by characteristics could be also applied for solving transient advection problems, that now, if we start from an arbitrary initial condition no steady long time solution is found, in a fully agreement with the theoretical analysis.

Solving the kinetic theory description of a short fiber suspension

The kinetic theory description of a short fiber suspension is given by the Fokker-Planck equation [START_REF] Azaiez | State-of-the-Art on numerical simulation of fiberreinforced thermoplastic forming processes[END_REF], that for a 2D orientation distribution (here considered for the sake of simplicity) can be written using the angular coordinate, as

dΨ dt + ∂ ∂ϕ Ψ(ϕ) dϕ dt -D r ∂ 2 Ψ ∂ϕ 2 = 0 (33)
or

dΨ dt + E 0 Ψ + E 1 ∂Ψ ∂ϕ + E 2 ∂ 2 Ψ ∂ϕ 2 = 0 ( 34
)
whose variational formulation in the conformation space results

2π 0 dΨ dt Ψ * dϕ + 2π 0 E 0 ΨΨ * dϕ + 2π 0 E 1 ∂Ψ ∂ϕ Ψ * dϕ - 2π 0 E 2 ∂Ψ ∂ϕ ∂Ψ * ∂ϕ dϕ = 0 (35)
subjected to the periodicity conditions Ψ(ϕ = 0) = Ψ(ϕ = 2π).

As the general form of the Fokker-Planck equation involves an advection term in the conformation coordinates we will consider a stabilized SUPG finite element approximation of the advection-diffusion differential operator.

After discretizing the conformation space, the following system of ordinary differential equations is obtained:

K 1 Ψ + K 2 Ψ = 0 (36)
from which it finally results:

Ψ + K -1 1 K 2 Ψ = 0 (37)
Now, for integrating in the physical space (driven-cavity domain) both strategies proposed in the previous section could be applied, as described in the next paragraphs.

Remark 3. If we assume a uniform mesh in the angular coordinate defined by the nodes

ϕ i , i = 1, • • • , N c , then the i-component of Ψ denotes Ψ(ϕ = ϕ i ).

Discontinuous Galerkin approximation.

Using the notation previously introduced, the fully discrete form of Eq. ( 37) results:

∂Ψ e ∂t + 1 |Ω e | i=n ed i=1 (Ψ e i -Ψ e ) Γ i e γ i (v • n) dΓ + 1 |Ω e | Ωe K -1 1 K 2 dΩ Ψ e = 0( 38
)
where the fact that the orientation distribution function is constant into each element has been taken into account.

Integration by characteristics.

In this case the integration along C(X) can be expressed as:

Ψ(t n+1 ) = Ψ(t n ) + ∆t K -1 1 K 2 Ψ(t n ) (39)
or

Ψ(t n+1 ) = I + ∆t K -1 1 K 2 M (n,n+1) Ψ(t n ) ( 40
)
and repeating the same reasoning we can finally write:

Ψ(t n+1 ) = M (0,1) • • • • • M (n,n+1) Ψ(t 0 ) = M (0,n+1) Ψ(t 0 ) (41)
Now, applying the periodicity constraint Ψ(t = T ) = Ψ(t = 0) and assuming T = m ∆t it results:

M 0,n+1 Ψ 0 = Ψ 0 ⇒ M (0,n+1) -I Ψ 0 = 0 (42)
that can be solved with the normality condition 2π 0 Ψ(ϕ)dϕ = 1 for obtaining the steady solution Ψ 0 at point X.

Numerical tests

We consider in this section some simple flows whose associated Fokker-Planck solution can be obtained exactly allowing to conclude about the accuracy of the discretization techniques just proposed.

For this purpose we consider a simple shear flow as well as a simple shear recirculating flow.

Simple shear flow

In this first example the flow kinematics is imposed and it remains uncoupled with the microstructure evolution:

v = Gy 0 (43)
In this particular flow the steady solution must be independent on the physical coordinates, and consequently the Fokker-Planck equation reduces to:

∂ ∂ϕ Ψ(ϕ) dϕ dt -D r ∂ 2 Ψ ∂ϕ 2 = 0 (44)
We can consider two situations according to the value of the diffusion coefficient:

• When diffusion vanishes (D r = 0): ∂ ∂ϕ Ψ(ϕ) dϕ dt = 0 → Ψ(ϕ) dϕ dt = A (45)
where dϕ dt is given by the Jeffery's equation:

φ = dϕ dt = G 2 (-1 + k cos 2ϕ) (46)
from which it directly results:

Ψ(ϕ) = A G/2(-1 + k cos 2ϕ) (47)
where the constant of integration is computed from the normality condition:

2π 0 Ψ(ϕ)dϕ = 1 (48) resulting finally Ψ(ϕ) = √ 1 -k 2 2π(1 -k cos 2ϕ) (49)
• For a non-zero diffusion coefficient the solution of Eq. ( 44) results:

Ψ(ϕ) = e z(ϕ) F 1 ϕ 0 e -z(x) dx + F 0 (50)
where

z(x) = x 0 φ Dr dϕ = G(k sin(2x) -2x) 4Dr (51)
F 0 = e z(π) 1 -e z(π) π 0 e -z(x) dx (52)
and

F 1 = 2π 0 e z(ϕ) ϕ 0 e -z(x) dx + F 0 dϕ (53)
For applying the technique based on the imposition of the periodicity condition in the present problem, that does not involve a recirculating flow, we use the fact that the steady orientation distribution in a simple shear flow does not evolve in space. Thus, one can apply that technique by imposing that the solution at point X must be identical to the one found in any point x located downstream with respect to X (x and X both located on the same streamline). On the other hand, for the application of the discontinuous Galerkin strategy we consider that the solution evolves in time but not in the physical space.

Figure 1 compares the solution obtained using both techniques for a unit shear rate, G = 1, and fibers characterized by the parameter k = 0.8: (i) the one based on the discontinuous Galerkin approximation and (ii) the one related to the integration by characteristics, both with respect to the exact solution. We can notice that both techniques are very accurate and that there are not significant differences between them.

Recirculating shear flow

In this section we consider a bit more complex shear flow which recirculates. The kinematics is now defined by:

v = u v = -y x 2 + y 2 x x 2 + y 2 (54)
being the fibers characterized by the parameter k = 0.8. The velocity field is assumed uncoupled with the orientation distribution in order to derive an exact solution that is then used to conclude about the performance of both discretization techniques. As previously, we consider two cases:

• D r = 0.
It is easy to prove that the steady solution verifies dΨ dt = ∂Ψ ∂ϕ ω, being ω the rotation velocity that for this flow results ω = x 2 + y 2 , being constant on each streamline. Thus, the solution of the resulting Fokker-Planck equation turn out:

∂ ∂ϕ (( φ -ω)Ψ) = 0 ⇒ Ψ(ϕ) = A φ -ω (55)
where the integration constant A is determined again from the normality condition.

• D r = 0. Now, the Fokker-Planck equation results

∂ ∂ϕ (( φ -ω)Ψ) -Dr ∂ 2 Ψ ∂ϕ 2 = 0 (56) which leads to ( φ -ω) Ψ -Dr ∂Ψ ∂ϕ = A
where A is the integration constant.

Eq. ( 56) is very close to Eq. ( 44), and then its solution can be written by considering φ -ω instead φ in the solution previously addressed:

     Ψ(ϕ) = e z(ϕ) F 1 ϕ 0 e -z(x) dx + F 0 F 0 = e z(π) 1-e z(π) π 0 e -z(x) dx F 1 = 2π 0 e z(ϕ) ϕ 0 e -z(x) dx + F 0 dϕ (57)
where

z(x) = x 0 φ -ω Dr dϕ (58)
In order to use the discontinuous Galerkin strategy a mesh was built in the circular domain (see figure 2). As the velocity field is tangent to the domain boundary, there are not fluxes of the orientation distribution across the element edges located on the domain boundary, and in consequence we don't need to prescribe any boundary condition (in advection problems boundary conditions can be only prescribed on the inflow boundary). We have selected three elements (identified in figure 2) where the solutions computed by using the different numerical strategies are compared.

The comparison of both techniques in the limit when D r = 0 is depicted in figure 2, where we can notice the extremely superior accuracy obtained when the integration by characteristics is performed. Moreover, as noticed in some of our former works (see [START_REF] Chaidron | On the periodicity of the extra-stress tensor in steady recirculating viscoelastic flows[END_REF] for instance) the quality of the stabilized finite element solution is degraded as we approach the center of rotation whose singularity diffuses towards its neighborhood along the normal direction to the streamlines. To overcome this problem very fine meshes must be considered in this region to capture the high gradients of the velocity (in particular its direction) and the orientation distribution fields. On the contrary the integration by characteristics is not affected but the presence of the singularity, but when it will be coupled with a finite element kinematics solver, finer meshes will be required in the center of rotation neighborhood to guarantee the solution accuracy.

The driven cavity flow

In this section the driven cavity flow problem is simulated by using the most accurate integration procedure of both just analyzed. It consists of an integration by characteristics (in the physical space) coupled with a standard SUPG finite element discretization in the conformation space. However, despite the right scheme accuracy, the proposed algorithm requires significant computing efforts. By this reason, we propose in this section an efficient model reduction, that has been successfully applied in some of our former works for simulating complex fluid flow problems [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF].

The cavity domain consists of a unit length square. The flow kinematics is computed in the mesh depicted in figure 4 assuming a Newtonian solvent fluid with unit viscosity. The fluid velocity is prescribed on the cavity boundary, vanishing on the lower and lateral walls, and taking a parabolic profile on the upper wall (avoiding velocity discontinuities at the upper corners):

v(x, y = 1) = 16x 2 (1 -x) 2 V max 0 ( 59
)
where V max was set to V max = 10D r . Figure 4 depicts the computed velocity field as well as the reconstructed stream function Ξ(x). The orientation field is depicted in figure 5. In that figure the ellipsoids axes are related to the second order orientation tensor eigenvalues, that represent the probability of finding the fiber aligned on the direction of the associated eigenvector.

Advanced strategies: techniques based on model reduction

In the previous section the driven cavity flow problem was simulated by using an integration by characteristics (in the physical space) coupled with a standard SUPG finite element discretization in the conformation space. However, despite the right scheme accuracy, the proposed algorithm requires significant computing efforts. By this reason, we propose in this section some new advanced model reduction strategies, that were successfully applied in some of our former works for simulating complex fluid flow problems.

The first strategy consists in computing the solution along one or several flow streamlines, and then extract the most significant characteristics of the solution in the conformation space by applying the proper orthogonal decomposition. This technique was successfully applied in [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF] and throughout the present section will be summarized and applied for solving the computing the fluid microstructure in the flow cavity problem.

The second strategy that we are analyzing lies in the application of a separated representation of the trial and test fields involved in the variational formulation related to the suspension microstructural state. Separated representations and tensor product approximation basis were successfully used in some of our former works (see [START_REF] Ammar | A new family of solvers for some clases of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] [6] and the references therein) for solving different kinetic theory models arising in complex fluids modelling.

In the next paragraphs we are summarizing the basic ideas related to the proper orthogonal decomposition that will be then applied for solving the kinetic theory model related to the microstructural state of the fiber suspension in a driven cavity flow.

The proper orthogonal decomposition revisited

We assume that the evolution of a certain field T (x, t) is known. For the sake of simplicity we assume in what follows that this field is scalar, but all the developments can be extended to general tensorial fields. In practical applications, this field is defined at the nodes of a spatial mesh x i and for some times

t m = m × ∆t, with i ∈ [1, • • • , N ] and m ∈ [0, • • • P ]. We introduce the notation T (x i , t m ) = T m (x i ) ≡ T m i (t m
). T m defines the vector containing the nodal degrees of freedom at time t m . The main idea of the proper orthogonal decomposition is to obtain the most typical or characteristic structure φ(x) among these T m (x), ∀m. The maximization of:

α = m=P m=1 i=N i=1 φ(x i )T m (x i ) 2 i=N i=1 (φ(x i )) 2 (60) leads to: m=P m=1 i=N i=1 φ(x i )T m (x i ) j=N j=1 φ(x j )T m (x j ) = α i=N i=1 φ(x i )φ(x i ); ∀ φ (61)
which can be rewritten in the form

i=N i=1 j=N j=1 m=P m=1 T m (x i )T m (x j )φ(x j ) φ(x i ) = α i=N i=1 φ(x i )φ(x i ); ∀ φ (62)
Defining the vector φ such that its i-component is φ(x i ), Eq. ( 62) results in the eigenvalue problem (63), whose eigenvectors related to the highest eigenvalues define the characteristic structure of T m (x),

φT c φ = α φT φ; ∀ φ ⇒ c φ = αφ (63)
where the two points correlation matrix is given by

c ij = m=M m=1 T m (x i )T m (x j ) ⇔ c = m=M m=1 T m (T m ) T (64)
which is symmetric and positive definite. If we define the matrix Q containing the discrete field history:

Q =      T 1 1 T 2 1 • • • T P 1 T 1 2 T 2 2 • • • T P 2 . . . . . . . . . . . . T 1 N T 2 N • • • T P N      (65)
then it is easy to verify that the matrix c in Eq. ( 63) results

c = Q Q T (66) 6.1.1 A posteriori reduced modelling.
To define a reduced model we start solving the eigenvalue problem defined by Eq. ( 63) and then selecting the eigenfunctions φ k associated with the eigenvalues belonging to the interval defined by the highest eigenvalue and that value divided by a large enough value (10 8 in our simulations). In practice n is much lower than N . Now, we could try to use these n eigenfunctions φ k for approximating the solution T m (x), ∀m. For this purpose we need to define the matrix

B = [φ 1 • • • φ n ] B =     φ 1 (x 1 ) φ 2 (x 1 ) • • • φ n (x 1 ) φ 1 (x 2 ) φ 2 (x 2 ) • • • φ n (x 2 ) . . . . . . . . . . . . φ 1 (x N ) φ 2 (x N ) • • • φ n (x N )     (67)
Now, if we consider the linear system of equations resulting from the semi-implicit discretization of a generic parabolic partial differential equation:

G m T m+1 = H m (68)
expressing

T m+1 = i=n i=1 ζ m+1 i φ i = B ζ m+1 (69)
equation ( 68) results

G m T m+1 = H m ⇒ G m B ζ m+1 = H m (70)
and by multiplying both terms by B T we obtain

B T G m B ζ m+1 = B T H m (71)
where the size of

B T G m B is n × n, instead N × N .
When n ≪ N , as is the case in numerous physical models, the solution of Eq. ( 71) is preferred because its reduced size.

Remark 4. Equation ( 71) can be also deduced by introducing the approximation (69) into the Galerkin variational form related to the partial differential equation. If trial and test functions in the vartiational formulation of a parabolic partial differential equation are approximated using the form (69) then Eq. ( 71) is easily deduced.

Remark 5. An alternative technique to reduce the size of the eigenvalue problem lies in the application of the snapshot proper orthogonal decomposition. This technique is based on finding the significant modes from the application of the POD, but assuming that those modes can be written as a linear combination of the M snapshots that were used to define the decomposition. The main advantage of this strategy is that theses modes result from the eigenproblem defined by

Q T Q Υ = λ Υ ( 72
)
whose size is P × P instead of N × N . Only the eigenvectors related to large enough eigenvalues are retained. From these eigenvectors the reduced approximation functions are computed using the fact that these functions are linear combination of the snapshots, i.e.:

φ i = Q Υ i (73)
Remark 6. In the previous analysis the reduced basis was built from the computed unknown field evolution that was carried out solving the discrete evolution problem. Thus, one could ask about the interest of a such approach. There are two kind of approaches widely considered. The first approach consists in solving the non-reduced model in a short time interval, allowing the extraction of the characteristic functions and then the definition of the reduced approximation basis that is then used to perform the solution of the reduced evolution model in larger time intervals with the associated computing time savings. The other approach consists in solving the non-reduced model in the whole time interval, whose solution allows defining the reduced approximation basis that then could be used for solving "similar" models, as the ones that involve for example slight variations in some material parameters or in the boundary conditions. Some recent advances in such approaches can be found in [START_REF] Park | The use of the Karhunen-Love decomposition for the modelling of distributed parameter systems[END_REF], [START_REF] Maday | The reduced basis element method: application to a thermal fin problem[END_REF], [START_REF] Ryckelynck | A priori hyper-reduction method: an adaptive approach[END_REF], [START_REF] Ryckelynck | An efficient "A Priori" model reduction for boundary element models[END_REF], [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF] [9], [START_REF] Burkardt | POD and CVT-based reduced-order modeling of Navier-Stokes flows[END_REF], [START_REF] Gunzburger | Reduced-order modeling of time-dependent PDEs with multiple parameters in the boundary data[END_REF] and the references therein.

Towards and adaptive strategy: Enriching the approximation basis.

Obviously, accurate simulations require an error evaluation as well as the possibility of adapting the approximation basis by introducing new functions able to describe the solution features. Ryckelynck proposed in [START_REF] Ryckelynck | A priori hyper-reduction method: an adaptive approach[END_REF] an adaptive procedure, able to construct or enrich the reduced approximation basis. For this purpose, he proposed to enrich the reduced approximation basis by adding some Krylov subspaces generated by the equation residual.

Despite the fact that this enrichment tends to increases the number of approximation function, when it is combined with a POD decomposition that continuously reduces this number, the size of the problems is quickly stabilized. This strategy, which was successfully used in [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] as well as in [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF], is summarized in the present section.

We assume the solution accurately described in the interval ]0, t s ] by using the reduced basis B. Now, the solution evolution is performed in the time interval ]t s , t s+1 ] solving Eq. (71) making use of the reduced approximation basis defined by matrix B:

ζ m+1 = B T G m B -1 B T H m (74)
When time t s+1 is reached, a control step is performed in order to evaluate the accuracy of the solution computed using the reduced basis. The control step is performed only at the end of each time interval.

We assume that t s+1 -t s = M × ∆t and consequently the residual at t s+1 , R M , can be computed from

R M = G M -1 B ζ M -H M -1 (75)
If the norm of the residual is small enough R M < ǫ B ζ M (being ǫ a small enough parameter) the computed solution can be assumed as good, and the time integration goes on in ]t s+1 , t s+2 ] using Eq. (74) without changing the reduced approximation basis.

On the contrary, if R M ≥ ǫ B ζ M , the approximation must be improved. For this purpose, we propose to enrich the reduced approximation basis by introducing the just computed residual (or some Krylov's subspaces generated by the residual, as proposed in [START_REF] Ryckelynck | A priori hyper-reduction method: an adaptive approach[END_REF]):

B ← [B R M ] ( 76
)
and now, the evolution is recomputed in the ]t s , t s+1 ] using Eq. ( 74) with the just updated reduced basis B. When the convergence criterion is satisfied (i.e. R M < ǫ B ζ M ) we apply a POD on the entire past history ]0, t s+1 ] for defining the reduced basis able to represent all the past evolution of the unknown field, which results in an updated B that can be considered as the optimal reduced basis to describe the evolution of the unknown field in ]0, t s+1 ]. Using the just updated reduced basis B the time evolution given by Eq. ( 74) is performed in ]t s+1 , t s+2 ] and the solution accuracy is checked at t = t s+2 . If a basis enrichment was needed at the end of the previous time interval, at t = t s+1 , then the length of the present time interval is reduced according to t s+2 -t s+1 = (t s+1 -t s )/2, and if not enrichment was needed, then the time interval length is increased according to t s+2 -t s+1 = 2 (t s+1 -t s ) (see [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] for more details on the time interval length adaptation).

The enrichment tends to increase the size of the reduced basis, but the POD reduces its size. The combination of both procedures allows to stabilize the size of the model as we noticed in some of our former works [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF].

Building-up the reduced kinetic theory model of the fiber suspension

We consider an arbitrary point in the flow domain, at which an isotropic fiber distribution is assumed. Now, the evolution of the orientation distribution is computed along the associated streamline using the "a posteriori" model reduction algorithm just proposed applied on Eq. [START_REF] Muller | A purely elastic transition in Taylor-Couette flow[END_REF], and the solution at different times t i are stoked in matrix Σ, whose i-th column contains the orientation distribution at each node used to discretize the conformation space at time t i . Now, we can extract the most significant information related to the orientation distribution evolution by applying a proper orthogonal decomposition. For this purpose we solve the eigenproblem defined by

(ΣΣ T ) V = λ V (77)
which leads to N c couples V j , λ j (N c being the number of nodes used in the discretization of the conformation space). Now, only the must significant eigenvectors are retained V 1 • • • V r , with r << N c according to the criterion:

λ i max λ i > 10 -4 (78)
Then, a good approximation basis can be built from these significant eigenvectors:

Ψ = Ψ 0 + n i=1 V i ξ i (79)
which can be written in a matrix compact form as:

Ψ = Bξ (80)
that is inserted into the discrete form of the Fokker-Planck equation

B T K 1 B K 1 ξ + B T K 2 B K 2 ξ = 0 (81)
that defines a linear system of size r instead N c .

Results

To conclude about the accuracy of the just proposed reduced modelling we depict in figure 6 the error associated with the second order orientation tensor, defined in Eq. ( 5), which is computed at each element of the flow domain mesh by:

E e = Tr((a e -a e red ) 2 ) ( 82
)
where a e is the solution computed at element Ω e using the global approximation basis, and a e red the one computed from the reduced model than only contains around 10 approximation functions.

We can also define the mean error from:

E = Ne e=1 E e N e (83)
that results for the error depicted in figure 6: E = 1.3 × 10 -5 . We can notice the excellent accuracy obtained using the reduced modelling, despite the significant reduction of the CPU time. The problem solution using model reduction is 600 times faster.

The solution can be improved by considering different flow streamlines (or the same with different initial conditions), which lead to different reduced bases: V (j) i , i ≤ r (j) , where (j) refers to the different computed reduced bases. Now, from these bases we can extract a global one

V 1 • • • V R by adding to the first one V 1 • • • V (1) r
only the vectors contained in the other reduced bases whose projection on the previous updated basis is lower than a certain tolerance (to maximize their linear independence).

Thus, when in the flow just analyzed another trajectory closer to the flow domain walls is considered, the solution is improved as we can notice in figure 7, being the mean error in this case E = 8 × 10 -8 (using 15 approximation functions in the reduced basis).

Empirical natural closure relation

As was indicated in the first section, the description of the microstructure state by means of the moments of the orientation distribution function (orientation tensors) is suitable from a computational point of view because the important CPU times savings that its use implies, but this approach requires the introduction of a closure relation. Today, it is widely accepted that the use of such a closure relation can have a significant, and sometimes unpredictable, impact on the computed solution. Numerous closure relations have been proposed, but there is not a universal or exact relation to be used in any flow conditions.

In this section we are defining an empirical and local natural closure relation that can be adapted locally to each flow conditions. This procedure was deeply described in [START_REF] Pruliere | Empirical natural closure relation for short fiber suspension models[END_REF] and its main ideas are revisited in the present section.

Fundamentals

Any closure approximation of the fourth order orientation tensor can be written using the Caley-Hamilton theorem as well as partial normalization and symmetry [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF] in the general form:

   A = i=6 i=1 B i β i in 3D A = i=3 i=1 B i β i in 2D (84)
where the different fourth order tensors involved in the former expression write:

             B 1 = S(I ⊗ I) B 2 = S(I ⊗ a) B 3 = S(a ⊗ a) B 4 = S(I ⊗ a 2 ) B 5 = S(a ⊗ a 2 ) B 6 = S(a 2 ⊗ a 2 ) ( 85
)
where S( ) refers to the symmetric component. Now, taking into account the normalization conditions, we obtain the relations that the coefficients β i verify in the three-dimensional case:   10 1 0 0 7 2 0 0 4

    β 1 β 2 β 3   =   1 -2D 4P 4P 0 1 -6D 4(P -D) 7 5 2(3 -4D)     β 4 β 5 β 6   +   0 6 0   (86)
and in two-dimensional:

8 1 0 1 β 1 β 2 = 4P -1 β 3 + 0 1 ( 87
)
where P = det(a) and D = 1 2 1 -Tr(a 2 ) . In the general case, we can write Eqs. ( 86) and (87) in the matrix form:

f β u = gβ d + h (88)
Thus, the matrix form of Eq. ( 84) becomes

A = Wβ = W β u β d = W f -1 g I β d + W f -1 h 0 = Wβ d + W (89)
where now A denotes de vector form of the fourth order orientation tensor. Now, considering the former equation for the n couples of tensor (A i , a i ) extracted via the Karhunen-Loève decomposition from the most significant distribution functions (obtained by solving the Fokker-Planck equation along some streamlines), one could be interested in determining the optimal value of β d , as originally considered in [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF] to derive the natural closure relation. However, in [START_REF] Dupret | Modelling the flow of fiber suspensions in narrows gaps[END_REF], authors considered only the orientation distribution related to an initial isotropic distribution of infinite slender fibers, without diffusion effects, immersed in a simple shear flow.

Computing the local natural closure relation

For computing the unknown beta coefficients we proceed by using a moving least square technique. For this purpose we define the approximation error

J(a, β d ) = n i=1 δ ǫ ( a -a i ) γ i A i -Wi β d -Wi 2 (90)
where the weighting coefficient γ i scales with the value of the eigenvalue related to the eigenfunction that served for computing tensors A i and a i . The window function δ ǫ ( ) is introduced for increasing the weight of closer couples in the error expression. In our simulations we have considered the following window function:

δ ǫ (d) = δ d ǫ ǫ with δ d ǫ = e -( d ǫ ) 2 √ π (91)
In the case of ǫ → ∞ we can write δ ǫ = 1 √ π . Thus, J becomes independent of a and the expression (90) reduces to the usual global least squares. In this case the vector β d does not require to be updated, which results in significant computing time savings.

In the general case, the minimization of the functional (90) leads to:

∂J ∂β d = 0 = 2 n i=1 δ ǫ ( a -a i ) γ i WT i A i -Wi β d -Wi (92)
which results in:

n i=1 δ ǫ ( a -a i ) γ i WT i Wi β d = n i=1 δ ǫ ( a -a i ) γ i WT i A i -Wi (93)
or

β d = n i=1 δ ǫ ( a -a i ) γ i WT i Wi -1 n i=1 δ ǫ ( a -a i ) γ i WT i A i -Wi (94)
that allows computing for each value of a the coefficients of the empirical closure relation β d , from which one could compute β u from Eq. (88). Both, β d and β u , define completely the closure relation:

A = W(a)β(a) (95)
The use of the just described procedure only requires the extraction of the couples (A i , a i )

Extracting the characteristic structure of the orientation distribution

For extracting the characteristic structure of the orientation distribution one could solve the Fokker-Planck equation (that does not involve any closure relation) by using one of the techniques previously described, along some flow streamlines.

Thus, one could consider different snapshots of the orientation distribution (stored in the matrix Σ), from which the characteristic structure of the solution evolution could be determined by direct application of the Karhunen-Loève decomposition. Now, we could define the couples (A i , a i ) using three different strategies:

1. The simplest alternative lies in applying the Karhunen-Loève decomposition to the solution history stored in Σ. Obviously, the resulting eigenfunctions do not have a fully physical meaning, and therefore they are not strictly positives. Thus, if we compute the second and fourth order moments (A i , a i ) associated with the model eigenfunctions, then these tensors have unusual properties (e.g. a i do not have a unit trace, ...). However, one could forget these "conceptual difficulties" and apply the procedure described in the previous section to the couples (A i , a i ).

2. Other possibility to avoid these modes without a full physical meaning consists of assuming that the orientation distribution can be written everywhere on the considered streamline as a linear combination of some orientation distributions (snapshots), for example, those were used for applying the Karhunen-Loève decomposition. In order to avoid the consideration of two distributions too close, one could start with the initial distribution, adding a new distribution to the reduced approximation basis only when it is far enough from all the functions defining the reduced approximation basis. The advantage of this approach is that all the modes used for approximating the orientation distribution evolutions are real physical distributions and consequently the associated second order orientation tensors are perfectly defined from all points of view.

3. Finally, one could consider a third possibility related to the snapshot Karhunen-Loève decomposition. This technique is based on finding the significant modes from the application of the Karhunen-Loève decomposition, but assuming that those modes can be written as a linear combination of the N snap snapshots that were used to define the decomposition. The main advantage of this strategy derives from the fact that it is easy to prove that theses modes result from the eigenproblem defined by

Σ T Σ Υ = λΥ ( 96
)
whose size is N snap ×N snap instead of N c ×N c (N c being the number of degrees of freedom used to define the approximation in the configuration space -the one related to the orientation -), and where only the eigenvectors Υ i related to large enough eigenvalues λi are retained. From those eigenvectors the reduced approximation modes are computed using the fact that those modes are linear combination of the snapshots, i.e.:

V i = ΣΥ i (97)

Numerical tests

In this section we are addressing two rheometric simple flows as well as the complex driven cavity flow problem.

• Simple shear flow.

Firstly, we consider the 2D flow defined by the velocity field v T = ( γy, 0) with γ = 1. The fiber aspect ratio was adjusted to have k = 0.8, being the diffusion coefficient D r = 0.1.

The 2D orientation distribution is computed using a stabilized SUPG finite elements discretization in the conformation space, on the mesh associated with the conformation space that consists of 200 nodes uniformly distributed in [0, 2π[.

When the Fokker-Plank equation is integrated along a trajectory from the isotropic initial state until reaching the steady state characterized by a small enough variation between two consecutive time steps, that is:

Ψ(t n+1 ) -Ψ(t n ) ≤ 10 -6 Ψ(t 1 ) -Ψ(t 0 ) ( 98
)
The three procedures described in the previous sections are applied to compute the characteristic structure of the solution (orientation distribution) evolution, that will be then used to compute the couples (A i , a i ). The functions defining the reduced approximation basis obtained by applying the three procedures are depicted in Figs. 8, 9 and 10 respectively. We can notice that only 5 functions are needed for approximating accurately the evolution of Ψ.

In order to compare these different strategies, we have computed the reference second order orientation tensor from the resolution of the Fokker-Planck equation, and then the different orientation tensors using the closures proposed in the previous section. The reference and approximated second order orientation tensors were compared at different times, for different values of the diffusion coefficients, using the following error measure:

E = a ref -a appr (99)
The table that follows groups the different computed errors at the end of the evolution process. We can notice that when the diffusion coefficient is small enough the standard natural closure relation works perfectly because we are close to the hypothesis under which it was derived. However, when the diffusion coefficient is increased, the solution accuracy is degraded. On the other hand, the adaptive natural closure that we have defined works perfectly for any diffusion coefficient, except for the too small values that induce numerical difficulties in the finite element discretization of the Fokker-Planck equation in the conformation space. We must mention that due to the fact that the Karhunen-Loève and the Snapshot Karhunen-Loève natural closures are built from modes (approximation functions) that have not a full physical meaning we have preferred using a global least square technique. Moreover, the global least squares fitting can be applied without a significant impact on the solution accuracy and with a considerable computing time reduction. There is no significative differences between the global and local (moving) least squares. Because the global least squares is suitable from the computational point of view (as was indicated above) in that follows we will consider ǫ = ∞.

• Extensional flow.

We consider a 2D extensional flow whose kinematics is defined by

gradv = ǫ 0 0 -ǫ (100)
being the elongation rate ǫ = 1. The fiber suspension is characterized by fibers whose shape factor leads to k = 0.8. The orientation distribution (assumed 2D) is discretized using the SUPG finite element method on the mesh associated with the conformation space that consists of 200 nodes uniformly distributed in [0, 2π[. Figure 11 depicts the time evolution of the error defined in Eq. ( 99) for D r = 0.01 and D r = 0.1 where only the snapshot natural closure is compared with the standard natural one. We can notice the superior accuracy of the empirical snapshot natural closure for both diffusion coefficients.

The table that follows groups the errors, computed at the end of the time simulation, computed by using the different natural closure relations. Again, the best results are obtained by using the empirical snapshot natural closure. • The driven cavity flow problem. Finally, we consider the complex flow generated in a driven cavity that involves a short fiber suspension characterized by the same material parameters that in the previous tests. The velocity is prescribed on the domain boundary according to: v(x = 0, y) = 0, v(x = 1, y) = 0, v(x, y = 0) = 0 and v(x, y = 1) = (16v max x 2 (1 -x) 2 , 0). The velocity field is then solved by assuming a Newtonian behavior and by applying a standard mixed finite element formulation where the velocity and pressure approximations verify the LBB stability condition. The Fokker-Planck equation governing the evolution of the fiber orientation distribution function is then solved along some closed streamlines, where the periodicity condition of that distribution function was imposed as described in Ammar and Chinesta [START_REF] Ammar | A particle strategy for solving the Fokker-Planck equation governing the fiber orientation distribution in steady recirculating flows involving short fiber suspensions[END_REF]. From the computed orientation distribution function, the characteristic modes are extracted by using the technique previously described based on the application of the Karhunen-Loève decomposition, allowing to fit the empirical snapshot natural closure strategy previously introduced. Now, the evolution equation associated with the second order orientation tensor is solved by assuming different closure relations: linear, quadratic, hybrid, natural and the empirical natural one.

If the Fokker-Planck equation is solved in the whole domain considering a fine enough mesh, the fiber distribution function can be computed everywhere without addressing any closure relation. Thus, the second order orientation tensor related to the just computed distribution function can be used as reference solution for comparison purposes.

The error between the computed second order orientation tensor (computed by using the most usual closure relations) and the reference solution coming from the solution of the Fokker-Planck equation is reported in the following In general, models arising from the kinetic theory descriptions here defined in moderate or even highly multidimensional spaces that include the physical space, time and the conformation coordinates whose number depends on the level of richness of the physical model. Thus, the dimension of models ranges from 6D for the simplest kinetic theory models describing for example the fiber orientation to thousands or millions of dimensions for those arising from the kinetic theory of systems involving macromolecules or the ones coming from quantum chemistry in which the dimension of the space in which the wavefunction is defined scales linearly with the number of particles involved in the system. The solution of models like the one just described needs for new strategies, because the standard ones exhibit curse of dimensionality. Some strategies have been recently proposed for solving models defined in multi-dimensional spaces. The sparse grid techniques [START_REF] Bungartz | Sparse grids[END_REF] are one of the most popular, but separated representations were also used for solving the models encountered in quantum mechanics (Hartree-Fock and post-Hartree-Fock techniques).

The sparse grid are restricted (at least in our knowledge) for treating models involving up to twenty dimensions [START_REF] Achdou | Computational methods for option pricing[END_REF]. On the other hand, separated representations like the one considered in the multi-configuration-self-consistent-fields (MCSCF), fix a priori the number of products (anti-symmetrized) and also assumes that those products are constructed by combining a certain number of one-dimensional functions.

Recently we proposed a new strategy able to solve highly multi-dimensional models circumventing the curse of dimensionality. This technique was successfully applied for treating some multi-dimensional models encountered in the kinetic theory description of complex fluids [START_REF] Ammar | A new family of solvers for some clases of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF] [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representation[END_REF]. It allows defining the optimal number of products containing the optimal one-dimensional functions. In what follows we are summarizing the main ideas that the aforementioned technique involves. For the sake of simplicity, we are illustrating the solution procedure by solving the Poisson problem defined in a space of dimension D.

Illustrating the separated representation of the 2D Poisson's problem

△T = -f (x 1 , x 2 , ..., x D ), ( 101
)
where T is a scalar function of (x 1 , x 2 , ..., x D ). Problem (101) is defined in the domain (x 1 , x 2 , ..., x D ) ∈ Ω =] -L, +L[ D with vanishing boundary conditions. The problem solution can be written in the form:

T (x 1 , x 2 , ..., x D ) = ∞ j=1 α j D k=1 F kj (x k ), (102)
where F kj is the j th basis function, with unit norm, which only depends on the k th coordinate.

It is well known that the solution of numerous problems can be accurately approximated using a finite (sometimes very reduced) number of approximation functions, i.e.:

T (x 1 , x 2 , ..., x D ) = Q j=1 α j D k=1 F kj (x k ). ( 103
)
The previous expression implies the same number of approximation functions in each dimension, but each one of these functions could be expressed in a discrete form using different number of parameters (nodes of the 1D grids). Now, an appropriate numerical procedure is needed for computing the coefficients α j as well as the Q approximations functions in each dimension.

The proposed numerical scheme consists of an iteration procedure that solves at each iteration n the following three steps:

Step 1: Projection of the solution in a discrete basis If we assume the functions F kj (∀j ∈ [1, ..., n]; ∀k ∈ [1, ..., D]) known (verifying the boundary conditions), the coefficients α j can be computed by introducing the approximation of T into the Galerkin variational formulation associated with Eq. (101):

Ω ∇T * • ∇T dΩ = Ω T * f dΩ. (104)
Introducing the approximation of T and T * :

T (x 1 , x 2 , ..., x D ) = n j=1 α j D k=1 F kj (x k ), (105)
and

T * (x 1 , x 2 , ..., x D ) = n j=1 α * j D k=1 F kj (x k ), (106)
we have

Ω ∇   n j=1 α * j D k=1 F kj (x k )   • ∇   n j=1 α j D k=1 F kj (x k )   dΩ = = Ω   n j=1 α * j D k=1 F kj (x k )   f dΩ. (107)
Now, we assume that f (x 1 , • • • , x D ) can be written in the form

f (x 1 , • • • , x D ) = m h=1 D k=1 f kh (x k ). (108)
Eq. ( 107) involves integrals of a product of D functions each one defined in a different dimension. Let D k=1 g k (x k ) be one of these functions to be integrated. The integral over Ω can be performed by integrating each function in its definition interval and then multiplying the D computed integrals according to:

Ω D k=1 g k (x k ) dΩ = D k=1 L -L g k (x k )dx k , (109)
which makes possible the numerical integration in highly dimensional spaces. Now, due to the arbitrariness of the coefficients α * j , Eq. (107) allows to compute the n-approximation coefficients α j , solving the resulting linear system of size n × n. This problem is linear and moreover rarely exceeds the order of tens of degrees of freedom. Thus, even if the resulting coefficient matrix is densely populated, the time required for its solution is negligible with respect to the one required for performing the approximation basis enrichment (step 3).

Step 2: Checking convergence From the solution of T at iteration n given by Eq. ( 105) we compute the residual Re related to Eq. ( 101):

Re = Ω (△T + f (x 1 , • • • , x D )) 2 T . (110)
If Re < ǫ (epsilon is a small enough parameter) the iteration process stops, yielding the solution T (x 1 , • • • , x D ) given by Eq. (105). Otherwise, the iteration procedure continues.

The integral in Eq. ( 110) can be written as the product of one-dimensional integrals by performing a separated representation of the square of the residual. E. Pruliere, A. Ammar, N. El Kissi, F. Chinesta

Step 3: Enrichment of the approximation basis From the coefficients α j just computed the approximation basis can be enriched by adding the new function N k=1 F k(n+1) (x k ). For this purpose we solve the non-linear Galerkin variational formulation related to Eq. ( 101):

Ω ∇T * • ∇T dΩ = Ω T * f dΩ, (111)
using the approximation of T given by:

T (x 1 , x 2 , ..., x D ) = n j=1 α j D k=1 F kj (x k ) + D k=1 R k (x k ). ( 112
)
The solution of Eq. ( 111) can be expressed from the stationarity of functional J(T ):

J(T ) = Ω 1 2 (∇T ) 2 dΩ - Ω T f dΩ, (113)
which results in:

δJ = Ω ∇δT • ∇T dΩ - Ω δT f dΩ = 0, (114)
that corresponds to Eq. ( 111) by putting δT = T * . Now, taking the variation of T according to its expression given by Eq. ( 112), where the variation of the known functions F kj vanishes, we have:

δT = δ(R 1 (x 1 ) • • • R D (x D )) = = δR 1 (x 1 )R 2 (x 2 ) • • • R D (x D ) + • • • + R 1 (x 1 )R 2 (x 2 ) • • • δR D (x D ), (115)
that can be written as

T * (x 1 , x 2 , ..., x D ) = = R * 1 (x 1 )R 2 (x 2 ) • • • R D (x D ) + • • • + R 1 (x 1 )R 2 (x 2 ) • • • R * D (x D ). (116)
This leads to a non-linear variational problem, whose solution allows to compute the N functions R k (x k ). Functions F k(n+1) (x k ) are finally obtained by normalizing, after convergence of the non-linear problem, the functions R 1 , R 2 , ..., R D .

To solve this problem we introduce a discretization of those functions R k (x k ). Each one of these functions is approximated using piecewise 1D linear functions. If we assume than p k nodes are used to construct the interpolation of function R k (x k ) in the interval [-L, L], then the size of the resulting discrete non-linear problem is k=D k=1 p k . The price to pay for avoiding a whole mesh in the multidimensional domain is the solution of a non-linear problem. However, even in high dimensions the size of the non-linear problems remains moderate and no particular difficulties have been found in its solution.

Concerning the computation time, even when the non-linear solver converges quickly, this step consumes the main part of the global computing time. Different non-linear solvers have been analyzed: fixed-point, Newton or one based on an alternating directions fixedpoint scheme.

First convergence analysis reveals that even when one uses piecewise linear 1D finite elements for approximation all the one dimensional functions involved by the separated representation, the convergence rate (in the energy norms) is fourth order. This results is a little bit strange because we expected the usual second order convergence. We are performing further analysis to check if there is some super-convergence phenomena. At the same time one could modify the one-dimensional approximations to enhance the approximation consistency and then the convergence rates.

Matrix form of the 2D Poisson's problem separated representation

In what follows we give the matrix form of both steps of the numerical procedure. From a practical point of view all the functions must be defined in a discrete form. Thus, functions F j (x) (respectively G j (y)) and R(x) (respectively S(y)) are defined using a 1D finite element interpolation that in our simulations is assumed piecewise linear. We denote by N (resp. M) the vector containing the value of the p (resp. q) shape functions N p (x) (resp. M q (y)). Finally F j , G j , R and S are the nodal description of those functions.

• Computing of the alpha coefficients: Projection stage.

At iteration n we can write T in the following matrix form :

T (x, y) = N T F 1 M T G 1 N T F 2 M T G 2 . . . N T F n M T G n     α 1 α 2 . . . α n     (117)
from which the gradient reads:

∂T ∂x ∂T ∂y = B α = = dN T F 1 M T G 1 dN T F 2 M T G 2 . . . dN T F n M T G n N T F 1 dM T G 1 N T F 2 dM T G 2 . . . N T F n dM T G n     α 1 α 2 . . . α n     , (118)
where dN (resp. dM) represents the vector containing the value of the shape functions derivatives with respect to the coordinate x (resp. y) at point x (resp. y).

The weighting function related to the Galerkin's variational formulation T * can be expressed in the same form

T * (x, y) = α * 1 α * 2 . . . α * n      F T 1 NG T 1 M F T 2 NG T 2 M . . . F T n NG T n M      , (119)
and its gradient reads

∂T * ∂x ∂T * ∂y = α * 1 α * 2 . . . α * n      F T 1 dNG T 1 M F T 1 NG T 1 dM F T 2 dNG T 2 M F T 2 NG T 2 dM . . . . . . F T n dNG T n M F T n NG T n dM      = = α * T B T . (120)
Now, we consider the variational formulation related to Eq. ( 101)

Ω GradT GradT * dΩ = Ω f (x, y) T * dΩ, (121)
where the fact that T * vanishes at the domain boundary has been introduced (both the trial and the test functions are in the Sobolev's space H 1 0 (Ω)). Moreover, we assume that f (x, y) can be also written in the form:

f (x, y) = m h=1 a h (x)b h (y). ( 122
)
This particular form will be justified later. Thus, Eq. ( 121) yields:

Ω [ ∂T * ∂x ∂T * ∂y ] ∂T ∂x ∂T ∂y dΩ = m h=1 Ω a h (x)b h (y)T * (x, y)dΩ. ( 123
)
Introducing the discrete form of the trial and test functions as well as their derivatives, we obtain:

Ω α * T B T B α dΩ = = m h=1 Ω α * 1 α * 2 ... α * n      a h (x)F T 1 N b h (y)G T 1 M a h (x)F T 2 N b h (y)G T 2 M . . . a h (x)F T n N b h (y)G T n M      dΩ, (124)
where the components of matrix B T B are given by:

(B T B) ij = F T i dNdN T F j G T i MM T G j + + F T i NN T F j G T i dMdM T G j . ( 125
)
Because all the terms in Eq. ( 124) are the product of functions of x and functions of y, the integral in the domain Ω can be separated as the product of integrals in x and y. Thus, we have:

α * T K α = α * T V ⇒ K α = V, (126)
where

K ij = F T i x=+L x=-L dNdN T dx F j G T i y=+L y=-L MM T dy G j + +F T i x=+L x=-L NN T dx F j G T i y=+L y=-L dMdM T dy G j , (127)
and

V i = m h=1 F T i x=+L x=-L N a h (x)dx G T i y=+L y=-L M b h (y)dy . (128)
The linear system (126) can be solved for computing the coefficients α. The inversion of matrix K is very simple due its reduced dimension. Of course the iteration procedure stops when α n+1 /max(α 1 , • • • , α n ) < ǫ (ǫ being a small enough parameter).

• Enrichment of the approximation basis:

At this stage, we are looking for a solution such that:

T (x, y) = j=n j=1 α j F j (x)G j (y) + R(x)S(y) T * (x, y) = R * (x)S(y) + R(x)S * (y) . (129)
In order to compute both R and S from the associated vartiational formulation (Eq. ( 121)), we start by writing both functions as well their gradients using a matrix form:

∂T ∂x ∂T ∂y = n j=1 α j dN T F j M T G j N T F j dM T G j + M T SdN T 0 (1,q) 0 (1,p) N T RdM T R S = = n j=1 α j D j + E R S , (130)
T * (x, y) = R * T S * T S T MN R T NM , (131)
∂T * ∂x ∂T * ∂y = R * T S * T S T MdN S T dMN R T dNM R T NdM = = R * T S * T F T . ( 132
)
Using these expressions the variational formulation becomes:

Ω R * T S * T n j=1 α j F T D j dΩ + Ω R * T S * T F T E R S dΩ = = m h=1 Ω R * T S * T S T Mb h (y)Na h (x) R T Na h (x)Mb h (y) dΩ, (133)
where

F T D j = dNdN T F j S T MM T G j + NN T F j S T dMdM T G j R T dNdN T F j MM T G j + R T NN T F j dMdM T G j , (134)
F T E = dNdN T S T MM T S NN T RS T dMdM T MM T SR T dNdN T R T NN T RdMdM T . (135)
Again, each integral in Eq. ( 133) can be separated in one integral in the coordinate x and another one in the coordinate y as illustrated in the following expressions:

Ω F T D j dΩ =   x=+L x=-L dNdN T dx F j S T y=+L y=-L MM T dy G j R T x=+L x=-L dNdN T dx F j y=+L y=-L MM T dy G j   + +   x=+L x=-L NN T dx F j S T y=+L y=-L dMdM T dy G j R T x=+L x=-L NN T dx F j y=+L y=-L dMdM T dy G j   , (136)
the integrals of the four submatrices of F T E being given by:

Ω (F T E) 11 dΩ = x=+L x=-L dNdN T dx S T y=+L y=-L MM T dy S, (137)
Ω (F T E) 12 dΩ = x=+L x=-L NN T dx R S T y=+L y=-L dMdM T dy , (138)
Ω (F T E) 21 dΩ = y=+L y=-L MM T dy S R T x=+L x=-L dNdN T dx , (139)
Ω (F T E) 22 dΩ = R T x=+L x=-L NN T dx R y=+L y=-L dMdM T dy , (140)
and the last term in Eq. ( 133)

m h=1   S T y=+L y=-L M b h (y) dy x=+L x=-L N a h (x) dx R T x=+L x=-L N a h (x) dx y=+L y=-L M b h (y) dy   . (141)
Eq. (133) can then be expressed in the compact form:

V 1 (R, S) + K(R, S) R S = V 2 (R, S). (142)
After imposing the vanishing of both functions R and S on their respective domain boundaries, the non-linear system (142) can be solved by applying a Newton strategy, the size of the linearized systems being p + q. As the solution of such problem is not unique, we can impose for example that the norm of S is one by introducing the associated Lagrange multiplier. To accelerate the convergence the first approximated field T is constructed by using the functions a h and b h .

When the solution R and S is obtained (after convergence of the Newton procedure) the new approximation functions can be determined by normalizing:

   F n+1 = R R T ( x=+L x=-L NN T dx ) R G n+1 = S S T ( y=+L y=-L MM T dy ) S . (143)
If the normality of S has been used, then the second expression of Eq. ( 143) reduces to G n+1 = S .

Matrix form of the multidimensional Poisson's problem

We summarize in this section the matrix form in the case of multidimensional Poisson problems.

• Projection stage

In this case Eqs. ( 127) and (128) result

K ij = N k=1 F T ki +L -L dN k dN T k dx k F kj × N l=1 l =k F T li +L -L N l N T l dx l F lj (144) V i = m h=1 N k=1 F T ki +L -L N k a kh (x k )dx k (145)
• Basis enrichment

We assume that matrix K and vectors V 1 and V 2 can be written as:

K =     χ 11 χ 12 χ 21 χ 22 . . . χ N N     , (146)
V 1 =   λ 1 . . . λ N   , (147)
V 2 =   γ 1 . . . γ N   . (148)
In those expressions χ k 1 k 2 is a matrix with dimensions (p k 1 , p k 2 ) and λ k 1 and γ k 1 are vectors of dimension (p k 1 , 1). Thus, is easy to verify that:

                       χ k 1 k 1 = +L -L dN k 1 dN T k 1 dx k 1 N l=1 l =k 1 R T l +L -L N l N T l dx l R l χ k 1 k 2 k 1 =k 2 = +L -L N k 1 N T k 1 dx k 1 R k 1 R T k 2 +L -L dN k 2 dN T k 2 dx k 2 × × N l=1 l =k 1 l =k 2 R T l +L -L N l N T l dx l R l (149) λ k 1 = n j=1 α j                  +L -L dN k 1 dN T k 1 dx k 1 F k 1 j × × N l=1 l =k 1 R T l +L -L N l N T l dx l F lj + + N k 2 =1 k 2 =k 1          +L -L N k 1 N T k 1 dx k 1 F k 1 j × ×R T k 2 +L -L dN k 2 dN T k 2 dx k 2 F k 2 j × × N l=1 l =k 1 l =k 2 R T l +L -L N l N T l dx l F lj                           , (150)
γ k 1 = m h=1     +L -L N k 1 a k 1 h (x k 1 ) dx k 1 × N l=1 l =k 1 R T l +L -L N l a lh (x l ) dx l     (151)
The enrichment functions are then computed by:

F k(n+1) = R k R T k +L -L N k N T k dx k R k . (152)

Building-up a separated representation of the kinetic model of the fiber suspension in complex flows

The Fokker-Planck equation associated to the fiber orientation can be written as a generic advection-diffusion equation whose solution can be searched using a separated representation and the algorithm just described. For the sake of simplicity we assume in what follows a steady sate 2D flow as well as a 2D steady state orientation distribution.

The kinetic theory model related to the fiber suspension was given in Eq. ( 34) that under the steady state assumption writes:

v • ∂Ψ ∂x + E 0 (ϕ, x)Ψ + E 1 (ϕ, x) ∂Ψ ∂ϕ + E 2 (ϕ, x) ∂ 2 Ψ ∂ϕ 2 = 0 (153)
where the components of velocity field depend on the spatial coordinates x, i.e. v T (x) = (u(x), v(x)), and

   E 0 = k( ∂v ∂y -∂u ∂x ) cos(2ϕ) -k( ∂u ∂y + ∂v ∂x ) sin(2ϕ)) E 1 = k 2 ( ∂v ∂y -∂u ∂x ) sin(2ϕ) + 1 2 ( ∂u ∂y + ∂v ∂x )(1 + k cos(2ϕ)) E 2 = -D r (154)
Functions E 0 and E 1 can be written in the general form:

E 0 = e 11 (ϕ
(x) E 1 = f 11 (ϕ) ∂u ∂x (x) + f 12 (ϕ) ∂u ∂y (x) + f 21 (ϕ) ∂v ∂x (x) + f 22 (ϕ) ∂v ∂y (x) (155)
In what follows complex flows will be addressed. In that cases, the solution of the kinetic theory model needs the prescription of the orientation distribution on the inflow boundary. For the sake of simplicity we assume an isotropic orientation distribution on the inflow boundary. Thus, to facilitate the prescription of the boundary conditions one could proceed to define a new variable Ψ, according to:

Ψ(ϕ, x) = Ψ(ϕ, x) -Ψ 0 (ϕ) (156)
where Ψ 0 = 1 2π to account the assumed isotropy. Thus, Eq. ( 153) becomes:

v • ∂ Ψ ∂x + E 0 (ϕ, x) Ψ + E 1 (ϕ, x) ∂ Ψ ∂ϕ + E 2 (ϕ, x) ∂ 2 Ψ ∂ϕ 2 = -E 0 (ϕ, x)Ψ 0 (157)
Now, introducing the separated representation of Ψ:

Ψ(ϕ, x) = i α i F i (ϕ)G i (x) (158)
As previously indicated the construction of this separated representation needs for an iteration procedure, each iteration consisting of two steps:

• Step 1: Projection step. At this step we assume the functions F i and G i known ∀i ∈ [1, ..., n]. We are looking for optimal approximation coefficients α i .

For that purpose, we consider the variational formulation of Eq. (157) that writes:

v• ∂ Ψ ∂x Ψ * + E 0 ΨΨ * + E 1 ∂ Ψ ∂ϕ Ψ * + E 2 ∂ Ψ ∂ϕ ∂Ψ * ∂ϕ = - E 0 Ψ 0 Ψ * (159)
Prior to define the discrete form related to this variational formulation, we must express functions F i et G i in a discrete form. For this purpose we define the vectors N et M containing the shape functions related to the one-dimensional linear piecewise finite element ϕ and x approximations, respectively. Thus, F i and G i can be written:

F i (ϕ) = N T (ϕ)F i G i (x) = M T (x)G i (160)
being F i and G i the vectors containing the nodal values of F i and G i respectively.

Using Eq. (160), Eq. (158) writes:

Ψ = N T F 1 M T G 1 N T F 2 M T G 2 • • • N T F n M T G n     α 1 α 2 . . . α n     = A T α (161)
If we define the vectors dN and dM containing the shape functions derivatives of N and M respectively, the gradient of Ψ results:

∂ Ψ ∂ϕ ∂ Ψ ∂x = dN T F 1 M T G 1 dN T F 2 M T G 2 • • • dN T F n M T G n N T F 1 dM T G 1 N T F 2 dM T G 2 • • • N T F n dM T G n     α 1 α 2 . . . α n     = = B ϕ T B x T α (162)
Now, the weighting functions are approximated as:

Ψ * = α * 1 α * 2 • • • α * n      F T 1 NG T 1 M F T 2 NG T 2 M . . . F T n NG T n M      = α * T A (163)
whose gradient writes:

∂Ψ * ∂ϕ ∂Ψ * ∂x = α * 1 α * 2 • • • α * n      dN T F 1 M T G 1 N T F 1 dM T G 1 dN T F 2 M T G 2 N T F 2 dM T G 2 . . . . . . dN T F n M T G n N T F n dM T G n      = = α * T B ϕ B x (164)
Introducing these expressions in the variational formulation we obtain:

α * T Av T B T x + E 0 AA T + E 1 AB ϕ T + E 2 B ϕ B T ϕ α = = - α * T E 0 Ψ 0 A (165) or α * T Av T B T x + E 0 AA T + E 1 AB ϕ T + E 2 B ϕ B T ϕ K α = = -α * T E 0 Ψ 0 A L ( 166
)
that taking into account the arbitrariness of α * leads to the linear system:

Kα = L (167)
The separated representation of each component of K and L simplifies the integration procedure, which writes:

K ij = F T i Ωϕ NN T dϕ F i × G T i Ωx MvdM T dx G i + +F T i Ωϕ Ne 11 N T dϕ F i × G T i Ωx M ∂v x ∂x M T dx G i + +F T i Ωϕ Ne 12 N T dϕ F i × G T i Ωx M ∂v x ∂y M T dx G i + +F T i Ωϕ Ne 21 N T dϕ F i × G T i Ωx M ∂v y ∂x M T dx G i +
• Step 2: Basis enrichment. Knowing the n functional products, we look for the an additional couple R(ϕ)S(x):

Ψ(ϕ, x) = n i=1 α i F i (ϕ)G i (x) + R(ϕ)S(x) (170)
whose discrete form writes:

Ψ(ϕ, x) = n i=1
N T (ϕ)F i M T (x)G i + N T (ϕ)RM T (x)S (171)

At this stage the unknowns are the new functions R and S. The problem becomes non-linear, needing for an appropriate linearization strategy. We have investigated in our former works the application of the Newton-Raphson algorithm, however a simple linearization based on a kind of alternation directions fixed point strategy seem robust and simple to implement. It consist of alternating the search of S (with R assumed known, in general it comes from the previous iteration) and the search of S using the just computed function R. In that follows we are illustrate the first step that consist of computing R by assuming S known:

Ψ * (ϕ, x) = R * (ϕ)S(x) = N T (ϕ)R * M T (x)S ( 172
)
The gradients of the unknown and weighting functions write:

∂ Ψ ∂ϕ ∂ Ψ ∂x = n i=1 dN T F i M T G i N T F i dM T G i + M T SdN T dM T SN T R (
The resulting linear system is solved under the constraint:

Ωϕ N T R = 0 (176)
that leads to:

α i F i (ϕ) G i (x) (179)
In this case instead of solving a problem in a space of dimension three (the two spatial dimensions x and y and the conformation coordinate ϕ, we look for a separated representation that defines the solution from functional products each one involving a 2D functions (defined in the physical space x) and a one dimensional function defined in the conformation space (the planar orientation in the example here addressed).

From these functional products we can reconstruct the orientation distribution function that allows computing the second order orientation tensor whose eigenvalues and eigenvectors lead to the ellipses depicted in figure 14 that is in perfect agreement with the solution computed using the method of characteristics in the physical space coupled with a discontinuous Galerkin representation in the conformation space depicted in Fig. 5.

Conclusions

In this work the solution of anisotropic fluids in steady recirculating flows (widely encountered in industrial applications) is revisited. The fluid model is described at the microscopic scale (kinetic theory framework) whereas the flow kinematics is modelled at the macroscopic scale (anisotropic Stokes).

Two standard techniques: the first one based on the application of a discontinuous Galerkin discretization and the second one exploiting the periodicity of the solution along the closed streamlines, are analyzed and compared. The second technique seems more accurate (for a given mesh) but it also involves higher computational efforts.

To alleviate the computational efforts two new efficient advanced strategies are described and analyzed. The first technique lies in the use of a model reduction based on the application of the so called proper orthogonal decomposition (or Karhunen-Loève decomposition). It allows significant computing savings for a given accuracy. Moreover, its use allows to define empirical closure approximation leading to accurate macroscopic modelling derived from the microscopic description.

The second advanced technique here described and applied to the solution of complex flows of complex fluids lies in the use of finite sums decompositions (separated representations) of the distribution function involved in the microscopic modelling of the suspension.

These strategies open new perspectives for the treatment of complex flows of complex fluids others than the short fiber suspensions. This natural extension will be addressed in some oncoming works.

ǫ

  = ∞ ǫ = 1 ǫ = 0.3 ǫ = 0.1 ǫ = 0.05 ǫ = 0.03 ǫ = 0

ΩϕN 7 . 4 . 1

 741 Figures 12 and 13 depict the different functions F i (ϕ) and G i (x) related to the driven cavity flow problem previously defined, built by the solution method just described. We can notice that the combination of 6 functional products is enough for accurately represent the searched orientation distribution function Ψ(x, ϕ), i.e.

Figure 2 .

 2 Figure 2. Shear recirculating flow domain: elements where the solution computed with the different numerical techniques are compared.

Figure 4 .FIGURESFigure 5 .

 45 FIGURES

Figure 6 .FIGURESFigure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 14 .

 67891014 FIGURES

  Now, we analyze the effect of moving least square by considering different values of ǫ for D r = 0.1, using the snapshot natural closure. The table that follows summarize the results.

	Diffusion Natural closure Karhunen-Loève	Snapshot	Snapshot-KL
			natural closure natural closure natural closure
	1	2.4 10 -4	2.3 10 -5	2.6 10 -5	2.9 10 -4
	0.1	0.0139	0.0095	0.0091	0.0355
	0.01	0.0046	0.0196	0.0113	0.0471

  table, where we can notice that the empirical snapshot natural closure remains again superior.

	Diffusion Linear Quadratic Hybrid Natural Snapshot
	1	0.0031	0.0704	0.0062	0.0019	0.0008
	0.1	0.1992	0.1516	0.08	0.0353	0.0289
	0.01	1.1286	0.1452	0.1048	0.1216	0.1126
	7 Advanced strategies: techniques based on separated representations

  ) ∂u ∂x (x) + e 12 (ϕ) ∂u ∂y (x) + e 21 (ϕ) ∂v ∂x (x) + e 22 (ϕ) ∂v ∂y

  Nf 11 dN T dϕ F i × G T Nf 12 dN T dϕ F i × G T Nf 21 dN T dϕ F i × G T

	and
	L i = F T i		Ωϕ	Ne 11 Ψ 0 dϕ × G T i		Ωv x	M	∂v x ∂x	dv x +
	+F T i	Ωϕ	Ne 12 Ψ 0 dϕ × G T i	Ωv x	M	∂v x ∂y	dv x +
	+F T i	Ωϕ	Ne 21 Ψ 0 dϕ × G T i	Ωv x	M	∂v y ∂x	dv x +
	+F T i	Ωϕ	Ne 22 Ψ 0 dϕ × G T i	Ωv x	M	∂v y ∂y	dv x
	+F T i	Ωϕ	Ne 22 N T dϕ F i × G T i		Ωx	M	∂v y ∂y	M T dx G i +
	+F T i	Ωϕ				i	Ωx	M	∂v x ∂x	M T dx G i +
	+F T i	Ωϕ				i	Ωx	M	∂v x ∂y	M T dx G i +
	+F T i	Ωϕ				i	Ωx	M	∂v y ∂x	M T dx G i +
	+F T i	Ωϕ	Nf 22 dN T dϕ F i × G T i	Ωx	M	∂v y ∂y	M T dx G i +
	+E 2 F T i		Ωϕ	dNdN T dϕ F i × G T i		Ωx	MM T dx G i	(168)

  Using the previous expressions, the variational formulation results:Ωϕ NN T S T Ωx vMdM T S+ Ωϕ e kl NN T S T Ωx (∇v) kl MM T S+ + 2 Ωϕ f kl NdN T S T Ωx (∇v) kl MM T S+ +E 2 Ωϕ dNdN T S T Ωx MM T S R = = n i=1 α i Ωϕ NN T F i S T Ωx vMdM T G i + Ωϕ e kl NN T F i S T Ωx (∇v) kl MM T G i + +α i Ωϕ f kl NdN T F i S T Ωx (∇v) kl MM T G i + +α i E 2 Ωϕ dNdN T F i S T Ωx MM T G i Ωϕ Ne kl Ψ 0 S T Ωx M(∇v) kl

	and	∂Ψ * ∂ϕ ∂Ψ * ∂x +α i = k=1 M T SdN T 2 l=1 k=1 2 l=1 2 k=1 2 l=1 2 k=1 2 l=1 dM + 2 2 2 -k=1 l=1	173)

T SN T R * (174)

List of Figures

Orientation distribution in a simple shear flow (G = 1) involving fibers characterized by the factor k = 0.8. Dr = 0 (up) and Dr = 0.1 (down). . . . . . [START_REF] Venkiteswaran | A QMC approach for high dimensional Fokker-Planck equations modelling polymeric liquids[END_REF] x y