On the role of standard and anomalous diffusion in suspensions of rigid rods

Dilute suspensions composed of rods are usually described by using the Jeffery’s model that only considers flow-induced orientation. When the concentration increases rods interaction cannot be neglected and the simplest way to take it into account is from a diffusion term that tends to recover an isotropic orientation distribution. However, when considering CNTs suspensions involving large interaction networks, fractional diffusion better describes linear viscoelastic tests. In this work we revisit the fractional diffusion model analyzing its behaviour when applied in nonlinear regimes.


Introduction
It is well known that the process-induced microstructure in short fiber composites and nano-composites determines the mechanical or functional properties of the final part.Thus, the development of accurate models and efficient computational solvers is crucial.Industrial applications usually involve semi-concentrated or concentrated short fiber suspensions in which particles interaction occur and their effects must be added to the flow induced orientation.
Suspensions of particles are accurately modeled by tracking each individual particle in the system.In the dilute regime, the motion of ellipsoidal particles immersed in a Newtonian fluid is accurately described by Jeffery's equation [17].In order to circumvent the difficulties (more computational than conceptual) related to simulations at the microscopic scale where too many particles are present, coarser models were introduced.The recent book [7] gives an overview of multi-scale approaches in computational rheology of rod suspensions.
Suspensions composed of rods have been described at the different scales: microscopic -the scale of the fiber -; mesoscopic -the scale of the population -; and the macroscopic -the scale of the part -.The main difficulties appear when the concentration increases, because interactions require appropriate models to be described.When the concentration is moderate those interactions tend to randomize the particles orientation and consequently can be accurately described by introducing a diffusion mechanism, as originally proposed by Folgar & Tucker [11].The interested reader can refer to many published works focusing in modelling [6] [12] [13] [14] [15], flows [3] [5] [26] and rheology [21] [23].
In this section we revisit the main modeling elements related to the Jeffery's and the Folgar & Tucker models, as well as the difficulties of the last one for describing some experimental rheological tests, that motivated in our former works the introduction of fractional diffusion [4].After revisiting the effects of fractional diffusion in linear viscoelasticity we move to the nonlinear regime to explore the effects of fractional diffusion on the orientation behavior and its associated rheological properties.

Revisiting Jeffery's model
The kinematics of an ellipsoidal particle oriented in the direction given by the unit vector p and immersed in a Newtonian fluid flow characterized by the velocity gradient ∇v is given by the Jeffery equation [17] ṗ where 2D = ∇v + (∇v) T , 2Ω = ∇v − (∇v) T , and the shape factor F = r 2 −1 r 2 +1 depends on the ellipsoid aspect ratio r (major to minor axes ratio).
Mesoscopic kinetic theory models result from the coarsening of microscopic descriptions.In kinetic theory models, the individuality of the particles is lost in favor of a statistical description that substitutes the microscopic entities with a series of conformation coordinates.In the case of suspension of rigid rods, the mesoscopic description consists in giving the fraction of rods that at position x and time t are oriented along direction p, ψ(x, t, p), whose evolution is given by ∂ψ ∂t that despite its apparent simplicity and linear character it suffers from the so-called curse of dimensionality because its multidimensional character.Finally, at the macroscopic scale, the distribution function is substituted with some of its moments [1,3], for example the second and fourth-order moments, a and A, that read respectively where S is the surface of the unit sphere where the orientation vector p is defined.
Here, the level of detail and the involved physics are sacrificed in favour of computational efficiency.The equations governing the time evolution of a are obtained by taking the time derivative of Eq. ( 3), expressing ψ from Eq. ( 2) and then integrating by parts.This procedure yields which, in view of Eq. ( 1), leads to This equation involves the fourth-order moment A. Similarly, the time derivative of the fourth-order moment involves the sixth-order moment A, and so on.Thus, an approximate closure relation is needed in order to express the fourth-order moment A as a function of the lower-order moment a. Different closure relations have been introduced and widely used [2] [10] [19] [25].In what follows we consider the hybrid closure relation, that constitutes a good compromise between implementation simplicity and solution accuracy.It reads: with f = 1 − d d det(a), d being the dimension of the considered space, d = 2 in 2D and d = 3 in 3D.Thus, the hybrid closure combines the linear closure A lin that results exact when the orientation distribution is isotropic, and the quadratic one A qua , at its turn exact in fully aligned suspensions -i.e.ψ(x, t, p) = δ(p − p(x, t))) -.Both are given respectively by and where I is the identity tensor.
For Newtonian suspending fluids, the different concentration regimes (dilute, semi-dilute, semi-concentrated and concentrated) have been extensively analyzed from the modelling, simulation and experimental viewpoints.
In the dilute regime, fibre-fibre interactions are neglected altogether.For the semi-dilute regime, these interactions are usually taken into account in the form of a phenomenological randomizing mechanism, i.e. one adds a diffusion term in the Fokker-Planck equation to obtain where D r is a diffusion coefficient.At the macroscopic scale, this leads to the Folgar & Tucker model [11], Remark 1.In the 2D case the Brownian term in the previous equation must be replaced by −4D r (a − I/2).
The associated constitutive equation reads: with τ given by Remark 2. In the 2D case the Brownian term in the previous equation must be replaced by βD r (a − I/2).

LVE modeling and anomalous diffusion
Linear viscoelastic -LVE -analysis were carried out in [21] to suspensions consisting of functionalized CNTs.LVE involves small amplitude oscillations applied to an essentially isotropic suspension (a iso ≈ I 3 ).In that case the linear closure relation A lin previously introduced is expected to be an accurate approximation for describing A.
We proved in [4] that in the simple shear flow involved in LVE analyses the linear closure approximation yields (A : D) 12 ≈ (A lin (a iso ) : D) 12 = γ 15 , from which the shear stress can be approximated in the general 3D case by leading to the real and imaginary components of the complex modulus and Thus, the loss modulus is found to scale linearly with the frequency ω of the applied oscillation in agreement with the experimental findings [21], however, the storage modulus G is expected scaling at small frequencies with ω 2 instead of ω 0.6 reported in [21], where authors propose the scaling D r ∝ ω p to control the slope of G (ω).
In [8] and [9] the flow-induced bending in the case of non-straight CNTs was proposed as possible mechanism.In other works, the poly-dispersity and the thermally-activated bending were also considered as possible mechanisms.
In [4] authors adopted a more phenomenological viewpoint and considered randomizing mechanism but instead of modeling it from a standard diffusion term, they proposed to consider an anomalous diffusion modelled using fractional derivatives, as described later.

Other anomalous behaviors
Other deviations were noticed when performing other rheological tests like the step strain, when after applying a step strain the stress relaxation results from Eq. ( 13), with the fluid at rest: where again τ = τ 12 and a = a 12 .The orientation evolution results Thus an exponential decay for a and the shear stress τ is expected from the model, however a power-law behaviour is observed experimentally, as reported in [21].

Fractional modeling
Anomalous sub-diffusion was revisited in [16] (the interested reader can also refer to the references therein), justifying the use of non-integer derivatives.In semi-concentrated suspensions the inter-particles interactions could be at the origin of those non-integer derivatives in the diffusion mechanisms.
Brownian diffusion in the configurational space (the one related to the particles orientation) implies ṗ| that leads to the Folgar & Tucker model.Its fractional counterpart reads For additional details on fractional derivatives the interested reader can refer to [18] [24].Now, the time derivative of the second order orientation tensor can be carried out in the usual way using the usual rotary velocity partition where ṗ| J is no more than the Jeffery's contribution whereas the Brownian one makes use of the fractional expression, that introduced into Eq.( 22) and proceeding as described in [4], results da dt with α = 1 − β, d = 2, 3 in 2 and 3 dimensions respectively, and ȧ| J given by Eq. ( 6).

LVE and step strain tests with fractional diffusion
When considering the LVE and step strain tests previous described but now considering the fractional orientation model (26) it results the following expressions of the storage and loss modulus [4] and At small frequencies G scales as ω 1−α = ω β .Thus, it suffices to select β = 0.6 (α = 0.4) to describe the observed experimental behavior.In the case of he stress relaxation after a step strain, with a = a 12 and τ = τ 12 τ = βD r a, where now the orientation time evolution results from When considering α = 0.4 we proved in [4] that both rheological tests fit perfectly with the model predictions.
However, it was never analyzed the impact of such fractional diffusion in nonlinear rheology, that is, in the transient evolution of orientation state when applying a steady state flow.

Nonlinear rheology with fractional diffusion
In this section, and for the sake of simplicity we consider 2D simple shear flows related to and the induced 2D orientation states.The orientation evolution is computed by integrating in time Eq. ( 26) with d = 2 considering different values of α ≤ 1.For α = 0 (β = 1) the results coincide with the ones associated to the Jeffery's model when D r = 0 and to the Folgar&Tucker model when D r = 0.
In this section we consider the simplest closure relation, the quadratic one A ≈ A qua expressed by (9).The initial orientation state is assumed isotropic that corresponds to a(t = 0) = I 2 .

Numerical discretization of the fractional orientation equation
The discrete version of the Grünwald-Letnikov formula [24] [18] for the fractional derivative of order α of function f , t ∈ [a, t], reads: where h is the discrete time step and g α m is given by with Γ(•) the gamma function.Expression (34) is not suitable for large values of m, for example when considering fine dicretizations, because the numerical issues related to the calculation of factorial and gamma functions.Thus, is is preferable using its equivalent form The numerical integration of Eq. ( 32), considering a quadratic closure linearized by considering a semi-implicit integration scheme, writes at time step n: Using the expression of the Grünwald-Letnikov fractional derivative Eq. ( 36) can be written in the discrete form: that accounts for the history of the orientation tensor a weighted by g α m .

Standard diffusion versus no-diffusion
In the present analysis standard diffusion mechanisms are assumed α = 0 (β = 1 and D r = 0) and compared with the Jeffery's solution characterized by the absence of diffusion D r = 0 in a simple shear flow; In this case the flow is characterized by ˙ = 0 and γ = 1.First we consider the standard diffusion model β = 1 (α = 0) for evaluating the impact of the diffusion coefficient D r by considering a null value (no diffusion at all) and D r = 0.
When applying a shear flow the rods are expected orienting on the flow direction, fact that validates the quadratic closure relation here considered.Thus, it is expected that the component a 11 evolves from its initial value a 11 (t = 0) = 0.5 towards a 11 (t → ∞) ≈ 1 when D r = 0 or towards a steady state value a st 11 < 1 (the orientation plateau) for D r = 0.In this simulation we would like to conclude if the introduction of randomizing effects (diffusion) affects or not the orientation process, by advancing or delaying the orientation process.
Figure 1 compares the orientation evolution a 11 (t) for two different values of the diffusion coefficient, D r = 0 and D r = 0.1.It can be noticed that the introduction of diffusion avoids a full alignment, with an orientation plateau a st 11 ≈ 0.85, however, the introduction of diffusion does not alter the orientation process that at the beginning of the orientation process follow the same evolution.
As the apparent viscosity enhancement depends fundamentally on the component a 12 of the orientation tensor, we compare in figure 2 its evolution in time.It is important to note that an overshoot in a 12 is accompanied of a similar overshoot of the apparent viscosity, measurable using adequate rheometers.As noticed in Fig. 2 the diffusion tends to reduce slightly the overshoot.In the limit, when D r → ∞, the isotropic state remains unaltered by the flow and consequently the overshoot disappears.It can be numerically observed that by increasing the diffusion coefficient D r the overshoot monotonically decreases.
In the present analysis fractional diffusion mechanisms are assumed α = 0 (β < 1 and D r = 0) and compared with the Jeffery's solution characterized by the absence of diffusion D r = 0.A simple shear flow is considered again.
In this case the flow is characterized by ˙ = 0 and γ = 1. Figure 3 compares the orientation evolution using the standard orientation model (D r = 0.05, α = 0) and the fractional one with(D r = 0.55, α = 0.5).It can be noticed that the introduction of the fractional diffusion delays the orientation process.Moreover, when focusing on the time evolution of a 12 (t) it can be observed that the use of fractional diffusion suppresses the overshoot.Thus, we can conclude that the use of fractional diffusion greatly affects the orientation process from both the quantitative and qualitative viewpoints.

Conclusions
In this paper we investigated the effects of fractional diffusion in nonlinear rheology involving suspensions composed of rods (short fibers, nano fibers, CNTs).The use of fractional diffusion was proved being a key element for describing the linear viscoelastic behavior of CNTs suspensions, proving the existence of large networks created very probably by the electrostatic effects originated by the CNTs functionalization.
In this paper we explored the impact that such fractional diffusion has when addressing nonlinear viscoelastic flow regimes.The main findings were that the use of fractional sub-diffusion delays the orientation process and at the same time suppresses the apparent viscosity overshoot.
This framework is more general and flexible that the one based on standard diffusion.The use of numerical simulations as the ones considered in [22] should be a valuable route for validating the different modeling scenarios.

Figure 1 .Figure 2 .
Figure 1.Orientation evolution with and without diffusion mechanisms and standard diffusion β = 1 (α = 0) in the case of a simple shear flow characterized by γ = 1.

Figure 3 .
Figure 3. Orientation evolution in a simple shear flow for both standard (broken line) and fractional diffusion (continuous line).