Kinetic Theory Microstructure Modeling in Concentrated Suspensions

Abstract

When suspensions involving rigid rods become too concentrated, standard dilute theories fail to describe their behavior. Rich microstructures involving complex clusters are observed, and no model allows describing its kinematics and rheological effects. In previous works the authors propose a first attempt to describe such clusters from a micromechanical model, but neither its validity nor the rheological effects were addressed. Later, authors applied this model for fitting the rheological measurements in concentrated suspensions of carbon nanotubes (CNTs) by assuming a rheo-thinning behavior at the constitutive law level. However, three major issues were never addressed until now: (i) the validation of the micromechanical model by direct numerical simulation; (ii) the establishment of a general enough multi-scale kinetic theory description, taking into account interaction, diffusion and elastic effects; and (iii) proposing a numerical technique able to solve the kinetic theory description. This paper focuses on these three major issues, proving the validity of the micromechanical model, establishing a multi-scale kinetic theory description and, then, solving it by using an advanced and efficient separated representation of the cluster distribution function. These three aspects, never until now addressed in the past, constitute the main originality and the major contribution of the present paper.

1. Introduction

Short fibers, nanofibers or carbon nanotubes—CNT—suspensions present different morphologies, depending on their concentrations. When the concentration is dilute enough, one can describe the microstructure by tracking a population of rods that move with the suspending fluid and orient depending on the velocity gradient according to the Jeffrey’s equation [1] that relates the orientation evolution with the flow velocity field. In that case, the motion and orientation of each fiber is assumed decoupled from the others. The representative population can involve too many fibers, and in that case, the computational efforts to track the population is unaffordable. Thus, the simple and well-defined physics must be sacrificed in order to derive coarser descriptions.

Kinetic theory approaches [2,3,4] describe such systems at the mesoscopic scale. Their main advantage is their capability to address macroscopic systems, while keeping the fine physics through a number of conformational coordinates introduced for describing the microstructure and its time evolution. At this mesoscopic scale, the microstructure is defined from a distribution function that depends on the physical space, the time and a number of conformational coordinates—the rods orientation in the case of slender body suspensions. The moments of this distribution constitute a coarser description, in general, used in macroscopic modeling [5].

At the macroscopic scale, the equations governing the time evolution of these moments usually involve closure approximations, whose impact in the results is unpredictable [6,7,8]. Other times, macroscopic equations are carefully postulated in order to guarantee the model objectivity and its thermodynamical admissibility.

In the case of dilute suspensions of short fibers, the three scales have been extensively considered to model the associated systems without major difficulties. However, as soon as the concentration increases, the difficulties appear. In the semi-dilute and semi-concentrated regimes, fiber-fiber interactions occur, but in general, they can be accurately modeled by introducing a sort of randomizing diffusion term [9]. There is a wide literature concerning dilute and semi-dilute suspensions, addressing modeling [10,11,12,13,14], flows [15,16,17,18] and rheology [19,20]. These models describe quite well the experimental observations. When the concentration increases, rods interactions cannot be neglected anymore and appropriate models addressing these intense interactions must be considered, as, for example, the one proposed in [21]. Recent experiments suggest that short fibers in concentrated suspensions align more slowly as a function of strain than models based on Jeffery’s equation predict [22]. For addressing this issue, Wang et al. [22] proposed the use of a strain reduction factor; however, this solution violates objectivity. Later, the same authors proposed an objective model by decoupling the time evolution of both the eigenvalues and the eigenvectors of the second order orientation tensor [23]. In [24], an anisotropic rotary diffusion is proposed for accounting for the fiber-fiber interactions, and the model parameters were selected by matching the experimental steady-state orientation in simple shear flow and by requiring stable steady states and physically realizable solutions.

The worst scenario concerns concentrated suspensions involving entangled clusters exhibiting aggregation/disaggregation mechanisms. A first approach in that sense was proposed in [25]. The first natural question is how to describe such systems? At the macroscopic scale, one could try to fit some power-law constitutive equation; however, this description does not allow one to describe the microstructure. At the microscopic scale, direct numerical simulations describing complex fiber-fiber interactions can be carried out in small enough representative volumes [26,27]. A natural candidate to be a reasonable compromise between (fine) micro and (fast) macro descriptions consists of considering, again, a kinetic theory description.

The main issue of such an approach lies in the fact that it must include two scales, the one involving the aggregates and the one related to the rods composing the aggregates. What are the appropriate conformational coordinates? How does one determine the time evolution of these conformational coordinates? How does one represent, simultaneously, both scales, the one related to the aggregates and the other related to the fibers? How does one derive the interaction mechanisms?... These questions will be addressed in this work.

In [28], the authors propose a first attempt to describe such clusters from a micromechanical model, but neither its validity nor the rheological effects were addressed. Later, in [29], the authors applied this model for fitting the rheological measurements in concentrated suspensions of carbon nanotubes (CNTs) by assuming a rheo-thinning behavior at the constitutive law level. However, three major issues were never addressed until now: (i) the validation of the micromechanical model by direct numerical simulation; (ii) the establishment of a general enough multi-scale kinetic theory description, taking into account interaction, diffusion and elastic effects; and (iii) propose a numerical technique able to solve the kinetic theory description. This paper focuses on these three major issues, proving the validity of the micromechanical model, establishing a multi-scale kinetic theory description and, then, solving it by using an advanced and efficient separated representation of the cluster distribution function. These three aspects, never until now addressed in the past, constitute the main originality and the major contribution of the present paper.

In this paper, we first propose two micro-mechanical models for describing the kinematics of rigid and deformable aggregates in Section 2 and Section 3, respectively. The double scale kinetic theory model will be, then, addressed in Section 4. Section 5 considers different interaction mechanics in order to derive an extended kinetic theory model. The issues related to its numerical solution will be deeply analyzed in Section 6. Section 7 presents some verification and validation tests of the proposed micro-mechanical models, and finally, first numerical results are given in Section 8.

2. Micro-Mechanical Description of the Kinematics of Rigid Clusters

As just discussed, when the rod concentration is high enough, rod interactions cannot be neglected anymore, and appropriate models addressing these intense interactions must be addressed. In the extreme case, a sort of cluster composed of entangled rods is observed [25,30]. When the suspension flows, those clusters seem animated by almost rigid motions.

Clusters can be viewed as entangled aggregates of rods, with these rods consisting of a rigid segment joining the two beads located at the segment ends, where it is assumed that forces apply. Microscopic simulations become unaffordable, due to the computational resources needed to address scenarios of real interest. On the other hand, macroscopic modeling is without any doubt too phenomenological for describing the main flow and microstructure features. Kinetic theory approaches could be an appealing compromise between accuracy and computational efficiency.

In what follows, we consider the cluster idealization depicted in Figure 1 in the case of a hypothetical 2D aggregate that does not involve size effects, but at least involves the orientation distribution. Obviously, this representation is a crude simplification, but it constitutes the simplest model and could be later enriched, without major difficulties, to include more complex and richer configurations.

Figure 1. Star representation of a cluster composed of rigid rods.

Figure 1. Star representation of a cluster composed of rigid rods.

First, we consider a 2D rigid cluster consisting of N rods of length, $2L$, oriented in the directions, ${\mathbf{p}}_i$, $i=1,\cdots ,N$, as sketched in Figure 2. Brownian effects can be neglected, and then, only flow induced hydrodynamic forces must be considered (the rod-rod hydrodynamic interactions are not considered in this work). Forces, ${\mathbf{F}}_i$, apply on each bead located at position, $L{\mathbf{p}}_i$ (Figure 2), and are proportional to the difference of velocity between the one flow unperturbed by the presence of the cluster at the bead location and the other one, the bead itself:

$${\mathbf{F}}_i=\xi (\nabla \mathbf{v}·{\mathbf{p}}_{i}L-{\dot{\mathbf{p}}}_{i}L)$$

(1)

where ξ is the friction coefficient and $\mathbf{v}$ the flow velocity field unperturbed by the cluster presence. We can notice that forces are, by construction, self-equilibrated.

Figure 2. Hydrodynamic forces applied on a rigid cluster composed of rods.

Figure 2. Hydrodynamic forces applied on a rigid cluster composed of rods.

Remark: In this paper, we consider the following tensor products:

• if $\mathbf{a}$ and $\mathbf{b}$ are first rank tensors, then the single contraction, “·”, reads $(\mathbf{a}·\mathbf{b})=a_j{b}_j$ (Einstein summation convention);
• if $\mathbf{a}$ and $\mathbf{b}$ are first rank tensors, then the dyadic product,“⊗”, reads ${(\mathbf{a}\otimes \mathbf{b})}_{jk}=a_j{b}_k$;
• if $\mathbf{a}$ and $\mathbf{b}$ are first rank tensors, then the cross product, “×”, reads ${(\mathbf{a}\times \mathbf{b})}_j={ϵ}_{jmn}{a}_m{b}_n$ (Einstein summation convention) with ${ϵ}_{jmn}$, the components of the Levi-Civita tensor, $\mathbf{ϵ}$ (also known as the permutation tensor);
• if $\mathbf{a}$ and $\mathbf{b}$ are, respectively, second and first rank tensors, then the single contraction, “·” reads ${(\mathbf{a}·\mathbf{b})}_j=a_{jm}{b}_m$ (Einstein summation convention);
• if $\mathbf{a}$ and $\mathbf{b}$ are second rank tensors, then the single contraction, “·”, reads ${(\mathbf{a}·\mathbf{b})}_{jk}=a_{jm}{b}_{mk}$ (Einstein summation convention);
• if $\mathbf{a}$ and $\mathbf{b}$ are second rank tensors, then the double contraction, “:” reads $(\mathbf{a}:\mathbf{b})=a_{jk}{b}_{kj}$ (Einstein summation convention).

The moment induced by forces applied on rod i is given by:

$${\mathbf{M}}_i=2L{\mathbf{p}}_i\times {\mathbf{F}}_i$$

(2)

By neglecting inertia effects, the resulting moment for the whole cluster must vanish:

$$\sum _{i=1}^{i=N}{\mathbf{M}}_i=\mathbf{0}$$

(3)

then, taking into account Equations (1) and (2) results in:

$$\sum _{i=1}^{i=N}{\mathbf{p}}_i\times (\nabla \mathbf{v}·{\mathbf{p}}_i)=\sum _{i=1}^{i=N}{\mathbf{p}}_i\times {\dot{\mathbf{p}}}_i$$

(4)

If we define the cluster angular velocity, $\mathbf{\omega }$, such that:

$${\dot{\mathbf{p}}}_i=\mathbf{\omega }\times {\mathbf{p}}_i$$

(5)

in the 2D case, where $\mathbf{\omega }$ and ${\mathbf{p}}_i$ are orthogonal, ${\mathbf{p}}_i\times {\dot{\mathbf{p}}}_i=\mathbf{\omega },\forall i$, and then, it results that:

$$\sum _{i=1}^{i=N}{\mathbf{p}}_i\times (\nabla \mathbf{v}·{\mathbf{p}}_i)=\sum _{i=1}^{i=N}{\mathbf{p}}_i\times {\dot{\mathbf{p}}}_i=N\mathbf{\omega }$$

(6)

Thus, the kinematics of the rigid 2D cluster can be defined from:

$$\mathbf{\omega }=\frac{{\sum }_{i=1}^{i=N}{\mathbf{p}}_i\times (\nabla \mathbf{v}·{\mathbf{p}}_i)}{N}$$

(7)

the angular velocity of each rod being j:

$${\dot{\mathbf{p}}}_j=\mathbf{\omega }\times {\mathbf{p}}_j=\frac{\left({\sum }_{i=1}^{i=N}{\mathbf{p}}_i\times (\nabla \mathbf{v}·{\mathbf{p}}_i\right)\times {\mathbf{p}}_j}{N}$$

(8)

when the cluster contains many rods, an alternative description consists of defining its orientation distribution, $\psi \left(\mathbf{p}\right)$. In this case, all the sums in the previous expressions must be substituted by the corresponding integrals weighted with the distribution function, ψ. Thus, in the 2D case, Equation (7) leads to:

$$\mathbf{\omega }={\int }_{\mathcal{S}(0,1)}\left(\mathbf{p}\times (\nabla \mathbf{v}·\mathbf{p})\right)\psi \left(\mathbf{p}\right)d\mathbf{p}$$

(9)

where $S(0,1)$ denotes the unit sphere surface. The angular velocity of any rod, $\mathbf{q}$, reads:

$$\dot{\mathbf{q}}=\mathbf{\omega }\times \mathbf{q}={\int }_{\mathcal{S}(0,1)}\left(\left(\mathbf{p}\times (\nabla \mathbf{v}·\mathbf{p})\right)\times \mathbf{q}\right)\psi \left(\mathbf{p}\right)d\mathbf{p}$$

(10)

Cross products can be written in a more compact form by using the third order Levi-Civita permutation tensor, $\mathbf{ϵ}$, that allows writing using the notation previously introduced, ${(\mathbf{u}\times \mathbf{v})}_i={ϵ}_{ijk}{u}_j{v}_k$ (Einstein summation convention) or $(\mathbf{u}\times \mathbf{v})=\mathbf{ϵ}:(\mathbf{v}\otimes \mathbf{u})$. Thus, Equation (9) writes:

$$\mathbf{\omega }=\mathbf{ϵ}:\left(\nabla \mathbf{v}·\mathbf{a}\right)$$

(11)

with $\mathbf{a}$, the second moment of the orientation distribution, $\psi \left(\mathbf{p}\right)$:

$$\mathbf{a}={\int }_{\mathcal{S}(0,1)}\mathbf{p}\otimes \mathbf{p}\psi \left(\mathbf{p}\right)d\mathbf{p}$$

(12)

where, for the sake of notational simplicity, we do not explicate the time dependence of the distribution function, ψ.

On the other hand, we can define the rotation tensor, $\mathbf{W}$, such that:

$$\mathbf{\omega }\times \mathbf{p}=\mathbf{W}·\mathbf{p}$$

(13)

that allows writing:

$$\dot{\mathbf{p}}=\mathbf{W}·\mathbf{p}$$

(14)

It is easy to prove from Equations (11)–(14) [28] that:

$$\mathbf{W}=\Omega +\mathbf{D}·\mathbf{a}-\mathbf{a}·\mathbf{D}$$

(15)

where $\mathbf{D}$ and Ω are, respectively, the symmetric and skew-symmetric components of the velocity gradient, $\nabla \mathbf{v}$.

Thus, a rigid cluster rotates with a velocity that only depends on the second moment of its orientation distribution, $\mathbf{a}$. Any rigid cluster having the same orientation tensor, $\mathbf{a}$, will have the same angular velocity. In this case, we should derive an expression given the time evolution of $\mathbf{a}$, $\dot{\mathbf{a}}$. For this purpose, we start from the definition of $\mathbf{a}$ in Equation (12) and consider its time derivative:

$$\dot{\mathbf{a}}={\int }_{\mathcal{S}(0,1)}\left(\dot{\mathbf{p}}\otimes \mathbf{p}+\mathbf{p}\otimes \dot{\mathbf{p}}\right)\psi \left(\mathbf{p}\right)d\mathbf{p}$$

(16)

where the fact that the time derivative of ψ vanishing is used. Now, taking into account Equation (14), it results that:

$$\dot{\mathbf{a}}={\int }_{\mathcal{S}(0,1)}\left((\mathbf{W}·\mathbf{p})\otimes \mathbf{p}+\mathbf{p}\otimes (\mathbf{W}·\mathbf{p})\right)\psi \left(\mathbf{p}\right)d\mathbf{p}=\mathbf{W}·\mathbf{a}+\mathbf{a}·{\mathbf{W}}^T$$

(17)

It is easy to prove the objectivity of expression (17) that corresponds to the Jauma’s convective derivative, as well as its null trace, i.e., $\text{tr}\left(\dot{\mathbf{a}}\right)=0$ [28].

3. Micro-Mechanical Description of the Kinematics of Deformable Clusters

In this section, we are considering a more realistic scenario consisting in deformable clusters. We consider two types of forces applying on each bead, one due to the fluid-rod friction, once more modeled from:

$${\mathbf{F}}_i^H=\xi (\nabla \mathbf{v}·{\mathbf{p}}_{i}L-{\dot{\mathbf{p}}}_{i}L)$$

(18)

where the superscript, $“H”$, refers to its hydrodynamic nature and the other ${\mathbf{F}}_i^C$, due to the rods entanglements. This last force is assumed scaling with the difference between the rigid motion velocity (the one that the bead would have if the cluster would be rigid) and the real one:

$${\mathbf{F}}_i^C=\mu (\mathbf{W}·{\mathbf{p}}_{i}L-{\dot{\mathbf{p}}}_{i}L)$$

(19)

when μ is a large enough force, ${\mathbf{F}}_i^C$ dominate the momentum balances, enforcing the cluster rotary velocity, $\mathbf{W}$, to each rod composing it (rigid behavior).

By adding both forces, it results that:

$${\mathbf{F}}_i={\mathbf{F}}_i^H+{\mathbf{F}}_i^C=L((\mu \mathbf{W}+\xi \nabla \mathbf{v})·{\mathbf{p}}_i-(\xi +\mu ){\dot{\mathbf{p}}}_i)$$

(20)

which can be rewritten as:

$${\mathbf{F}}_i=L(\xi +\mu )\left(\frac{\mu \mathbf{W}+\xi \nabla \mathbf{v}}{\xi +\mu }·{\mathbf{p}}_i-{\dot{\mathbf{p}}}_i\right)$$

(21)

Now, if we define the equivalent traceless gradient, $\mathbf{G}$:

$$\mathbf{G}=\frac{\mu \mathbf{W}+\xi \nabla \mathbf{v}}{\xi +\mu }$$

(22)

We can generalize all the results obtained for suspensions of dilute rods, by considering an equivalent friction coefficient, $\xi +\mu$, and an equivalent gradient of velocities, $\mathbf{G}$. In particular, the rod rotary velocity results in:

$${\dot{\mathbf{p}}}_i=\mathbf{G}·{\mathbf{p}}_i-\left(\mathbf{G}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)\right){\mathbf{p}}_i$$

(23)

The term $\mathbf{G}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)$ deserves some additional comment. Both $\mathbf{W}$ and Ω being skew-symmetric, it results that: $\mathbf{W}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)=0$ and the same in the case of considering Ω, implying $\mathbf{G}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)=\frac{\xi }{\xi +\mu }\mathbf{D}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)$, and then:

$${\dot{\mathbf{p}}}_i=\mathbf{G}·{\mathbf{p}}_i-\frac{\xi }{\xi +\mu }\left(\mathbf{D}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)\right){\mathbf{p}}_i$$

(24)

which can be rewritten as:

$${\dot{\mathbf{p}}}_i=\frac{\mu }{\xi +\mu }\mathbf{W}·{\mathbf{p}}_i+\frac{\xi }{\xi +\mu }\nabla \mathbf{v}·{\mathbf{p}}_i-\frac{\xi }{\xi +\mu }\left(\mathbf{D}:({\mathbf{p}}_i\otimes {\mathbf{p}}_i)\right){\mathbf{p}}_i=$$

$$=\frac{\xi }{\xi +\mu }{\dot{\mathbf{p}}}_i^J+\frac{\mu }{\xi +\mu }\mathbf{W}·{\mathbf{p}}_i=\frac{\xi }{\xi +\mu }{\dot{\mathbf{p}}}_i^J+\frac{\mu }{\xi +\mu }{\dot{\mathbf{p}}}_i^R$$

(25)

where ${\dot{\mathbf{p}}}_i^J$ represents the hydrodynamic contribution in the absence of rod-rod interactions (dilute regime described by the Jeffery’s equation) and ${\dot{\mathbf{p}}}_i^R$, the one coming from the rods’ entanglements, which results in a rigid-like cluster kinematics.

We can notice that when $\xi \gg \mu$, hydrodynamic effects are preponderant, and the rod kinematics is governed by Jeffery’s equation, i.e., ${\dot{\mathbf{p}}}_i\approx {\dot{\mathbf{p}}}_i^J$. In the opposite case, $\mu \gg \xi$, the cluster is too rigid, and the rods adopt the velocity dictated by the rigid cluster kinematics, ${\dot{\mathbf{p}}}_i\approx {\dot{\mathbf{p}}}_i^R$.

If we assume the orientation of the rods of a deformable cluster given by the orientation distribution, $\psi \left(\mathbf{p}\right)$, the time derivative of its second order moment, $\dot{\mathbf{a}}$ results in:

$$\dot{\mathbf{a}}={\int }_{\mathcal{S}(0,1)}\left(\dot{\mathbf{p}}\otimes \mathbf{p}+\mathbf{p}\otimes \dot{\mathbf{p}}\right)\psi \left(\mathbf{p}\right)d\mathbf{p}$$

(26)

By considering the expression of the microscopic velocity, $\dot{\mathbf{p}}$:

$$\dot{\mathbf{p}}=\frac{\xi }{\xi +\mu }{\dot{\mathbf{p}}}^J+\frac{\mu }{\xi +\mu }{\dot{\mathbf{p}}}^R$$

(27)

it results that:

$$\dot{\mathbf{a}}=\frac{\xi }{\xi +\mu }{\dot{\mathbf{a}}}^J+\frac{\mu }{\xi +\mu }{\dot{\mathbf{a}}}^R$$

(28)

with ${\dot{\mathbf{a}}}^J$ and ${\dot{\mathbf{a}}}^R$ being the resulting expressions when Equation (27) is introduced into Equation (26) by assuming $\mu =0$ and $\xi =0$, respectively:

$$\left\{\begin{array}{c}{\dot{\mathbf{a}}}^J=\Omega ·\mathbf{a}-\mathbf{a}·\Omega +\mathbf{D}·\mathbf{a}+\mathbf{a}·\mathbf{D}-2\mathbf{A}:\mathbf{D}\hfill \\ {\dot{\mathbf{a}}}^R=\mathbf{W}·\mathbf{a}+\mathbf{a}·{\mathbf{W}}^T\hfill \end{array}\right.$$

(29)

with $\mathbf{A}$ being the fourth order moment defined by:

$$\mathbf{A}={\int }_{\mathcal{S}(0,1)}\mathbf{p}\otimes \mathbf{p}\otimes \mathbf{p}\otimes \mathbf{p}\psi \left(\mathbf{p}\right)d\mathbf{p}$$

(30)

The objectivity of the resulting evolution equation of $\dot{\mathbf{a}}$ was proven in [28]. Usually, the fourth order moment, $\mathbf{A}$, is expressed from the second order one by considering any of the numerous closure relations proposed in the literature [31,32,33,34]. We do not address this issue here, and from now on, we consider that $\mathbf{A}=\mathcal{L}\left(\mathbf{a}\right)$, with $\mathcal{L}$, an appropriate algebraic operator.

Until now, the deformation mechanisms considered were purely viscous; however, it seems natural that the cluster can experience a elastic deformation with a fading-memory, that is, the reference configuration is a combination of the recent experienced deformations.

3.1. Fading Elasticity

Clusters immersed in a flowing suspension experience two types of transformations, a rigid-like one associated with rigid translations and rotations and a second one implying the cluster stretching. When the cluster conformation is described from the second order moment, $\mathbf{a}$, of its orientation distribution, ψ, the simplest way of quantifying stretching is by evaluating its eigenvalues change. In what follows, we consider that clusters try to recover a configuration that combines the recent ones by a sort of fading-memory, but only in what concerns the cluster stretching, that is, the memory does not concern the rigid movements.

We assume the orientation tensor, ${\mathbf{a}}_{\tau }$, known at times, $\tau \le t$. Tensor ${\mathbf{a}}_{\tau }$ can be diagonalized according to:

$${\mathbf{a}}_{\tau }={\mathbf{R}}_{\tau }^T·{\Lambda }_{\tau }·{\mathbf{R}}_{\tau }$$

(31)

The reference configuration, ${\Lambda }_t^r$, that the cluster tries to recover at time, t, reads:

$${\Lambda }_t^r=\underset{0}{\overset{t}{\int }}m(t-\tau ){\Lambda }_{\tau }d\tau$$

(32)

where the memory function, $m(t-\tau )$, decreases as $t-\tau$ increases, with $m(t-\tau )=0$ for $t-\tau \ge \theta$ and:

$$\underset{0}{\overset{t}{\int }}m(t-\tau )d\tau =1$$

(33)

As soon as the reference configuration at time t, ${\Lambda }_t^r$, is available, the reference second order moment at time t, ${\mathbf{a}}_t^r$, can be obtained by applying:

$${\mathbf{a}}_t^r={\mathbf{R}}_t^T·{\Lambda }_t^r·{\mathbf{R}}_t$$

(34)

If the characteristic time of the elastic recover is $\mathcal{T}$, then the equation governing the evolution of the second order tensor reads:

$$\dot{\mathbf{a}}=\frac{\xi }{\xi +\mu }{\dot{\mathbf{a}}}^J+\frac{\mu }{\xi +\mu }{\dot{\mathbf{a}}}^R+\frac{{\mathbf{a}}^r-\mathbf{a}}{\mathcal{T}}$$

(35)

Memory effects disappear when θ vanishes. In this case, the memory function approaches the Dirac’s delta distribution, $m(t-\tau )=\delta (t-\tau )$. Thus, the reference configuration corresponds to the present one, i.e., ${\mathbf{a}}_t^r={\mathbf{a}}_t$, implying the cancellation of the last term in Equation (28). Elasticity effects are also associated to Brownian effects, as discussed in [28,35,36], as well as to the rods’ bending [37,38].

4. Kinetic Theory Description

If we consider a microstructure consisting of a population of clusters and the kinematics of these clusters only depends on the second order moment of the orientation distribution of the rods that constitute them, the microstructure could be fully described from the clusters’ conformation distribution function, $\Psi (\mathbf{x},t,\mathbf{a})$, given the fraction of clusters that at position $\mathbf{x}$ and time t have a configuration described by $\mathbf{a}$, the last related to the orientation of the rods composing it, which was noted by $\psi \left(\mathbf{p}\right)$.

This approach constitutes a two scale kinetic theory description, with the coarsest scale described by Ψ, depending on tensor $\mathbf{a}$ and tensor $\mathbf{a}$ at the finest scale depending on the rods’ orientation distribution $\psi \left(\mathbf{p}\right)$.

Function $\Psi \left(\mathbf{a}\right)$ is subjected locally (at each position, $\mathbf{x}$, and time t) to the normality condition:

$${\int }_{\mathcal{A}}\Psi \left(\mathbf{a}\right)d\mathbf{a}=1$$

(36)

The balance related to Ψ leads to the Fokker-Planck equation:

$$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\mathbf{v}\Psi \right)+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)=0$$

(37)

where in the 2D case, ${\nabla }_x^T=\left(\frac{\partial }{\partial x_1},\frac{\partial }{\partial x_2}\right)$ and ${\nabla }_a^T=\left(\frac{\partial }{\partial a_{11}},\frac{\partial }{\partial a_{12}}\right)$, and where $\dot{\mathbf{a}}$ is given by Equations (28) and (29):

$$\dot{\mathbf{a}}=\frac{\xi }{\xi +\mu }\left(\Omega ·\mathbf{a}-\mathbf{a}·\Omega +\mathbf{D}·\mathbf{a}+\mathbf{a}·\mathbf{D}-2\mathbf{A}:\mathbf{D}\right)+\frac{\mu }{\xi +\mu }\left(\mathbf{W}·\mathbf{a}-\mathbf{a}·\mathbf{W}\right)+\frac{{\mathbf{a}}^r-\mathbf{a}}{\mathcal{T}}$$

(38)

with $\mathbf{W}$ given by Equation (15): $\mathbf{W}=\Omega +\mathbf{D}·\mathbf{a}-\mathbf{a}·\mathbf{D}$; $\mathbf{A}=\mathcal{L}\left(\mathbf{a}\right)$ and ${\mathbf{a}}^r$, defined in Equation (34).

5. Introducing Interaction Mechanisms

In concentrated regimes, clusters interact, modifying their configuration. We can distinguish two limit cases, the ones associated with soft and rigid clusters. Here, we do not consider the aggregation and breaking mechanisms that were addressed in [25,28].

5.1. Soft Clusters

In the case of soft clusters, we can expect that the orientation configuration tends to reach an isotropic distribution, i.e., $\psi \left(\mathbf{p}\right)=1/2\pi$ in 2D. The second order moment related to such an isotropic distribution results in ${\mathbf{a}}^{iso}=\frac{1}{2}\mathbf{I}$, where $\mathbf{I}$ denotes the unit tensor, ${\mathbf{I}}_{ij}={\delta }_{ij}$, with δ, the Kronecker’s delta. This fact implies that the clusters’ conformation distribution should evolve towards a Dirac distribution, i.e., ${\Psi }^{iso}\left(\mathbf{a}\right)=\delta (\mathbf{a}-{\mathbf{a}}^{iso})$.

This mechanism could be modeled by introducing a kind of collision term, ${\dot{\Psi }}^C$, in the right-hand member of the Fokker-Planck equation (37):

$${\dot{\Psi }}^C=\frac{\xi }{\xi +\mu }\frac{{\Psi }^{iso}-\Psi }{\Theta (\dot{\gamma },\varphi )}$$

(39)

with Θ, a characteristic time depending on the flow intensity, $\dot{\gamma }$, the last described by the second invariant of the strain rate tensor, $\mathbf{D}$, and, also, on the clusters’ concentration, $\varphi$. In the previous expression, the pre-factor, $\xi /(\xi +\mu )$, quantifies the cluster softness.

5.2. Rigid Clusters

In the case of rigid clusters, interactions can modify the cluster orientation, but not its intrinsic orientation distribution. Thus, the second order orientation tensor can rotate, but its stretching is prevented by the cluster rigidity.

Thus, the orientation tensor can change, but their eigenvalues must remain unchanged. It is easy to verify that the eigenproblem involves the characteristic equation:

$${\lambda }^2-\lambda \text{tr}\left(\mathbf{a}\right)+\text{det}\left(\mathbf{a}\right)=0$$

(40)

where taking into account the unit trace of tensor $\mathbf{a}$, $\text{tr}\left(\mathbf{a}\right)=a_{11}+a_{22}=1$, proves that all the orientation tensors having the same determinant have the same eigenvalues.

The 2D conformational domain is defined by $(a_{11},a_{12})\in {\Omega }^c=\mathcal{B}(\mathbf{C},R)$, which represents the ball centered at point $\mathbf{C}=(0.5,0)$, having a radius, $R=0.5$. By developing the determinant, $\text{det}\left(\mathbf{a}\right)=a_{11}·a_{22}-a_{12}^2=a_{11}-a_{11}^2-a_{12}^2$, it can be noticed that the determinant is constant along the circumferences centered at point $\mathbf{C}$. Thus, diffusion takes place along the circumferential direction defined by $\mathbf{t}$, $\mathbf{t}$ being the unit vector defining the tangent direction to the circles centered at $\mathbf{C}$:

$${\mathbf{t}}^T=\frac{(-a_{12},a_{11}-0.5)}{\sqrt{a_{12}^2+{(a_{11}-0.5)}^2}}$$

(41)

The diffusive flux, ${\mathbf{q}}_a$, reads:

$${\mathbf{q}}_a=-D_a\left({\nabla }_a\Psi ·\mathbf{t}\right)\mathbf{t}=-{\mathbf{D}}_a·{\nabla }_a\Psi$$

(42)

with the diffusivity tensor, ${\mathbf{D}}_a$, given by:

$${\mathbf{D}}_a=D_a\left(\mathbf{t}\otimes \mathbf{t}\right)$$

(43)

If we assume that the diffusion coefficient, $D_a$, scales with the cluster aspect ratio, according to $D_a={\mathcal{D}}_a\left(a_{12}^2+{(a_{11}-0.5)}^2\right)$, the diffusivity tensor reads:

$${\mathbf{D}}_a={\mathcal{D}}_a\left(\begin{array}{cc}{a}_{12}^2& -a_{12}\left(a_{11}-0.5\right)\\ -a_{12}\left(a_{11}-0.5\right)& {\left(a_{11}-0.5\right)}^2\end{array}\right)$$

(44)

where ${\mathcal{D}}_a$ will depend on the flow intensity, $\dot{\gamma }$, and the clusters’ concentration, $\varphi$.

The diffusive flux given by ${\mathbf{q}}_a=-{\mathbf{D}}_a·{\nabla }_a\Psi$ writes:

$${\mathbf{q}}_a=-{\mathcal{D}}_a\left(\begin{array}{cc}{a}_{12}^2& -a_{12}\left(a_{11}-0.5\right)\\ -a_{12}\left(a_{11}-0.5\right)& {\left(a_{11}-0.5\right)}^2\end{array}\right)·\left(\begin{array}{c}\frac{\partial \Psi }{\partial a_{11}}\hfill \\ \frac{\partial \Psi }{\partial a_{12}}\hfill \end{array}\right)$$

(45)

This diffusion mechanism is introduced into the right-hand member of the Fokker-Planck Equation (37) from the contribution, ${\dot{\Psi }}^{D_a}$ :

$${\dot{\Psi }}^{D_a}=-\frac{\mu }{\xi +\mu }{\nabla }_a·{\mathbf{q}}_a$$

(46)

where the pre-factor quantifies the cluster rigidity.

In the case of rigid clusters, it is easy to prove that both, advection and diffusion, operate in the circumferential direction, that is, along the circles centered at point $\mathbf{C}$, having a radius, $r\le R=0.5$. Thus, one could think that a better representation could consist of considering a more appropriate choice for the conformational coordinates, allowing a simple expression of diffusion effects. One possibility lies in the choice (in the 2D case) of λ and θ, the first being the highest eigenvalue of tensor $\mathbf{a}$, which takes values in $\lambda =[0.5,1]$, and the second one, the orientation of the principal direction related to the highest eigenvalue, $\theta \in [0,2\pi ]$. Even if such a representation seems better for representing the diffusion terms, it is not the case for describing the conformational velocities in the general case (deformable aggregates).

5.3. Diffusion in the Physical Space

Because of the clusters’ interactions, we can also expect the diffusion of the cluster distribution, Ψ, in the physical space. This flux will be noted by ${\mathbf{q}}_x$ and will be modeled by an isotropic diffusion coefficient, ${\mathcal{D}}_x$, that could depend on the flow intensity, $\dot{\gamma }$, as well as on the clusters’ concentration, ϕ:

$${\mathbf{q}}_x=-{\mathcal{D}}_x{\nabla }_x\Psi$$

(47)

In the non-isotropic case, it suffices to replace ${\mathcal{D}}_x$ by the diffusion tensor, ${\mathbf{D}}_x$, operating in the physical space. In what follows, we are considering an isotropic diffusion in the physical space.

The diffusion mechanism is introduced into the right-hand member of the Fokker-Planck equa-tion (37) from the contribution, ${\dot{\Psi }}^{D_x}$ :

$${\dot{\Psi }}^{D_x}=-{\nabla }_x·{\mathbf{q}}_x$$

(48)

5.4. Extended Fokker-Planck Equation

When all these mechanisms apply, the Fokker-Planck equation (37) becomes:

$$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\mathbf{v}\Psi \right)+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)=$$

$$={\nabla }_x·\left({\mathcal{D}}_x{\nabla }_x\Psi \right)+\frac{\mu }{\xi +\mu }{\nabla }_a·\left({\mathbf{D}}_a·{\nabla }_a\Psi \right)+\frac{\xi }{\xi +\mu }\frac{{\Psi }^{iso}-\Psi }{\Theta }$$

(49)

6. Numerical Methods

The main issue related to kinetic theory descriptions lies in its multidimensionality, because generalized Fokker-Planck equations involve the physical coordinates (space and time) and a number of conformational coordinates describing the microstructure at the desired level of detail. Classically, the curse of dimensionality was circumvented by different ways, the most usual being:

• the solution of the stochastic counterpart. However, it is well known that the control of the statistical noise and variance reduction are the main tricky issues, making possible the calculation of the distribution moments, but making difficult the reconstruction of the distribution itself [4];
• the derivation of macroscopic approaches, sometimes purely phenomenological, other times inspired or derived from finer descriptions, but even in the last case, the use of closure approximations is generally compulsory;
• introducing a first discretization by substituting conformational coordinates by a series of populations. Thus, for example, starting from:

  $$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\mathbf{v}\Psi \right)+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)=0$$

  (50)

  we could introduce different populations related to particular values of ${\mathbf{a}}_i$: $a_{{11}_i}$ and $a_{{12}_i}$,$i=1,\cdots ,\mathcal{N}$ and describe the microstructure of each one from its distribution function, ${\Psi }_i(\mathbf{x},t)=\Psi (\mathbf{x},t;{\mathbf{a}}_i)$, whose evolution, after discretizing with respect to the coordinate that was substituted by the populations, $\mathbf{a}$, in the present case, is given by the coupled lower dimensional Fokker-Planck equations:

  $$\frac{\partial {\Psi }_i}{\partial t}+{\nabla }_x·\left(\mathbf{v}{\Psi }_i\right)+\sum _{j=1}^{\mathcal{N}}{c}_j^i{\Psi }_j=0;i=1,\cdots ,\mathcal{N}$$

  (51)

  This approach was considered in the kinetic description of systems of active particles in order to circumvent the curse of dimensionality [39]; however, as soon as the dimensionality increases, the number of populations explodes.
• the use of separated representations in order to ensure that the complexity scales linearly with the model dimensionality. This technique, introduced some years ago in [40] for steady simulations and in [41] for transient solutions, consists of expressing the unknown field as a finite sum of functional products, i.e., expressing a generic multidimensional function, $u(x_1,\cdots ,x_d)$, as:

  $$u(x_1,\cdots ,x_d)\approx \sum _{i=1}^{i=N}{F}_i^1\left(x_1\right)\cdots F_i^d\left(x_d\right)$$

  (52)

  The interested reader can refer to [42,43,44,45,46,47,48,49,50] and the references therein for a deep analysis of this technique and its applications in computational rheology.

When the partial differential equation is purely advective, the method of particles constitutes an appealing technique for solving it. As soon as diffusion terms become dominant, separated representations can be efficiently applied. At present, the main issue related to the use of separated representations techniques for solving purely advective equations is the necessity of considering accurate stabilizations [51].

6.1. Method of Particles for Solving Advection-Dominanted Problems

This technique, deeply described in [52,53], consists in approximating the initial distribution, $\Psi (\mathbf{x},t=0,\mathbf{a})$ from $\mathcal{M}$, Dirac’s masses, ${\mathbf{a}}_i^0$, at each one of the $\mathcal{Q}$ positions, ${\mathbf{x}}_j^0$:

$$\Psi (\mathbf{x},t=0,\mathbf{a})=\sum _{j=1}^{\mathcal{Q}}\sum _{i=1}^{\mathcal{M}}{\alpha }_i^j\delta (\mathbf{a}-{\mathbf{a}}_i^0)\delta (\mathbf{x}-{\mathbf{x}}_j^0)$$

(53)

which represents a sort of approximation based on $\mathcal{Q}\times \mathcal{M}$ pseudo-particles, ${\mathcal{P}}_{ij}$, with initial positions and conformations given by:

$$\left\{\begin{array}{c}{\mathbf{x}}_{ij}^0={\mathbf{x}}_j^0,i=1,\cdots ,\mathcal{M};j=1,\cdots ,\mathcal{Q}\hfill \\ {\mathbf{a}}_{ij}^0={\mathbf{a}}_i^0,i=1,\cdots ,\mathcal{M};j=1,\cdots ,\mathcal{Q}\hfill \end{array}\right.$$

(54)

and whose position and conformation will be evaluated all along the flow simulation, from which the distribution will be reconstructed.

6.1.1. Advection Equation

When considering the purely advective balance equation:

$$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\mathbf{v}\Psi \right)+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)=0$$

(55)

the time evolution of position and conformation of each particle, ${\mathcal{P}}_{ij}$, is calculated by integrating:

$$\left\{\begin{array}{c}{\mathbf{x}}_{ij}\left(t\right)={\mathbf{x}}_{ij}^0+{\int }_{\tau =0}^{\tau =t}\mathbf{v}\left({\mathbf{x}}_{ij}\left(\tau \right)\right)d\tau \hfill \\ {\mathbf{a}}_{ij}\left(t\right)={\mathbf{a}}_{ij}^0+{\int }_{\tau =0}^{\tau =t}{\dot{\mathbf{a}}}_{ij}({\mathbf{a}}_{ij}\left(\tau \right),{\mathbf{x}}_{ij}\left(\tau \right))d\tau \hfill \end{array}\right.$$

(56)

As the position update only depends on the velocity field, which itself only depends on the position, it can be stressed that particles, ${\mathcal{P}}_{ij}$, $i=1,\cdots ,\mathcal{M}$, are following the same trajectory in the physical space, having as a departure point, ${\mathbf{x}}_j^0$.

Now, the orientation distribution at time t can be reconstructed from:

$$\Psi (\mathbf{x},t,\mathbf{a})=\sum _{j=1}^{\mathcal{Q}}\sum _{i=1}^{\mathcal{M}}{\alpha }_i^j\delta (\mathbf{a}-{\mathbf{a}}_{ij}\left(t\right))\delta (\mathbf{x}-{\mathbf{x}}_{ij}\left(t\right))$$

(57)

Obviously, smoother representations can be obtained by considering appropriate regularizations of the Dirac’s distribution, as the one usually performed within the SPH (Smooth Particles Hydrodynamics) framework [52,54].

When balance equations involve diffusion terms, there are two main routes based on the use of particles methods, one of a stochastic nature, the other fully deterministic.

6.1.2. Advection-Diffusion Equation

For illustrating the procedure when models involve diffusion terms, we consider the Fokker-Planck equation:

$$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\mathbf{v}\Psi \right)+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)=-{\nabla }_x·{\mathbf{q}}_x-{\nabla }_a·{\mathbf{q}}_a$$

(58)

where ${\mathbf{q}}_x$ and ${\mathbf{q}}_a$ are two diffusive fluxes operating in the physical and conformational spaces, respectively, both modeled from a Fick’s type law:

$$\left\{\begin{array}{c}{\mathbf{q}}_x=-{\mathcal{D}}_x{\nabla }_x\Psi \hfill \\ {\mathbf{q}}_a=-{\mathbf{D}}_a·{\nabla }_a\Psi \hfill \end{array}\right.$$

(59)

• Within the stochastic framework diffusion, terms can be modeled from appropriate random variables within a Lagrangian or a Eulerian description, the last one known as Brownian Configurations Fields (BCF). Both approaches were considered in our former works on the solution of Fokker-Planck equations [53,55]. Within the Lagrangian stochastic framework and starting from the initial cloud of pseudo-particles, ${\mathcal{P}}_{ij}$, representing the initial cluster distribution, $\Psi (\mathbf{x},t=0,\mathbf{a})$, the simplest particles updating writes:

  $$\left\{\begin{array}{c}{\mathbf{x}}_{ij}\left(t_{n+1}\right)={\mathbf{x}}_{ij}\left(t_n\right)+\mathbf{v}\left({\mathbf{x}}_{ij}\left(t_n\right)\right)\Delta t+{\mathcal{R}}_x(\Delta t)\hfill \\ {\mathbf{a}}_{ij}\left(t_{n+1}\right)={\mathbf{a}}_{ij}\left(t_n\right)+{\dot{\mathbf{a}}}_{ij}({\mathbf{a}}_{ij}\left(t_n\right),{\mathbf{x}}_{ij}\left(t_n\right))\Delta t+{\mathcal{R}}_a(\Delta t)\hfill \end{array}\right.$$

  (60)

  where $\Delta t$ is the time step and both random updates, ${\mathcal{R}}_x$ and ${\mathcal{R}}_a$, depend on the chosen time step (see [56] for more details and for advanced stochastic integrations).
  Obviously, because the random effects operate in the physical space, the $\mathcal{M}$ particles initially located at each position, ${\mathbf{x}}_j^0$, $j=1,\cdots ,\mathcal{Q}$, will follow different trajectories in the physical space along the simulation time. In order to obtain accurate enough results, we must consider a rich enough representation, that is, a large population of particles. For this purpose, we must consider large enough $\mathcal{M}$ and $\mathcal{Q}$. The large number of particles to be tracked seems a disadvantage of the approach at first view, but it must be noticed that the integration of each particle is completely independent of all the others, making possible the use of HPCon massive parallel computing platforms.
• The technique introduced for treating purely advective equations can be extended for considering diffusion contributions, as was described in [54], and which we revisit in what follows, within a fully deterministic approach.
  Equation (58) can be rewritten as:

  $$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\left(\mathbf{v}+\frac{{\mathbf{q}}_x}{\Psi }\right)\Psi \right)+{\nabla }_a·\left(\left((\dot{\mathbf{a}}+\frac{{\mathbf{q}}_a}{\Psi }\right)\Psi \right)=0$$

  (61)

  or

  $$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\tilde{\mathbf{v}}\Psi \right)+{\nabla }_a·\left(\dot{\tilde{\mathbf{a}}}\Psi \right)=0$$

  (62)

  where the effective velocities, $\tilde{\mathbf{v}}$ and $\dot{\tilde{\mathbf{a}}}$, are given by:

  $$\left\{\begin{array}{c}\tilde{\mathbf{v}}=\mathbf{v}-\frac{{\mathcal{D}}_x}{\Psi }{\nabla }_x\Psi \hfill \\ \dot{\tilde{\mathbf{a}}}=\dot{\mathbf{a}}-\frac{1}{\Psi }{\mathbf{D}}_a·{\nabla }_a\Psi \hfill \end{array}\right.$$

  (63)

  Now, the integration schema (56) can be applied by replacing material and conformation velocities, $\mathbf{v}$ and $\dot{\mathbf{a}}$, by their effective counterparts, $\tilde{\mathbf{v}}$ and $\dot{\tilde{\mathbf{a}}}$:

  $$\left\{\begin{array}{c}{\mathbf{x}}_{ij}\left(t\right)={\mathbf{x}}_{ij}^0+{\int }_{\tau =0}^{\tau =t}\tilde{\mathbf{v}}\left({\mathbf{x}}_{ij}\left(\tau \right)\right)d\tau \hfill \\ {\mathbf{a}}_{ij}\left(t\right)={\mathbf{a}}_{ij}^0+{\int }_{\tau =0}^{\tau =t}{\dot{\tilde{\mathbf{a}}}}_{ij}({\mathbf{a}}_{ij}\left(\tau \right),{\mathbf{x}}_{ij}\left(\tau \right))d\tau \hfill \end{array}\right.$$

  (64)

  This fully deterministic particle description requires much less particles than its stochastic counterpart, but as noticed in Equation (63), the calculation of the effective material and conformational velocities require the derivative of the cluster distribution, Ψ, with respect to both the physical and the conformational coordinates, and for that purpose, this distribution must be reconstructed all along the simulation (at each time step), a fact that constitutes a serious drawback for its implementation in massive parallel computing platforms. Moreover, to make possible the calculation of the distribution derivatives, the Dirac’s distribution must be regularized in order to ensure its derivability.

6.1.3. Discussion

It is important to note:

• Many times, the initial configuration reduces to a very small subdomain within ${\Omega }^c$. For example, when the initial configuration (associated with the suspension at rest) corresponds to isotropic clusters, the cluster distribution reads, $\Psi (\mathbf{x},t=0,\mathbf{a})=\delta (\mathbf{a}-0.5\mathbf{I})$. In other cases, the initial configuration can be expressed from some configurations without needing an exhaustive approximation in the whole domain, ${\Omega }^c$;
• In the case of homogeneous flows, the solution is the same everywhere in the physical domain, and the distribution becomes independent of $\mathbf{x}$, that is, equivalent to considering $\mathcal{Q}=1$ in the previous analysis;
• In the case of homogeneous flows, in absence of diffusion effects and when the initial configuration reduced to a single point, ${\mathbf{a}}^0$, in the conformational domain, ${\Omega }^c$ (e.g., isotropic initial configuration), a single particle suffices for tracking the microstructure evolution, i.e., $\mathcal{Q}=1$ and $\mathcal{M}=1$;
• In the case of homogeneous flows, but in the presence of diffusion effects, when the initial condition can be expressed from a single point, ${\mathbf{a}}^0$, in the conformational domain, ${\Omega }^c$, stochastic simulations need considering, $\mathcal{Q}=1$, but $\mathcal{M}\gg 1$. Random contributions will produce different trajectories, ${\mathbf{a}}_i\left(t\right)$, all of them starting at the same initial configuration, ${\mathbf{a}}^0$;
• In the situation described in the previous item, but when considering a deterministic particles strategy, the initial location of the $\mathcal{M}$ particles must be perturbed slightly for locating all of them in a close neighborhood of ${\mathbf{a}}^0$, because when all of them are located exactly at configuration ${\mathbf{a}}^0$, the gradient of the smoothed distribution vanishes, and then, all the trajectories will coincide;
• If the initial configuration explores the whole domain, ${\Omega }^c$ (that in the 3D case involves a 5D configurational domain), we are faced with the redoubtable curse of dimensionality;
• The number of particles to be considered depends on the richness of the distribution to be represented in the physical and conformational spaces throughout the simulation.

6.2. Separated Representations for Solving Diffusion-Dominated Problems

When the diffusion effects are dominant, techniques presented in the previous section become inefficient, because they require an excessive number of particles to produce accurate enough results, in particular, for reconstructing the distribution. In this case, standard mesh-based discretizations seem a better choice. However, as discussed before, mesh-based discretizations fail when addressing highly dimensional models, as is the case when addressing the solution of the previous introduced Fokker-Planck equations. In that case, separated representations seem the most appealing choice.

If we consider again the Fokker-Planck equation:

$$\frac{\partial \Psi }{\partial t}+{\nabla }_x·\left(\mathbf{v}\Psi \right)+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)={\nabla }_x\left({\mathcal{D}}_x{\nabla }_x\Psi \right)+{\nabla }_a\left({\mathbf{D}}_a·{\nabla }_a\Psi \right)$$

(65)

There are many separated representation choices. The most natural consists of separating time, physical and conformational spaces, i.e.:

$$\Psi (\mathbf{x},t,\mathbf{a})\approx \sum _{i=1}^N{X}_i\left(\mathbf{x}\right)·T_i\left(t\right)·A_i\left(\mathbf{a}\right)$$

(66)

Thus, when proceeding with the Proper Generalized Decomposition—PGD—constructor [28], we must solve on the order of N 2D or 3D (depending on the dimension of the physical space) boundary value problems—BVP—for calculating functions, $X_i\left(\mathbf{x}\right)$, the same number of 1D initial value problems—IVP—for calculating functions, $T_i\left(t\right)$, and, finally, the same number of 2D or 5D (the number of components of tensor $\mathbf{a}$ in the 2D or 3D case, where the symmetry and its unit trace have been taken into account) BVP for calculating functions, $A_i\left(\mathbf{a}\right)$.

This choice has a major difficulty: the problems related to the computation of $A_i\left(\mathbf{a}\right)$ are defined in 5D in the general three-dimensional case. To circumvent this difficulty, one could imagine a fully separated representation of the configuration space according to:

$$\left\{\begin{array}{cc}{A}_i\left(\mathbf{a}\right)=A_i^{11}\left(a_{11}\right)·A_i^{12}\left(a_{12}\right)\hfill & \text{in}\text{2D}\hfill \\ A_i\left(\mathbf{a}\right)=A_i^{11}\left(a_{11}\right)·A_i^{22}\left(a_{22}\right)·A_i^{12}\left(a_{12}\right)·A_i^{13}\left(a_{13}\right)·A_i^{23}\left(a_{23}\right)\hfill & \text{in}\text{3D}\hfill \end{array}\right.$$

(67)

In what follows, for the sake of simplicity, we focus on the 2D case in which $A_i\left(\mathbf{a}\right)=A_i^{11}\left(a_{11}\right)·A_i^{12}\left(a_{12}\right)$. If we consider such a separated representation, the conformational domain results in $(a_{11},a_{12})\in {\tilde{\Omega }}^c=[0,1]\times [-0.5,0.5]$, even if, as reported before, the real conformation domain is the ball, ${\Omega }^c=\mathcal{B}(\mathbf{C},R)$, with center, $\mathbf{C}=(0.5,0)$, and radius, $R=0.5$. When proceeding in the extended domain, ${\tilde{\Omega }}^c$, due to the consistent physics considered, the solution vanishes in ${\tilde{\Omega }}^c-{\Omega }^c$.

6.3. Hybrid Strategy

A hybrid schema can be defined by using grouped conformational coordinates, that is, $A_i\left(\mathbf{a}\right)$ in Equation (66) defined in ${\Omega }^c=\mathcal{B}(\mathbf{C},R)$. In 2D, problems associated with functions, $A_i\left(\mathbf{a}\right)$, are 2D, and then, they can be solved by using any standard discretization technique, like finite elements with suitable convective stabilization when required. In the 3D case, the problems related to the calculation of functions, $A_i\left(\mathbf{a}\right)$, are defined in 5D, and then, mesh-based discretizations fail. Full separated representations, such as the ones given in Equation (67), are the most appealing choice when diffusion dominates. However, when advection dominates, a valuable alternative consists of using a method of particles, such as the ones previously described. However, it is important noticing that methods of particles are not suitable for directly computing steady state solutions (they result from the long-time solutions). Thus, an alternative separated representation reads:

$$\Psi (\mathbf{x},t,\mathbf{a})\approx \sum _{i=1}^N{X}_i\left(\mathbf{x}\right)·A_i^t(\mathbf{a},t)$$

(68)

where the time is grouped with the conformational coordinates, and then, particle-based techniques can be used for computing the transient solutions, $A_i^t$, in the conformational domain.

When performing the separated representation indicated in Equation (68), the problems to be solved for calculating functions, $A_i^t$, do not depend on $\mathbf{x}$, which is equivalent to the solution of a problem in a homogenous flow. All the remarks addressed in that sense when considering the methods of particles (summarized in Section 6.1.3.) apply in the present case.

7. Validating the Micro-Mechanical Models

After having proposed the micro-mechanical model, as well as the kinetic theory description, and widely discussing its numerical treatment, it is time to validate the proposed micro-mechanical model, before employing it to describe complex microstructures involved in concentrated flowing suspensions.

First of all, we are verifying the pertinence of the micro-mechanical model from a theoretical point of view. For this purpose, we first consider two particular configurations of rigid clusters for which the exact solution is known:

• Rigid cluster composed of a single rod. In this case, the exact solution corresponds to the particularization of the Jeffery’s solution [1] for an ellipsoidal particle of infinite aspect ratio, which writes:

  $$\dot{\mathbf{p}}=\nabla \mathbf{v}·\mathbf{p}-\left(\nabla \mathbf{v}:(\mathbf{p}\otimes \mathbf{p})\right)\mathbf{p}$$

  (69)

  In the case of a single rod, $N=1$, whose orientation is defined by $\mathbf{p}$, Equation (8) reduces to:

  $$\dot{\mathbf{p}}=\mathbf{\omega }\times \mathbf{p}=\left(\mathbf{p}\times (\nabla \mathbf{v}·\mathbf{p})\right)\times \mathbf{p}$$

  (70)

  and using the vector triple product formula, $\mathbf{a}\times (\mathbf{b}\times \mathbf{c})=\mathbf{b}(\mathbf{a}·\mathbf{c})-\mathbf{c}(\mathbf{a}·\mathbf{b})$, it results that:

  $$\dot{\mathbf{p}}=\nabla \mathbf{v}·\mathbf{p}-\left(\nabla \mathbf{v}:(\mathbf{p}\otimes \mathbf{p})\right)\mathbf{p}$$

  (71)

  which corresponds with the Jeffery’s equation for rods.
• We are now proving that the kinematics of a rigid ellipsoidal cluster composed of rods can be described by the Jeffery solution associated with an ellipsoid.
  In [28], we proved that for describing non-infinite aspect ratio 2D particles, it suffices to consider a rigid system composed of two rods, aligning perpendicularly one with respect to the other, and having lengths, $2L_1$ and $2L_2$, $L_1\ge L_2$, such that the parameter, r, in the general Jeffery’s equation:

  $$\dot{\mathbf{p}}=\Omega ·\mathbf{p}+r\mathbf{D}·\mathbf{p}-r\left(\mathbf{D}:(\mathbf{p}\otimes \mathbf{p})\right)\mathbf{p}$$

  (72)

  which depends on the ellipsoid aspect ratio, ν (length-diameter ratio):

  $$r=\frac{{\nu }^2-1}{{\nu }^2+1}$$

  (73)

  resulting from $\nu =\frac{L_1}{L_2}$.
  We proved that the kinematics of a rigid cluster only depends on the second order orientation tensor, $\left(\mathbf{a}\right)$, associated with the orientation distribution of its rods, $\psi \left(\mathbf{p}\right)$. Moreover, because of the symmetry of $\mathbf{a}$, it can be diagonalized, its eigenvalues in the 2D case being ${\lambda }_1$ and ${\lambda }_2$ (${\lambda }_1\ge {\lambda }_2$). We can evaluate the angular velocity of the rigid cluster characterized by the distribution moment, $\mathbf{a}$, according to Equation (14). It is easy to prove that this velocity is exactly the one related to the Jeffery’s one for an ellipsoid of aspect ratio, $\nu ={\lambda }_1^2/{\lambda }_2^2$.
• In the case of a rigid cluster with an isotropic configuration, which if $\mathbf{a}=0.5\mathbf{I}$ (in 2D), it can be deduced from Equations (12) and (15) that as expected, the cluster rotates with the velocity dictated by the flow vorticity, Ω.

The previous verifications prove the consistency of the micromechanical model previously elaborated and whose objectivity was also proven in [28]; however its validity for describing real physical situations needs further analyses.

For validating the micromechanical model proposed in Section 3, we are considering a direct numerical simulation consisting of a population of 500 rigid rods with a given initial orientation distribution, which evolves because of the enforced flow, as well as the interaction between the rods constituting the cluster. A dense enough cluster is needed in order to ensure rigid behaviors that need a sufficient number of contact points. Using less than 500 rods does not allow for describing rigid enough clusters. The behavior does not depend on the rod length, because the velocity gradient is assumed constant at the cluster level. The contact law is assumed viscous, and its intensity is modified for ranging from deformable to rigid clusters. For additional details on the direct simulation at the microscopic scale (at the scale of the rods), the interested readers can refer to [8,26,27] and the references therein.

Figure 3 depicts the initial isotropic configuration, where the yellow points represent the contact points at which different rods interact. Then, we enforce the simple shear flow defined by the velocity field, $\mathbf{v}=(x_2,0)$, and adjust the intensity of the contact viscous forces applied in the microscopic direct calculation at the contact points, which allow simulating aggregates of different softness, ranging from rigid to deformable. In all cases, during the cluster configuration evolution, the second order orientation tensor is calculated by considering the star configuration depicted in Figure 1, which results in joining the cluster center of gravity to each one of the beads that it involves.

We compare the cluster configuration from its second order orientation tensor obtained from the direct simulation and the one resulting from the model described in Section 3.

Figure 3. Direct microscopic simulation: isotropic rods population and contact points at the initial time. A dense enough cluster ensuring many contacts between the different rods is needed in order to represent rigid enough clusters.

Figure 3. Direct microscopic simulation: isotropic rods population and contact points at the initial time. A dense enough cluster ensuring many contacts between the different rods is needed in order to represent rigid enough clusters.

Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 compare model predictions versus direct simulation calculations when starting from an isotropic initial configuration for different cluster softness (ranging from very rigid to very soft). As we consider model (28), there is a single tunable parameter: $\mu /\xi$, $\xi \ne 0$. Thus, we could consider $\xi =1$ and look for the value of μ that allows for fitting the evolution calculated from the direct simulation.

With a single tunable parameter, the comparisons allow one to conclude the excellent accuracy (qualitative and quantitative) of the proposed micromechanical model, concerning both the orientation intensity and the evolution period.

Figure 4. Model versus direct simulation: identified $\mu =60$.

Figure 4. Model versus direct simulation: identified $\mu =60$.

Figure 5. Model versus direct simulation: identified $\mu =13$.

Figure 5. Model versus direct simulation: identified $\mu =13$.

Figure 6. Model versus direct simulation: identified $\mu =2.45$.

Figure 6. Model versus direct simulation: identified $\mu =2.45$.

Figure 7. Model versus direct simulation: identified $\mu =1.1$.

Figure 7. Model versus direct simulation: identified $\mu =1.1$.

Figure 8. Model versus direct simulation: identified $\mu =0.6$.

Figure 8. Model versus direct simulation: identified $\mu =0.6$.

8. Numerical Results

In this section, we are addressing preliminary results concerning the transient solution of the 2D Fokker-Planck equation (74) in a simple shear and homogeneous flow:

$$\frac{\partial \Psi }{\partial t}+{\nabla }_a·\left(\dot{\mathbf{a}}\Psi \right)=\frac{\mu }{\xi +\mu }{\nabla }_a·\left({\mathbf{D}}_a·{\nabla }_a\Psi \right)+\frac{\xi }{\xi +\mu }\frac{{\Psi }^{iso}-\Psi }{\Theta }$$

(74)

with $\Theta =1$, $\dot{\mathbf{a}}$ given by Equation (28) and:

$${\mathbf{D}}_a={\mathcal{D}}_a\left(\begin{array}{cc}{a}_{12}^2& -a_{12}\left(a_{11}-0.5\right)\\ -a_{12}\left(a_{11}-0.5\right)& {\left(a_{11}-0.5\right)}^2\end{array}\right)$$

(75)

for different values of ${\mathcal{D}}_a$ and $\mu /\xi$. The initial condition consists of an isotropic and normalized Gaussian distribution around the isotropic state.

Figure 9, Figure 10 and Figure 11 show the long-time cluster distributions for an effective diffusion coefficient, ${\mathcal{D}}_a=1$, and different values of the cluster softness, $\mu /\xi =100$, $\mu /\xi =1$ and $\mu /\xi =0.1$, respectively. The first corresponds almost to rigid clusters, whereas the last one corresponds to the softest ones. In the last case, we can notice that, as soon as clusters deform, they tend to align in the shear direction, $\theta =\pi /4$, and when the deformation progresses (larger radius), clusters tend to align in the flow direction, $\theta =0$. This mechanism is tempered by the diffusion that applies along the circumferential direction.

The previous solutions were calculated by applying a standard finite element discretization with a fine enough mesh of ${\Omega }^c$. When the problem is solved by considering a fully separated representation (67) in the extended conformational domain, ${\tilde{\Omega }}^c$, the solution related to rigid clusters, $\mu /\xi =100$, involved less than 10 modes, the number of terms involved in the solution associated with the softest aggregates, $\mu /\xi =0.1$, being around 40.

Figure 9. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mu /\xi =100$.

Figure 9. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mu /\xi =100$.

Figure 10. Long-time cluster distribution (left) and contour levels (right) $\Psi (t\to \infty ,\mathbf{a})$ for $\mu /\xi =1$.

Figure 10. Long-time cluster distribution (left) and contour levels (right) $\Psi (t\to \infty ,\mathbf{a})$ for $\mu /\xi =1$.

Figure 11. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mu /\xi =0.1$.

Figure 11. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mu /\xi =0.1$.

In order to point out the diffusion effects, we consider clusters with equal softness, $\mu /\xi =0.1$, but subjected to different values of the effective diffusion. Figure 12 and Figure 13 show the long-time cluster distributions for effective diffusion coefficients, ${\mathcal{D}}_a=0.1$ and ${\mathcal{D}}_a=100$, respectively. When the diffusion increases, we can notice that the flow deforms the clusters at the same time that diffusion homogenizes their orientations. For the largest diffusion, the distribution is almost uniform on the circles centered at $\mathbf{C}$ and, consequently, the deformation mechanism seems less effective.

Figure 12. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mathcal{D}=0.1$.

Figure 12. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mathcal{D}=0.1$.

Figure 13. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mathcal{D}=100$.

Figure 13. Long-time cluster distribution (left) and contour levels (right), $\Psi (t\to \infty ,\mathbf{a})$, for $\mathcal{D}=100$.

9. Conclusions

We proposed in this paper a two scale kinetic theory model of concentrated suspensions. Its main ingredients are:

• A micromechanical model describing the configuration evolution of deformable aggregates, whose kinematics only depends on the second order moment of the orientation distribution of their rods;
• The introduction of a cluster distribution that depends on the conformation of the aggregates involved in the suspension;
• The consideration of many mechanisms at both scales, the finest related to the aggregates (fading-memory elasticity) and the coarsest associated with the suspension itself, taking into account the clusters interactions.

Because the resulting kinetic theory model is defined in a highly dimensional space, the numerical strategies available for solving it were discussed deeply.

The consistency and pertinence of the micromechanical model was verified and validated by reproducing some well established results and by comparing its predictions with direct numerical simulations considering a population of rods and taking into account, explicitly, all the interactions between them.

First simulations of the kinetic theory model were finally addressed, proving the interest of the proposed approach for addressing a complex and recurrent modeling issue, the one concerning the description of rich, evolving microstructures encountered in flowing concentrated suspensions of rods.

Deeper analyses from both the modeling and the simulation point of views are in progress.