Numerical analysis of the coupling between the flow kinematics and the fiber orientation in eulerian simulations of dilute short fibers suspensions flows

of the flow the fiber the


Numerical Analysis of the Coupling Between the Flow Kinematics and the Fiber Orientation in Eulerian Simulations of Dilute Short Fiber Suspensions Flows
Francisco Chinesta 1 and Arnaud Poitou- Firstly, the momentum balance equations, without inertia and mass terms: N umerical modelling of non-Newtonian flows involves usually the coupling between equations of motion, which define an elliptical problem, and the fluid constitutive equation, which introduces an advection problem related to the fluid history. In short fiber suspensions (SFS) models, the extra-stress tensor depends on the fiber orientation whose evolution can be modeled from a transport problem. In all cases the flow kinematics and the fiber orientation are coupled: the kinematics of the flow governs the fiber orientation, and the presence and orientation of the fibers modify the flow kinematics. For example, in a contraction flow of a dilute suspension, large recirculating areas appear (Lipscomb et aI., 1988).
If one uses SFS flows in material forming processes, the final fiber orientation state depends on the process (induced anisotropy). In this case we need to compute the coupled model (kinematics-fiber orientation) in order to predict the final mechanical properties of the conformed pieces (which depend strongly on the fiber orientation). Of major industrial interest is process optimization to align the fibers in structural pieces in the main stress directions.
where Q is the stress tensor.
Secondly, the incompressibility condition: Di~=Q where y represents the velocity field. (1) This paper focuses on an accurate evaluation of short fibers suspensions models coupling the flow kinematics with the fiber orientation evolution. In coupled models the flow kinematics is usually solved using the finite element method, where the fiber orientation is introduced in the constitutive equation through its value in some points (nodes or integration points). In this paper we will compare in a simple steady shear flow, the exact solutions of the extrastresses associated with the fibers' presence with the numerical simulations obtained using both the method of characteristics and the discontinuous Galerkin's method to solve the equation governing the generalized gradient evolution, in order to avoid the introduction of any closure relation. The error introduced if a quadratic closure relation is considered in the constitutive equation will be also quantified.
Thirdly, the constitutive equation for a dilute suspension of high aspect-ratio particles is given by: From a physical point of view, we can consider that the eigenvalues of the second-order orientation tensor (g) represent the probability of finding the fiber in the direction of the corresponding eigenvectors.
And lastly, if we consider spheroidal fibers, we can describe the orientation evolution with the Jeffery equation: where g is the vorticity tensor, k is a constant that depends on the fiber aspect ratio r (fiber length to fiber diameter ratio): where p denotes the pressure, 1 the unit tensor, " the viscosity (which depends on the chosen model as explained in Meslin, 1999), g the strain rate tensor, N p a scalar parameter depending on both the fiber concentration and the fiber aspect ratio, ":" the double dyadic product, i.e. [~: glij = 0ijkPkl and~the fourth-order orientation tensor defined by: (4) k=,l-l/,l+l A solution of this equation is given by: where P is the unit vector aligned in the fiber axis direction, ® denotes the dyadic product (l.e, (.e ® .e)ij = P'pj)' and \jf(.e) is the orientation distribution function satisfying the normality condition: (12) where .eo is the initial orientation Po= .e (t =0) and the transport is the solution of the following transport equation: The evolution equation of the second-order orientation tensor results: If \jf(p) = o(p -13), with 00the Dirac's distribution, all the orientation probabilitY is concentrated in the direction defined byp, and the corresponding orientation tensor will be tJ. =P® P® P-® p.
We can also define the second order orien1ation tensor as: It is easyto verify that if \jf(p) = o(pp), the fourth-order orientation tensor can be written as the tensorial product of the second order orientation tensor, i.e.: and then, the second order orientation tensor related to a planar isotropic orientation state is: For general expressions of \jf(p) the previous relation is not exact, and in that case Equation (7) becomes a closure approximation known as quadratic closure relation. However, other closure relations are usually applied (Advani and Tucker, 1990).
In a planar case, the isotropic orientation state is defined by the uniform distribution function: (8) A similar equation can be derived for the evolution of the fourth-order orientation tensor, which in this case involves the sixth-order orientation tensor.
The solution of Equation (14) with a quadratic closure relation (~=~®~) can be written as: where~o is the initial second order orientation tensor Qo = Q(t = 0).
Remark 1: Slender body theories (Batchelor, 1971;Dmh and Armstrong, 1984) assume that the fibers have an infinite aspect ratio, i.e. k = 1 and the suspension viscosity is equal to the solvent viscosity. In that case, to describe the increase of the shear viscosity it is necessary to introduce diffusion effects to avoid the full alignment of the fibers with the flow (Folgar and Tucker, 1984). The orientation equation becomes: The introduction of the diffusion term Dr (~-1/2) seems no longer necessary if we consider the spheroidal particles model because in that case the fact that k *' 1 induces a continuous rotation of the fibers which can describe the rise in shearviscosity (Meslin, 1999).
In the numerical simulation of dilute fiber suspensions there are three kinds of numerical models: model is the necessity of proceeding with a fully Eulerian description, whereas the fixed point coupled models can proceed with a Eulerian-Lagrangian mixed description, using, as previously described, the method of characteristics (the most accurate technique solving advection equations) in the solution of the orientation equation.
from which the expression for the second order orientation tensor Q results: 2. Coupled models using the Evans hypothesis solve the flow kinematics assuming that at each point of the fluid domain the fibers are aligned with the flow (local alignment hypothesis) (Lipscomb et al., 1988;Chiba et al., 1990). In this case, the flow model is given by Equations (1), (2) and (3), where the fiber orientation is given by: 1. Uncoupled models solve the kinematics without considering the presence of fibers. In this case, the equations of motion result in the Stokes's problem defined by Equations (1), (2) and (3) with N p =O. From the velocity field obtained by the solution of the Stokes's problem, we can compute the fiber orientation evolution solving any of the models described previously. Givler et al. (1983) solved the Jeffery equation along the flow streamlines. To take into account different initial orientation distributions, a great number of fibers with different initial orientation are considered, and the averaging of their orientations at each point allows us to compute the different orientation tensors.
The first type of model does not allow prediction of the strong variations in the flow patterns that appear for example in contraction or expansion flows, even for dilute regimes (Lipscomb et al., 1988). The second type of model, extensively used, allows us to predict the existence of recirculating flow areas near the contraction corners (Lipscomb et aL, 1988;Chiba et al., 1990). Moreover, as observed in several works (Towsend and Walters, 1993), and theoretically proven in Poitou et al. (2000) and Chinesta and Chaidron (2001), for fibers with infinite aspect ratio the local alignement of the fibers with the flow is the only solution in general recirculating flows. However, the local alignment cannot be applied either to divergent flows for which this solution is not stable, or to fibers with finite aspect ratio whose orientation changes continuously along their trajectory. As expected, the best modeling, widely used, consists of coupling flow kinematics and fiber orientation, as described for the third type of models. Different numerical strategies are used to solve the orientation equation: the method of characteristics, the SUPG or the discontinuous finite element methods or the discontinuous finite volumes method, using the second-or the fourth-order orientation tensors, as well asthe generalized gradient~asvariables. One of the main difficulties in using Eulerian discretisation techniques to compute the evolution of the tensorial orientation fields is the low accuracy of standard finite element interpolations of these orientation tensors, as discussed in Chinesta et al. (1999) and Chinesta et aL (2000).
On the other hand, the ellipticity of the equations of motion, when the orientation field is taken at the previous iteration, allows us to use a standard Galerkin formulation in finite elements. In general, this technique is used in most of the related works. Thus, in each element of the finite element mesh, the fiber orientation is considered only at some points (integration points). When the fibers are assumed of spheroidal shape and finite aspect ratio, it is well known that in a simple shear flow they cannot be fully aligned in the flow direction, rotating continuously. The higher the fiber aspect ratio the longer the fibers remain close to the flow direction, and faster is their rotation. We can expect that considering only the orientation state in some points we can predict an orientation state, which does not correspond with the physical reality. Thus, the probability of finding, at the integration points, the fibers aligned in the flow direction increases with the fibers' aspect ratio. In this case Np'~:Q .. 0 and the numerical model cannot predict the rise in shtiar viscosity experimentally observed. An alternative way to avoid this over-alignment through modeling is to consider the slender body theory, where, in spite of the fact that the fibers are assumed to have an infinite aspect ratio, the diffusion term introduced in Equation (16) avoids the full alignment of the fibers with the flow.
This paper evaluates the accuracy of coupled models using the spheroidal particles theory. We will analyse a simple shear flow, for which we can derive theoretically the rise in shear viscosity. These results will be compared with the different numerical simulations, without using any closure relation, in order to derive a conclusion about their accuracy, 3. Coupled models take into account both the dependence of the kinematics on the fiber orientation and the orientation induced by the flow kinematics. Usually the coupled models are solved by means of a fixed point strategy. In this case, at each iteration the flow kinematics results from Equations (1), (2) and (3) assuming the fiber orientation field at the previous iteration. From the kinematics just computed, the fiber orientation is updated solvlnq the advection equation governing the evolution of the second-or the fourth-order orientation tensors, but in these cases the consideration of a closure relation is required. These advection equations can be integrated by using any accurate numerical technique for hyperbolic equations: the method of characteristics, SUPG or discontinuous finite element techniques, discontinuous finite volumes, ... (Rosenberg et al., 1990;Ausias, 1991;Altan et aL, 1992;Souloumiac, 1996;Azaiez et al., 1997;Chiba and Nakamura, 1998;Poitou et al., 2000). Coupled models solving simultaneously the flow kinematics and the fiber orientation (fully coupled models) are rare in literature. The main difficulty of using fully coupled models is the different character of the model equations, which requires specific numerical techniques. Another particularity of this kind of To take into account the effects of the fiber orientation on the flow kinematics, an accurate description of the fibers' orientation state in the equations of motion is required. This paper focuses on this area. The aim of this work is not the computation of coupled solutions, because, as commented previously, these calculations do not introduce additional difficulties if both an accurate strategy for solving advection equations and an accurate description of the orientation field into the equations of motion are performed. For this reason, only simple shear flows are considered in order to evaluate the capabilities of different numerical strategies to introduce accurately the fiber orientation distribution into the motion equations. In spite of the simplicity of such flows, the conclusions can be generalized to more general flows, where an effective coupling between the orientation and motion equations must be carried out for simulating industrial applications.
In the same way, other effects such as the fiber-fiber or wall-fiber interactions, which can be important in semi-concentrated or concentrated flows regimes, are not considered in this work, because to introduce accurately these effects more sophisticated mechanical models are required. In a first approximation, the introduction of diffusion effects in the equation governing the evolution of the fiber orientation distribution could describe the fiber-fiber interaction in semi-concentrated flows. In this case, the fiber orientation field is smoother than the one obtained without considering diffusion effects, and consequently, it can be accurately introduced into the equations of motion to compute coupled solutions.
Finally, this work concerns only planar orientation distributions. However, the analysis proposed in this paper can be applied to 3D orientation distributions as well as to general 3D flows.

Orientation Solver
In this section we will consider the velocity field as a datum. As just mentioned, we will use a numerical strategy to compute the fiber orientation field without the introduction of any closure relation. For this purpose, we propose to solve the evolution equation related to the generalized gradient fk: where fl is the unit outwards vector defined on the domain boundary an. Now, we must integrate the orientation equation along the characteristics(streamlines for steady flows) from the boundary condition fk~=fk (~E r -.J =J until we return to the starting point~. A fourth-order Runge-Kutta scheme with a control step has been used in the numerical integration, which leadsto very accurate solutions. The main limitation of this strategy is the necessity of reconstructing each trajectory associated with each point where the orientation solution must be computed. In general, this procedure is highly time consuming for the CPU. A semi-Lagrangian strategy that combines an integration by characteristics into each element and a Lagrangian interpolation on the mesh skeleton, reducing drastically the CPU time, was proposed by the authors in Chinesta et al. (2000). Another fully Eulerian strategy applies the discontinuous Galerkin method to solve Equation (19), as we will describe in the next section.

A First-order Discontinuous Finite Element Discretisation
The simplest discontinuous Galerkin formulation lies in taking a value of fk constant into each element of a finite element mesh. We will denote by g(t) the value of the generalized gradient fk in the element n e at time t. The initial condition is given by g(t = 0) =1, VeE[l, ... ,N e ], where N e is the number of elements in the mesh, i.e, n = u~~~ene. The integration of an advection equation requires also a boundary condition on the inflow boundary 1_, which in our case results in fk~e r -.J = J.
The first-order time discretization is given by: ve e [1, ..., N e ] where an e is the element boundary and 11 is the unit-outwards vector defined on the element boundary. Asg is not defined on the element boundaries, we take on the inflow boundary an:, The main advantage of this explicit strategy is a significant reduction in the CPU time.
(32) (31 ) To obtain a from 8 we use the Jefferyequation: where the symmetry relation 01212 = 01212 has been introduced, and A' J1 represents the rise in shear viscosity, which depends on°1

212'
Remark 2: When the fibers have a quasi-infinite aspect ratio, k = 1, the steady solution is given by the local alignment of the fiber with the flow, i.e. 'II fJt" a) =l>(A, 8 -0), from which it results 01212 = 0, and consequently A' J1 =O.

Computin9 the Fiber Orientation Tensors and the Rise in Shear Viscosity
Computing the Orientation Tensors from the Generalized Gradient With the value of &e computed at different points with the method of characteristics, or in each element using the justdescribed first-order discontinuous Galerkin method, the expression of the orientation tensors can be computed applying the numerical procedure described in this section.
With the orientation distribution on the inflow boundary 'II oC4e r _, e) done, the expression of f! and f! can be obtained from their definition expressions: = from the upstream element~and q+ is the flux leaving the element {le. Taking into account the incompressibility, we obtain: We will be able to evaluate the previous expressions as soon as the fiber orientation distribution at any point~, '¥~, p), is computed from the initial one, '11 0 (known on the inflow boundary), applying the following balance equation: where a defines the fiber orientation at point~, which, when the fiber reaches point~(both~and~are on the same trajectory), its orientation is defined by the angle 8. Moreover, all the fibers whose initial orientation at point~is in the interval [a, a + 00.] will be in the oriented interval [8, 8 +de] when they reach the position~•

1-k 2
Analysis of the Solution Along a Streamline for Fibers Aligned in the Flow Direction on the Inflow Boundary In this case, at each point x, the fibers will be perfectly aligned in the direction given by the Jefferyequation, i.e.: From this equation we can come to the conclusion about the periodicity in the fibers' orientation. The spatial fiber orientation period X is in our example X = 41t/~= 40.24. We will considerthe inflow boundary located at x = 0, where in a first time, the fibers are aligned with the flow. Later, an isotropic fiber orientation on the inflow boundary will be considered. (30) The  (19), an average of°1 212 will be obtained implicitly. In this situation we can define: The average value in the spatial period X can be considered, using the ergodic principle, as the exact steady local solution: where 01212 is the value found in the element n e , which is assumed to be constant into the element. The numerical application gives us the value aftfl = 0.058 very close again to the exact solution.
The main advantage of using a discontinuous finite element method to solve the equation governing the evolution of~, Equation (19), is that it does not require the resolution of any linear system, and moreover, it reduces drastically the CPUtime with respect to an integration by characteristics.
where 0l c 1q1denotes the value of 01212 at the center of gravity of the element of.
When the number of elements increases, the numerical solution aftf.2 converges to the exact one, Le. aftf.2~a,e~rtf = 0.059 (in our numerical example) when N e~0 0.

Analysis of the Solution in the Whole Domain Using an Eulerian Discretisation
In this section we consider the 2D flow model defined by the following kinematics: which is far from the exact value. Until now the generalized gradient~has been obtained exactly from Equation (33) at the center of gravity of each which is very close to the exact value. Now, we can evaluate°1 2 at the center of gravity of each element 0f¥,i as well as the mean value of 01212 assuming a quadratic closure relation a,ntn qcr:

Analysis of the Solution Along a Streamline for an Isotropic Fiber Orientation Distribution on the Inflow Boundary
In this case, the variation of°1 212 is smoother than the previous one, because even if the fibers rotate very quickly, not all the fibers rotate at the same point, due to their different initial alignments. Thus, for the mesh considered in the previous example, we obtain: If we taker =1, then the period X of~depends on the y coordinate in the following form: (42)

X(y)=Ñ
ow, if we define a uniform mesh in the fluid domain [0, X(y = 1)] x [0,1], the characteristic element size near the upper wall, y = 1, will be lower than the fiber orientation period X(y), whereas in the area close to the lower boundary, y =0, the characteristic element size is greater than the fiber orientation period.
For this scenario, very close to actual general flow situations, we would like to compare results obtained by an accurate computation of tensor~at the center of gravity of the mesh elements (by means of the method of characteristics, for example) with the one obtained solving the evolution equation of~, Equation (19), by means of the discontinuous Galerkin method (which assumes a value of~constant into each element, as described previously).
Using a mesh of 940 elements, N e = 940, and imposing the local alignment of the fibers with the flow on the inflow boundary, i.e. eti =(1,0), we obtain the following results: Now, if we assume on the inflow boundary an isotropic fiber orientation, i.e, '¥~E r -> =1/21t we obtain: N.
--""",-,212-' Nomenclature accuracy is not significantly degraded when a considerable reduction in the number of degrees of freedom (in the order of 10) is performed. Even if the result obtained with the integration by characteristics seems slightly better, this technique requires much more CPU time than the second numerical technique. In order to improve the accuracy of the finite element technique, a local remeshing is needed. As in this technique we use a discontinuous approximation of the unknown field~, we can proceed to non-structured local remeshings.

Conclusions
In this paper we have shown that in spite of the fast rotation of the fibers with finite aspect ratio in shear flows, the extrastresses are correctly evaluated in considering the fiber orientation only in a reduced number of points. If we use uniform meshes, the characteristic element length will be, in general, greater or smaller than the fiber rotation period. Even in that case, the solution accuracy is excellent, mainly when an isotropic fiber distribution is assumed on the inflow boundary.
We have also noticed that very accurate solutions may be obtained by using the discontinuous Galerkin method, and they open new perspectives in the use of fully Eulerian coupled techniques. The highest deviations were found when the extrastress tensor was computed using a quadratic approximation for the fourth-order orientation tensor.

Other Symbols
I8i dyadic product.