# One and Two-Fiber Orientation Kinetic Theories of Fiber Suspensions

Miroslav Grmela, Amine Ammar, Francisco Chinesta

► To cite this version:

Miroslav Grmela, Amine Ammar, Francisco Chinesta. One and Two-Fiber Orientation Kinetic Theories of Fiber Suspensions. Journal of Non-Newtonian Fluid Mechanics, 2013, 200, pp.17-33. 10.1016/j.jnnfm.2012.10.009. hal-01007165

HAL Id: hal-01007165

https://hal.science/hal-01007165

Submitted on 2 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# One and two-fiber orientation kinetic theories of fiber suspensions

Miroslav Grmela<sup>a</sup>, Amine Ammar<sup>b</sup>, Francisco Chinesta<sup>c</sup>

<sup>a</sup>École Polytechnique de Montréal, C.P. 6079, succ. Centre-ville, Montréal, Québec, Canada H3C 3A7

<sup>b</sup>Arts et Métiers ParisTech, 2 Boulevard du Ronceray, 49035 Angers cedex 01, France

<sup>c</sup>EADS Corporate Foundation International Chair, GEM Institute, École Centrale de Nantes, BP 92101, 44321 Nantes cedex 3, France

The morphology influencing rheological properties of suspensions of rigid spheres constitutes the flow induced collective ordering of the spheres characterized by two or more sphere distribution functions. When the rigid spheres are replaced by rigid fibers, the collective order in the position of the spheres is replaced by the flow induced orientation of the fibers that suffices to be characterized by one-fiber orientation distribution function. A flow induced collective ordering of fibers (both in position and orientation), that can only be characterized by two or more fiber distribution functions, can still however constitute an important part of the morphology. We show that two types of interaction among fibers, one being the Onsager-type topological interaction entering the free energy and the other the hydrodynamics interaction entering the dissipative part of the time evolution, give indeed rise to a collective order in the orientation influencing the rheology of fiber suspensions.

Keywords:

Fiber orientation

Kinetic theory

Hamiltonian dynamics

Generic

## 1. Introduction

When the suspended particles are rigid spheres, imposed flows induce in them a collective organization that then transforms the suspension into a rheologically complex fluid. The two-sphere distribution function is needed to characterize such collective order. When the suspended particles are fibers, it is the flow induced collective orientation of the fibers that takes a dominant role in determining the rheology. The fiber orientation can be well described by the one-fiber distribution function. It is however still possible that imposed flows induce, in addition to the overall orientation of the fibers, also a finer structure in both the position and the orientation of the fibers, that contributes then to the rheological complexity of fiber suspensions. Such fine structure can only be described by two or more fiber distribution functions. The objective of this paper is to investigate at least some aspects of the presence of the flow induced collective order in fiber suspensions and of its role in rheology. We identify two types of interaction among fibers that indeed give rise to the collective ordering of fibers influencing rheological properties. The first such interaction is the topological interaction entering the free energy. This type of interaction was introduced first by Onsager [1] in his equilibrium analysis of liquid crystals. The second such interaction arises in the consideration of fiber-fluid friction-type interactions. We derive a two-fiber kinetic equation that takes into account these two interactions and a corresponding to it formula expressing

the extra stress tensor in terms of the two-fiber distribution function. By solving numerically the kinetic equation we obtain a two-fiber correlation function describing the collective ordering of the fibers and rheological predictions in imposed simple shear flows.

The investigation reported in this article relates to previous investigations in Stokesian dynamics [2,3] and in kinetic theories [4–6] of semidilute fiber suspensions. The passage from the governing equations of Stokesian dynamics of *N* fibers to kinetic equations consists of, first, writing the Liouville equation that corresponds to the time evolution equations of *N* fibers, and second, reducing it to equations governing one or a small number fiber distribution functions. We shall not attempt in this article to make such passage. Koch and Shaqfeh in [5] investigate corrections (see Eq. (38) in [5]) to Jeffery's one-fiber kinetic equation (see [4]) that arise due to the presence of other fibers in the suspension. The isotropic medium in which the fiber whose orientation is followed in Jeffery kinetic theory is suspended is replaced with a linear anisotropic medium characterized by a unit vector **n**. In our investigation the vector **n** becomes the orientation vector of the second fiber followed in the two-fiber kinetic theory.

Our analysis presented in this article contributes to the study of semidilute fiber suspensions in the following four aspects:

(i) The arguments on the basis of which we derive one and two-fiber kinetic equations and the corresponding to them extra stress tensors are top-down. This means that our starting point is an overall structure of the governing equations guaranteeing certain important properties of their solutions. The top-down approach is complementary to the more
frequently used bottom-up approach that begins with a detailed microscopic analysis of fiber dynamics. We shall comment about both approaches in more detail in Section 2.2.

(ii) Both nondissipative (i.e. Jeffery's or modified Jeffery's type) and compatible with them dissipative parts are included in the one and two-fiber kinetic equations discussed in this article.

(iii) The Onsager topological interaction among fibers, originally introduced by Onsager [1] only in the equilibrium theory of fiber suspensions, is taken into account in kinetic equations.

(iv) Both morphological and rheological consequences of the one and two-fiber kinetic equations are investigated in this article.

The analysis remains incomplete in particular in the following three aspects.

First, we restrict our investigation to the order in the orientation and ignore the possible fine-structure ordering in the positions of the centers of mass of the fibers. We make this limitation for two reasons. The overall orientation of the fibers is known to play a dominant role in the rheology of fiber suspensions. We thus expect that also a possible fine structure in the orientation, that can only be seen in two or more fiber distribution functions, will play a more important role in rheology than a possible fine structure in the position of the fibers. We recall that the presence of migration (i.e. formation of shear induced spatial inhomogeneities in fiber suspensions), that may be related also to the formation of a fine spatial structure, has been observed and discussed in [7,6]. The second reason for the limitation to the orientation order is the simplicity. However, the general method with which we discuss below the two-fiber orientation kinetic theory can also be applied to the two-fiber position and orientation kinetic theory. The way to proceed to make such extension is indicated during the development of the two-fiber orientation kinetic theory.

Second, we formulate both the kinetic theory and the corresponding to it rheology in the top-down manner without complementing it with bottom-up considerations. However, some results obtained in [5] in the context of microhydrodynamic investigations of semidilute fiber suspensions do complement our top-down derivation. The advantage of the bottom-up approach is the clarity of the microscopic physics involved, its disadvantage lies in the difficulties to transform the microscopic physics into the macroscopic physics of suspensions. In the top-down approach the disadvantages and advantages are interchanged.

Third, we offer only rheological observations as an indirect experimental access to the orientation ordering of fibers. We do not investigate alternative, as for example optical, methods.

The paper is organized as follows. We begin (in Section 2.1) with definitions of the distribution functions serving as the state variables characterizing the morphology of fiber suspensions. Then, in Section 2.2, we introduce general features of the top-down method that we use throughout this paper to develop kinetic theories. In Section 2.3, we present the method by re-deriving with it Jeffery's one-fiber orientation kinetic theory [4]. In Section 2.4 we then extend Jeffery's theory to a two-fiber orientation kinetic theory. Section 3 is devoted to presenting the numerical method (the Proper Generalized Decomposition (PGD) method that is recalled in Appendix) used to solve the kinetic equations and to presenting the results.

## 2. Orientation kinetic theory

In this section we develop one and two-fiber orientation kinetic theory of fiber suspensions. We begin by defining the fiber

distribution functions and then continue by specifying the equations governing their time evolution and compatible with them expressions for the extra stress tensor.

### 2.1. State variables

States of a single fiber are characterized by

$$ (\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m}) \qquad (1) $$

where $(\mathbf{r}, \mathbf{v})$ is the position vector and momentum of the center of mass of the fiber, $\mathbf{p}$ is the unit vector along the fiber, and $\mathbf{m}$ its angular momentum. Similarly, states of two fibers are characterized by

$$ (\mathbf{r}_1, \mathbf{r}_2, \mathbf{v}_1, \mathbf{v}_2, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2) \qquad (2) $$

where the index "1" and "2" denote respectively the fiber one and the fiber two. Alternatively, we rewrite (2) in new coordinates

$$ (\mathbf{r}, \mathbf{v}, \mathbf{R}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2) \qquad (3) $$

where

$$ \begin{gathered} \mathbf{r} = \frac{1}{2}(\mathbf{r}_1 + \mathbf{r}_2); \quad \mathbf{R} = \mathbf{r}_2 - \mathbf{r}_1 \qquad (4) \\ \mathbf{v} = \mathbf{v}_1 + \mathbf{v}_2; \quad \mathbf{w} = \frac{1}{2}(\mathbf{v}_2 - \mathbf{v}_1) \end{gathered} $$

We recall that the angular velocity $\omega$ and the angular momentum $\mathbf{m}$ are related by $\omega = E_m$, where $E(\mathbf{p}, \mathbf{m})$ is the energy and $E_m$ is a shorthand notation (that we shall use throughout this paper) for $\partial E/\partial \mathbf{m}$. By involving explicitly the angular momentum (or equivalently the angular velocity) in the formulation of fiber dynamics we are taking into account the inertial of the fibers (see also the text following Eq. (12) below). The angular momentum will eventually disappear in the process of deriving the kinetic equations but we absolutely need it to start the process.

Now we consider all fibers (let there be $k$ of them) in the suspension. We characterize their states by $k$-fiber distribution function (see e.g. Eq. (4) in [5]) $f_k(\mathbf{r}_1, \mathbf{v}_1, \mathbf{p}_1, \mathbf{m}_1, \dots, \mathbf{r}_k, \mathbf{v}_k, \mathbf{p}_k, \mathbf{m}_k)$. By integrating $f_k$ over coordinates of $k-1$ resp. $k-2$ fibers we arrive at one-fiber distribution function

$$ f_1(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m}) \qquad (5) $$

resp. two-fiber distribution function (we use the coordinates (4))

$$ f_2(\mathbf{r}, \mathbf{v}, \mathbf{R}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2) \qquad (6) $$

Now we pass from fibers themselves to suspensions of fibers. We need to make such step since we want to discuss both morphology and rheology (i.e. we need also an expression for the extra stress tensor). We choose to regard the suspension as an incompressible isothermal fluid involving an internal structure (also called a microstructure). States of the fluid are characterized by the overall momentum field $\mathbf{u}(\mathbf{r})$ and the internal structure is characterized by either (5) or (6). Consequently, state variables of the suspension are

$$ (\mathbf{u}(\mathbf{r}), f_1(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m})) \qquad (7) $$

in the setting of one-fiber kinetic theory, and

$$ (\mathbf{u}(\mathbf{r}), f_2(\mathbf{r}, \mathbf{v}, \mathbf{R}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2)) \qquad (8) $$

in the setting of two-fiber kinetic theory. We recall that in the case of rigid spheres, (7) reduces to $(\mathbf{u}(\mathbf{r}), f_1(\mathbf{r}, \mathbf{v}))$ and (8) to $(\mathbf{u}(\mathbf{r}), f_2(\mathbf{r}, \mathbf{v}, \mathbf{R}))$. The inclusion of the overall momentum field to the state variables is essential for obtaining, as a part of a single analysis, expression for the extra stress tensor $\sigma$ (see (45) and (57)) and time evolution equations for the microstructure. The intrinsic compatibility of $\sigma$ and the microstructural equations is guaranteed.
We continue by turning our attention to reduced descriptions of the microstructure in which moments of the distribution functions replace the distribution functions themselves.

First, we replace the microscopic momentum **v** and **m** in the one-fiber distribution function by the macroscopic momentum **V**(**r**, **p**) and **M**(**r**, **p**) defined by

$$ 
\begin{aligned}
\mathbf{V}(\mathbf{r}, \mathbf{p}) &= \int d\mathbf{v} \int d\mathbf{m} v f_1(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m}) \\
\mathbf{M}(\mathbf{r}, \mathbf{p}) &= \int d\mathbf{v} \int d\mathbf{m} m f_1(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m})
\end{aligned}
\quad (9) $$

and the microscopic moments ($\mathbf{v}$, $\mathbf{w}$, $\mathbf{m}_1$, $\mathbf{m}_2$) in the two-fiber distribution function by the macroscopic moments

($\mathbf{V}(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)$, $\mathbf{W}(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)$, $\mathbf{M}_1(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)$, $\mathbf{M}_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)$) defined by

$$ 
\begin{aligned}
\mathbf{V}(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) &= \int d\mathbf{v} \int d\mathbf{w} \int d\mathbf{m}_1 \int d\mathbf{m}_2 \mathbf{v} f_2(\mathbf{r}, \mathbf{R}, \mathbf{v}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2) \\
\mathbf{W}(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) &= \int d\mathbf{v} \int d\mathbf{w} \int d\mathbf{m}_1 \int d\mathbf{m}_2 \mathbf{w} f_2(\mathbf{r}, \mathbf{R}, \mathbf{v}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2) \\
\mathbf{M}_1(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) &= \int d\mathbf{v} \int d\mathbf{w} \int d\mathbf{m}_1 \int d\mathbf{m}_2 \mathbf{m}_1 f_2(\mathbf{r}, \mathbf{R}, \mathbf{v}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2) \\
\mathbf{M}_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) &= \int d\mathbf{v} \int d\mathbf{w} \int d\mathbf{m}_1 \int d\mathbf{m}_2 \mathbf{m}_2 f_2(\mathbf{r}, \mathbf{R}, \mathbf{v}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2)
\end{aligned}
\quad (10) $$

The one and two-fiber distribution functions ($\psi_1(\mathbf{r}, \mathbf{p})$, $\psi_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)$) become

$$ \psi_1(\mathbf{r}, \mathbf{p}) = \int d\mathbf{v} \int d\mathbf{m} f_1(\mathbf{r}, \mathbf{v}, \mathbf{p}) \quad (11) $$

and

$$ 
\begin{aligned}
\psi_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) &= \int d\mathbf{v} \int d\mathbf{w} \int d\mathbf{m}_1 \\
&\quad \times \int d\mathbf{m}_2 f_2(\mathbf{r}, \mathbf{R}, \mathbf{v}, \mathbf{w}, \mathbf{p}_1, \mathbf{p}_2, \mathbf{m}_1, \mathbf{m}_2)
\end{aligned}
\quad (12) $$

We note that, due to the identity of all the fibers, $\mathbf{M}_1(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, -\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) = \mathbf{M}_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_2, \mathbf{p}_1)$. Also, we note that even if the inertia of single fibers is negligible and thus the microscopic angular momentum **m** is very small, the “collective” angular momentum **M** is not necessarily very small.

Consequently, the state variables characterizing states of fiber suspensions become

$$ (\mathbf{u}(\mathbf{r}), \psi_1(\mathbf{r}, \mathbf{p}), \mathbf{V}(\mathbf{r}, \mathbf{p}), \mathbf{M}(\mathbf{r}, \mathbf{p})) \quad (13) $$

in the one-fiber kinetic theory and

$$ (\mathbf{u}(\mathbf{r}), \psi_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{V}(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{W}(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_1(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_2(\mathbf{r}, \mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)) \quad (14) $$

in the two-fiber kinetic theory.

Later on (see Eq. (43) in the context of one-fiber kinetic theory and Eq. (55) in the context of two-fiber kinetic theory), we shall eliminate the angular momenta fields from the state variables but we need them (as will become clear in Sections 2.3.1 and 2.4) to begin the construction of the time evolution equations agreeing with mechanics and thermodynamics.

Next, we turn our attention to the coordinate **r** describing the position of the microstructure. Formation of spatial inhomogeneities induced by imposed flows in fiber suspensions have been observed and discussed (see e.g. [7,6]). Nevertheless, in this article we shall exclude from our consideration migration of fibers or formation of agglomerates, and we shall limit ourselves to spatially homogeneous suspensions. We thus omit **r** in the distribution functions and also omit the field **V** in both (13) and (14). Our sets of state variables thus simplify to

$$ (\mathbf{u}(\mathbf{r}), \psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p})) \quad (15) $$

and

$$ (\mathbf{u}(\mathbf{r}), \psi_2(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{W}(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_1(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_2(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)) \quad (16) $$

The one-fiber kinetic theory will be developed in Section 2.3 with the state variables (15). Its two-fiber extension will be however derived with a set of state variables that is still simpler than (16). We omit the coordinate **R** and the field **W**. These omissions yield the following set of state variables:

$$ (\mathbf{u}(\mathbf{r}), \psi_2(\mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_1(\mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_2(\mathbf{p}_1, \mathbf{p}_2)) \quad (17) $$

How do we pass from (16) to (17)? We simply integrate (16) over **R**. But here we need a word of caution.

In order to see more clearly the role that the position and the orientation correlations play in fiber suspensions we imagine the suspended particles to change continuously from spheres to infinitely long fibers. In other words we consider infinite number of fiber suspensions with the aspect ratio $a$ of the particles (i.e. length divided by its thickness) ranging from one to infinity. In the one extreme, when the particles are spheres (i.e. when $a = 1$), the imposed flow induces position correlations and the orientation as well as its correlations are, of course, totally absent. In the other extreme, when particles are very long fibers (i.e. when $a \to \infty$), the suspension becomes an interconnected network of fibers. In such suspensions we expect that the imposed flow induces orientation correlations more likely than the position correlations. The correlation function $\psi_2(\mathbf{p}_1, \mathbf{p}_2)$ is the probability that one fiber chosen anywhere (as far as the position coordinate is concerned) in the suspension (i.e. in the interconnected network of fibers in this case) has the orientation $\mathbf{p}_1$ provided another fiber chosen again anywhere (as far as the position coordinate is concerned) in the suspension has the orientation $\mathbf{p}_2$.

Now we move from suspensions with $a = 1$ to suspensions with $a$ that is larger than 1 but that still remains close to 1. In this family of suspensions the orientation starts to play an important role in the responses to imposed flows but the correlations will still be induced more likely in the position than in the orientation. In any case, any possible correlation in the orientation will be limited to fibers that are spatially close one to the other. On the other hand, when we move from suspensions with $a \to \infty$ to suspensions with smaller $a$ but still remaining very large, we expect that the dominant role of the orientation and the correlation in the orientation holds till the length of the fibers is comparable to the average distance between the fibers. Any change in the orientation of one chosen fiber in such suspensions is felt, due to the fiber-fiber contact interactions and also due to the strong hydrodynamic interactions among the fibers, throughout all the suspension or at least throughout a cluster (a local interconnected network of fibers) of fibers to which the chosen fiber belongs. Formation of clusters of entangled fibers has indeed been observed for example in [8–12] in suspensions of carbon nanotubes in polymer liquids. In the case of the emergence of clusters (17) are defined by averaging (16) over $\mathcal{D}_R$:

$$ 
\begin{aligned}
\psi_2(\mathbf{p}_1, \mathbf{p}_2) &= \frac{1}{V_R} \int_{\mathcal{D}_R} d\mathbf{R} \, \psi_2(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) \\
\mathbf{M}_1(\mathbf{p}_1, \mathbf{p}_2) &= \frac{1}{V_R} \int_{\mathcal{D}_R} d\mathbf{R} \, \mathbf{M}_1(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2) \\
\mathbf{M}_2(\mathbf{p}_1, \mathbf{p}_2) &= \frac{1}{V_R} \int_{\mathcal{D}_R} d\mathbf{R} \, \mathbf{M}_2(\mathbf{R}, \mathbf{p}_1, \mathbf{p}_2)
\end{aligned}
\quad (18) $$

where $\mathcal{D} \subset \mathbb{R}^3$ is the region confining the cluster, $V_R$ is its volume. However, we have to emphasize here that the emergence of clusters is the emergence of a new larger scale component of the morphology that will react with imposed flows and will influence the rheology. This new component of morphology and its influence on
rheology is absent in kinetic equations derived below. What is included in them is the morphology described in terms of (18) and the rheology associated with it. The second fiber “2” whose orientation we are comparing with the fiber “1” in (18) is in a neighborhood of the fiber “1”. The neighborhood can be the whole suspension (in the case of a global network of fibers) or the cluster. Any information about D and about the cluster morphology can, in principle, arise in a more detailed analysis in which the starting point is the choice of (16) as state variables. In this paper we limit ourselves to (17) as state variables and expect, on the basis of the physical picture described above, that our results apply to suspensions of fibers in which length of the fibers is of the same order as the average distance among them. As for the investigation with (16) as state variables, we shall limit ourselves only to indicating (in particular at the end of Section 4) the changes that would have to be made in the analysis if (17) is replaced by (16).

## 2.2. Time evolution

There are essentially two strategies that can be used to derive equations governing the time evolution of the above state variables. One, that we shall call bottom-up, begins with microscopic details of the motion of one or two fibers, the other, called top-down, begins with an overall structure of the equations guaranteeing certain properties of their solutions. Both strategies are complementary and ideally they both should be followed simultaneously. We shall now briefly recall Jeffery's bottom-up derivation of the one-fiber kinetic equation and describe the structure serving as the starting point of the top-down derivation. Subsequently, we shall re-derive with the top-down method the one-fiber orientation kinetic equation (in Section 2.3) and derive a two-fiber orientation kinetic equation (in Section 2.4). By formulating and solving the Stokes problem for a fiber immersed in a fluid, Jeffery in [4] has arrived at an equation governing the time evolution of p. A more detailed analysis of one or several fibers suspended in a fluid, developed in [2,3], has been used as a basis for direct computer simulations. These types of microhydrodynamic considerations can also serve as a basis for developing a kinetic theory of fiber suspensions. Indeed, the Jeffery equations transformed into the corresponding to them Liouville equation and supplied with an appropriate dissipative term, becomes the Jeffery kinetic equation governing the time evolution of $\psi_1(p)$. Similarly, an appropriately modified Liouville equation corresponding to equations arising in [2,3] would lead to more complete one and two-fiber kinetic equations. The bottom-up arguments based specifically on the slender-body approach to solving the Stokes problem have indeed been used to construct extensions of the Jeffery kinetic equation to semidilute suspensions in [5]. It is remarkable that the same paper also evokes (in Section 3 of [5]) some type of top-down arguments leading to the same result.

In general, the main advantage of the bottom-up approach is that the parameters entering the kinetic equations have a clear meaning in the microscopic description of the fiber-fluid and fiber-fiber interactions. Its main disadvantage is that the dissipative term in the kinetic equation and the expression for the extra stress tensor needed to transform solutions of the kinetic equation into rheological properties require separate microhydrodynamic type derivations. The mutual compatibility of the approximations used in the three separate microhydrodynamic derivations is not guaranteed. The advantages and disadvantages of the top-down approach followed below are reversed.

The top-down derivations are particularly well known in the context of the classical fluid mechanics. The required properties of solutions of the governing equations are conservations of the total mass, the total momentum, and the total energy. This requirement (together with the assumption of locality) provides a general

framework for the time evolution equations. The vector field generating the time evolution is a divergence of a flux. Specification of the fluxes (called constitutive relations) represents then the second step in the top-down derivation. In this second step the framework obtained in the first step is filled with the particular physics of the system under consideration. In order to possibly improve this top-down approach, we can ask the question as to whether there are still other (beside the global conservation) general requirements on the equations of fluid mechanics. Two such requirements have been proposed. First it is the requirement that the governing equations reflect classical mechanics of the microscopic particles composing the fluids, the second is that they reflect the equilibrium thermodynamics providing a good description of fluids at equilibrium states. As shown by Clebsch in [13] and later also by Arnold in [14] (see also [15]), the first requirement is conveniently expressed in the Hamiltonian structure of the time reversible part of the governing equations. The second requirement has been extensively studied in nonequilibrium thermodynamics [27]. Roughly speaking, this investigation introduces another field (having the physical interpretation of entropy) which is required to be a function, satisfying certain properties, of the fields constituting the state variables. The vector field generating its time evolution is required to be again divergence of a flux but now supplemented with a positive source term. In fluid mechanics itself, but in particular then when trying to formulate a structure applicable also to a larger class of mesoscopic theories, it turns out to be advantageous to put in front the Hamiltonian and the thermodynamic structures. The conservation requirements can be included as particular properties (degeneracies) of these two structures. The derivation of one and two fiber kinetic theory of fiber suspensions developed below illustrates this top-down method.

Before proceeding to the illustrations, we recall the structure expressing mathematically the compatibility with mechanics and thermodynamics. Let $x$ denote the state variables (e.g. the fields (15) or (16)). The time evolution of $x$ is governed by

$$ \frac{dx}{dt} = L E_x + \frac{\partial \Xi}{\partial S_x} \qquad (19) $$

The first term on the right hand side of (19) represents the Hamiltonian dynamics; $L$, called a Poisson bivector, expresses kinematics of $x$, $E(x)$ is a potential (a real valued function of $x$) having the physical interpretation of energy, by $E_x$ we denote the derivative of $E$ with respect to $x$. The compatibility with thermodynamics is expressed in the second term on the right hand side of (19); $\Xi$, called a dissipation potential, is a real valued function of $S_x$, where $S(x)$ is another potential having the physical interpretation of entropy. All the quantities $L$, $\Xi$, $E$, $S$ appearing in (19) are required to satisfy certain properties. Instead of listing them in the abstract setting (19), we shall present them below in the particular setting of kinetic equations discussed in this paper.

The one and two-fiber kinetic equations will be introduced in this paper as particular realizations of (19). It is useful to recognize in (19) four modules: (i) state variables $x$, (ii) the Poisson bivector $L$, (iii) the dissipation potential $\Xi$, and (iv) the potentials $E$ and $S$. To construct a particular realization of (19) means to specify the four modules. In both one and two-fiber kinetic theories we shall proceed systematically from the first to the fourth modulus.

Eq. (19), called in [16,17] GENERIC, has emerged gradually as a common structure extracted from well established (i.e. well tested with experimental observations) mesoscopic dynamical theories like for instance the classical fluid mechanics and the Boltzmann kinetic theory in [18–24]. It has been shown in [25,26] that the abstract structure (19) arises also as a natura extension of thermodynamics to mesoscopic dynamics (the time evolution generated by (19) can be seen as a continuous sequence of Legendre transformations).
## 2.3. One-fiber orientation kinetic theory

We now begin to construct particular realizations of (19). First, we turn to one-fiber kinetic theory. The first modulus, namely the state variables $x$, has already been discussed in Section 2.1. The state variables have been chosen to be (7) or alternatively (15).

Next, we turn to the second modulus, i.e. to the Hamiltonian kinematics. We begin with the time reversible and non dissipative Hamiltonian dynamics of a single fiber (in Section 2.3.1). In Section 2.3.2, we then complete the dynamics by providing its time irreversible dissipative part. In Section 2.3.3 we arrive, by eliminating the rapidly varying angular momentum field $M$, to Jeffery's equation.

We recall that a time evolution is called time reversible if there exists a transformation $\mathcal{I}$ of the state variables, called a parity transformation, that satisfies the following two properties: (i) $\mathcal{I}$ is an involution (i.e. $\mathcal{I} \circ \mathcal{I} = \text{identity transformation}$), and (ii) $\mathcal{I}$ compensates the inversion of time (i.e. the time evolution equation is invariant with respect to the application of both $\mathcal{I}$ and the inversion of time). In all the dynamical systems that will arise in this paper the parity transformation $\mathcal{I}$ consists of changing signs of all velocity-type state variables (i.e. the fluid momentum $\mathbf{u}$, the angular momentum $\mathbf{m}$, and the angular momentum field $M$ of the fibers). The exact meaning of “dissipation” is explained in Section 2.3.2.

### 2.3.1. Hamilton's dynamics

It is well known [28] that the equations governing the time evolution of (1)

$$ 
\begin{aligned}
\dot{r}_i &= E_{v_i} \\
\dot{v}_i &= -E_{r_i} \\
\dot{p}_i &= -(\mathbf{p} \times E_{\mathbf{m}})_i \\
\dot{m}_i &= -(\mathbf{m} \times E_{\mathbf{m}})_i - (\mathbf{p} \times E_{\mathbf{p}})_i
\end{aligned} 
\quad (20) 
$$

of solid mechanics can be cast into the Hamiltonian form

$$ \dot{A} = \{A, B\}^{(r, v)} + \{A, E\}^{(p, m)} \quad \text{holds for all } A \quad (21) $$

By $E(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m})$ we denote energy (a real valued function of $(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m})$); $E_{\mathbf{x}}$ is a shorthand notation for $\frac{\partial E}{\partial \mathbf{x}}$ for $\mathbf{x} = \mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m}$; $A$ is a real valued function of $(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m})$; $\times$ denotes the vector product: $(\mathbf{a} \times \mathbf{b})_i = \epsilon_{ijk}a_j b_k$, $\mathbf{a}$ and $\mathbf{b}$ are vectors in $\mathbb{R}^3$ and $\epsilon$ is the alternating tensor; $\mathbf{x}$ is the derivative of $\mathbf{x}$ with respect to the time $t$. The Poisson brackets $\{A, B\}^{(r, v)}$ and $\{A, B\}^{(p, m)}$ are given by

$$ 
\begin{aligned}
\{A, B\}^{(r, v)} &= A_{r_i} B_{v_i} - B_{r_i} A_{v_i} \\
\{A, B\}^{(p, m)} &= m_i (A_m \times B_m)_i + p_i (A_m \times B_p - B_m \times A_p)_i
\end{aligned} 
\quad (22) 
$$

We use hereafter the summation convention over the repeated indices. From the physical point of view, the Poisson bracket expresses kinematics of the state variables $(\mathbf{r}, \mathbf{v}, \mathbf{p}, \mathbf{m})$. Due to our limitation to spatially homogeneous suspensions (see Section 2.1 above) we omit hereafter the coordinates $(\mathbf{r}, \mathbf{v})$.

The verification of the equivalence of (20) and (21) proceeds as follows. The left hand side of (21) equals $A_{p_i} \dot{p}_i + A_{m_i} \dot{m}_i$, the right hand side $-A_{p_k} \epsilon_{ijk} p_i m_j + A_{m_k} [\epsilon_{ijk} m_i E_{m_k} + \epsilon_{ijk} p_i E_{p_k}]$. We then rewrite (21) into the form $A_{p_i}(\bullet)_i + A_{m_i}(\bullet \bullet)_i = 0$ which then, due to the requirement that (21) holds for all $A$, implies $(\bullet) = 0$ and $(\bullet \bullet) = 0$ that are Eq. (20).

Next, we note that also the Liouville equation

$$ \frac{\partial f_1(\mathbf{p}, \mathbf{m})}{\partial t} = -\frac{\partial}{\partial p_i}(f_1(\mathbf{p}, \mathbf{m})\dot{p}_i) - \frac{\partial}{\partial m_i}(f_1(\mathbf{p}, \mathbf{m})\dot{m}_i) \quad (23) $$

corresponding to (20) can be cast into the Hamiltonian form

$$ \dot{A} = \{A, E\}^{(f_1)} \quad \text{holds for all } A \quad (24) $$

with the Poisson bracket

$$ \{A, B\}^{(f_1)} = \int d\mathbf{p} \int d\mathbf{m} f_1 \{A_{f_1}, B_{f_1}\}^{(pm)} \quad (25) $$

The symbols $A$ and $B$ now stand for real valued functions of the distribution function $f_1(\mathbf{p}, \mathbf{m})$. The symbol $A_{f_1(p, m)}$ denotes the functional derivative of $A$ with respect to $f_1(\mathbf{p}, \mathbf{m})$. The energy $E(f_1)$ appearing in (24) is given by

$$ E(f_1) = \int d\mathbf{p} \int d\mathbf{m} f_1 E^{(\text{fiber})}(\mathbf{p}, \mathbf{m}) \quad (26) $$

where $E^{(\text{fiber})}(\mathbf{p}, \mathbf{m})$ denotes the fiber energy appearing in (21). The equivalence of (24) and (23) can be verified in the same way as the equivalence of (20) and (21) was verified in the preceding paragraph. The bracket (25) is a Poisson bracket since it depends linearly on $A_{f_1}$ and $B_{f_1}$, since $\{A, B\}^{(f_1)} = -\{B, A\}^{(f_1)}$, and since it satisfies the Jacobi identity

$$ \{\{A, B\}^{(f_1)}, C\}^{(f_1)} + \{\{B, C\}^{(f_1)}, A\}^{(f_1)} + \{\{C, A\}^{(f_1)}, B\}^{(f_1)} = 0. \text{ The first two properties are manifestly displayed, only the third property is not obvious. It can either be verified by direct (but tedious) calculations or by using a general result [28,29] relating Lie algebra structure to the Poisson bracket on its dual.} $$

Now we proceed to the reduced description that uses $(\psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p}))$ as state variables. Our objective is to use the projection (9) to derive from (23) equations governing the time evolution of $(\psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p}))$. What we require is that the projected equation represents again the Hamiltonian system. In other words, we require that the Hamiltonian nature of the time evolution is preserved in the projection. There are two routes that we can take. We can either apply the projection (9) directly on Eq. (23) or we can apply it separately on the Poisson bracket (25) and the energy (26). Taking the first route, we find that the projected equation contains, beside the moments $(\psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p}))$, also higher order moments for which we need a closure (i.e. we need to express the higher order moments in terms of the moments $(\psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p}))$. We thus face the problem of finding a closure with which the closed equation represents Hamiltonian dynamics. The second route is both easier and physically more meaningful. It turns out that the Poisson bracket (25) projects exactly without a need for closure. The energy does not, in general, projects exactly but the energy has a very clear physical meaning and we can, first of all postpone its specification, and, if it is needed, we can discuss its specification by trying to express in terms of $(\psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p}))$ the physics that is expressed in $E(f_1)$ in terms of $f_1(\mathbf{p}, \mathbf{m})$. We take below the second route.

We begin with the kinematics expressed in the Poisson bracket (25). We restrict in (25) the functions $A(f_1)$ and $B(f_1)$ to those that depend on $f_1(\mathbf{p}, \mathbf{m})$ only through their dependence on the moments $(\psi_1(\mathbf{p}), \mathbf{M}(\mathbf{p}))$. Consequently, we replace in (25) the functional derivatives $A_{f_1(p, m)}$ with $A_{\psi_1(p)} + m_i A_{M_i(p)}$ and similarly we rewrite the derivatives $B_{f_1(p, m)}$. Simple calculations then lead to

$$ \{A, B\}^{(\psi_1 M)} = \int d\mathbf{p} \left[ M_i(A_M \times B_M)_i + \psi_1 p_i \left( A_M \times \frac{\partial}{\partial \mathbf{p}}(B_{\psi_1}) - B_M \times \frac{\partial}{\partial \mathbf{p}}(A_{\psi_1}) \right)_i \right. \\ \left. + M_i p_j \left( A_M \times \frac{\partial}{\partial \mathbf{p}}(B_{M_i}) - B_M \times \frac{\partial}{\partial \mathbf{p}}(A_{M_i}) \right)_j \right] \quad (27) $$

We note that this bracket does not involve higher moments and we thus do not need any closure. Also, we note that the bracket (27) is certainly a Poisson bracket since it has arisen from the Poisson bracket (25) just by restricting the class of functions $A$ and $B$. We do not need to verify the Jacobi identity.

Hamilton's equations corresponding to the bracket (27) take the form
$$
\begin{align}
\frac{\partial \psi_1}{\partial t} &= \frac{\partial}{\partial p_i} (\psi_1 \mathbf{p} \times E_M)_i \nonumber \\
\frac{\partial M_i}{\partial t} &= -\psi_1 \left( \mathbf{p} \times \frac{\partial}{\partial \mathbf{p}} (E_{\psi_1}) \right)_i + \frac{\partial}{\partial p_j} (M_i \mathbf{p} \times E_M)_j - (\mathbf{M} \times E_M)_i \tag{28} \\
&\quad - M_j \left( \mathbf{p} \times \frac{\partial}{\partial \mathbf{p}} (E_{M_j}) \right)_i \nonumber
\end{align}
$$

where $E(\psi_1, \mathbf{M})$ is energy that we leave at this point unspecified.

So far, we have considered only fibers and not the fluid in which they are suspended. The coupling between the time evolution of the fibers and the fluid, caused by fluid-fiber interactions, will show up in our analysis only in the dissipation discussed below in Section 2.3.2. The equation governing the time reversible part of the evolution of the fluid momentum field $\mathbf{u}(\mathbf{r})$ is thus independent of $f(\mathbf{p}, \mathbf{m})$, Following [13,14,29], its time evolution is governed by Hamilton's equation

$$ \dot{A} = \{A, E\}^{(u)} \text{ holds for all } A \quad (29) $$

with the Poisson bracket

$$ \{A, B\}^{(u)} = \int d\mathbf{r} u_i (\partial_j (A_{u_i}) B_{u_j} - \partial_j (B_{u_i}) A_{u_j}) \quad (30) $$

where $\partial_i$ is a shorthand notation for $\frac{\partial}{\partial u_i}$. Written explicitly, the equation governing the time evolution of $\mathbf{u}(\mathbf{r})$ has the familiar form

$$ \frac{\partial u_i}{\partial t} = -\partial_j (u_i E_{uj}) - \partial_i \pi \quad (31) $$

where $\pi(\mathbf{r})$, the hydrostatic pressure, is given by $\pi(\mathbf{r}) = -e(\mathbf{u}; \mathbf{r}) + u_j E_{uj}$, $e(\mathbf{u}; \mathbf{r})$ is the energy density, i.e. $E(\mathbf{u}) = \int d\mathbf{r} e(\mathbf{u}; \mathbf{r})$.

Finally, for the whole suspension, the time reversible part of the time evolution is governed by

$$ \dot{A} = \{A, E\}^{(u,f_1)} \text{ holds for all } A \quad (32) $$

where

$$ \{A, E\}^{(u,f_1)} = \{A, E\}^{(f_1)} + \{A, E\}^{(u)} \quad (33) $$

when we use $(\mathbf{u}(\mathbf{r}), f_1(\mathbf{p}, \mathbf{m}))$ as state variables and by

$$ \dot{A} = \{A, E\}^{(u, \psi_1, M)} \text{ holds for all } A \quad (34) $$

where

$$ \{A, E\}^{(u, \psi_1, M)} = \{A, E\}^{(\psi_1, M)} + \{A, E\}^{(u)} \quad (35) $$

when we use (15) as state variables. In the explicit form, Eq. (32) becomes the pair of Eqs. (23) and (31) and Eq. (34) the pair of Eqs. (28) and (31). The energy $E$ has been left so far undetermined. We have to make only one requirement, namely that $E$ is invariant with respect to the parity transformation $\mathcal{I}$. This is because only with such energy the above equations are time reversible.

### 2.3.2. Dissipation

The time evolution equations arising in the previous section are compatible with mechanics (they possess the Hamiltonian structure) and with the conservation requirements (the energy is conserved since $\{E, E\} = 0$ – see (32) or (34) – and the total fluid momentum $\int d\mathbf{r} \mathbf{u}(\mathbf{r})$ is conserved since the right hand side of (31) is divergence of a flux). What remains now is to make the above time evolution equations compatible also with equilibrium thermodynamics. Since the suspensions under consideration are assumed to be isothermal, we achieve it simply by replacing in Eqs. (23) and (31) (or alternatively in Eqs. (28) and (31)) the energy $E$ with the free energy $\Phi$ and by supplementing these equations with new terms that make the free energy to decrease during the time evolution. As we have done it with the energy $E$, we also postpone as much as possible the specification of the free energy $\Phi$.

Here we only require that the free energy $\Phi$, as well as the energy $E$, is invariant with respect to the parity transformation $\mathcal{I}$.

From now on, we shall continue to develop the time evolution equations only for the state variables (15).

Following nonequilibrium thermodynamics, the new term representing the time irreversible and dissipative part of the time evolution is

$$
\begin{align}
\left(\frac{\partial u_i}{\partial t}\right)_{\text{diss}} &= -\frac{\partial \Xi}{\partial \Phi_{u_i}} \nonumber \\
\left(\frac{\partial \psi_1}{\partial t}\right)_{\text{diss}} &= -\frac{\partial \Xi}{\partial \Phi_{\psi_1}} \nonumber \\
\left(\frac{\partial M_i}{\partial t}\right)_{\text{diss}} &= -\frac{\partial \Xi}{\partial \Phi_{M_i}} \tag{36}
\end{align}
$$

where $\Xi$, called a dissipation potential, is required to satisfy the following properties:

$$
\begin{align}
\Xi \text{ is a real valued function of } (\Phi_u, \Phi_{\psi_1}, \Phi_M) \nonumber \\
\Xi \text{ reaches its minimum at } (0,0,0) \tag{37} \\
\Xi \text{ is a convex function in a neighborhood of } (0,0,0) \nonumber
\end{align}
$$

Indeed, it follows from (28), (31) and from (36) that $\dot{\Phi} = -\Phi_{u_i} \frac{\partial \Xi}{\partial \Phi_{u_i}} - \Phi_{\psi_1} \frac{\partial \Xi}{\partial \Phi_{\psi_1}} - \Phi_{M_i} \frac{\partial \Xi}{\partial \Phi_{M_i}}$ which is negative in virtue of (37).

The remaining problem is thus to construct the dissipation potential $\Xi$. Being inspired again by nonequilibrium thermodynamics we proceed by, first, introducing thermodynamic forces driving the suspension to equilibrium, and second, by constructing from them the dissipation potential.

The thermodynamic forces, denoted hereafter by the symbol $X$ (or later in Section 2.2 also by the symbol $Y$), are required to satisfy the following properties:

(i) All the thermodynamic forces must disappear at equilibrium states. Since $\Phi_u$, $\Phi_{\psi_1}$, and $\Phi_M$ equal zero at equilibrium states (because $\Phi$ reaches its minimum at such states), the thermodynamic forces will be therefore proportional to $\Phi_u$, $\Phi_{\psi_1}$, and $\Phi_M$. We shall not consider hereafter the strongly nonequilibrium situations in which higher order terms in $\Phi_u$, $\Phi_{\psi_1}$, and $\Phi_M$ would play a significant role.

(ii) In order that (36) represents the time irreversible part of the time evolution, signs of the forces change when the signs of $\mathbf{u}$ and $\mathbf{M}$ change (provided $\Xi$ is a quadratic function of the forces – see (39) below). This then means that the thermodynamic forces are proportional to $\Phi_u$, $\Phi_M$ and are independent of $\Phi_{\psi_1}$. In addition, this also means that the higher order terms that we have decided to neglect in the previous point begin with the third order.

(iii) In order to conserve the total overall momentum $\int d\mathbf{r} \mathbf{u}(\mathbf{r})$ also in the time irreversible time evolution, the first equation in (36) has to be a divergence of a flux. This then means that the thermodynamic forces cannot depend on $\Phi_u$ itself but only on $\nabla \Phi_u$. This is because $\frac{\partial \Xi}{\partial \Phi_{u_i}} = -\nabla \frac{\partial \Xi}{\partial \nabla \Phi_{u_i}}$.

(iv) Consequently, the thermodynamic forces are constructed from the following vectors: $\Phi_M$, $\Phi_u$, $\nabla$, $\mathbf{p}$. Moreover, the forces have to be vectors of the same type as $\Phi_M$ (i.e. pseudovectors). This then means that the forces constructed from $\Phi_u$, $\nabla$ and $\mathbf{p}$ have to involve the vector product. If we restrict ourselves to the vectors that depend linearly on $\Phi_u$ (see the point (ii) above), then there are two vectors that satisfy this requirement: $\nabla \times \Phi_u$ and $\mathbf{p} \times D\mathbf{p}$, where $D_{ij} = \frac{1}{2}(\partial_j(\Phi_{u_i}) + \partial_i(\Phi_{u_j}))$.

The above considerations lead to the following two forces
$$ \mathbf{X}^{(M)} = \Phi_M + \frac{1}{2} \nabla \times \Phi_u \qquad (38) $$

$$ \mathbf{X}^{(D)} = \mathbf{p} \times \mathbf{Dp} $$

These are the forces that we are proposing to be the forces driving the dissipation in the kinetic theory of fiber suspensions that uses (15) as state variables. We see indeed that $\mathbf{X}^{(M)} = 0$ implies that the rotation of the fibers follows exactly the rotation of the whole suspension, and $\mathbf{X}^{(D)} = 0$ (together with the $\mathrm{tr}\mathbf{D} = 0$ that is a consequence of the overall incompressibility) implies that the velocity gradient of the suspension (and thus also its rotation) disappears.

Next, we turn to the dissipation potential $\Xi$. Our problem is to construct from the two forces (38) a scalar satisfying the properties (37). In view of our limitation to lower orders in $\Phi_u$ and $\Phi_M$ and thus also in $\mathbf{X}^{(M)}, \mathbf{X}^{(D)}$ – see (38), the dissipation potential that we look for has the following form:

$$ \Xi = \int d\mathbf{p} (\mathbf{X}^{(M)}, \mathbf{X}^{(D)}) \psi_1 \Lambda \begin{pmatrix} \mathbf{X}^{(M)} \\ \mathbf{X}^{(D)} \end{pmatrix} \qquad (39) $$

The matrix

$$ \Lambda = \begin{pmatrix} \Lambda^{(MM)} & \Lambda^{(MD)} \\ \Lambda^{(MD)} & \Lambda^{(DD)} \end{pmatrix} \qquad (40) $$

in general a function of the fields (15), has to be positive definite in order to satisfy (37). Otherwise, our top-down arguments do not provide any more information about $\Lambda$. We can be either content with considering $\Lambda$ as material parameters (i.e. the parameters, beside those entering the free energy $\Phi$, in which the individual nature of the suspension under consideration is expressed) whose values are found by comparison with results of experimental observations, or we can turn to the comparison with Jeffery's bottom-up derivation where all parameters entering the kinetic equation have a clear meaning inside the microscopic analysis of fluid-fiber interactions.

We are now in position to write explicitly the governing Eqs. (28), (31), and (36) for the state variables (15) (in order to avoid complex notation, we shall replace the matrices $\Lambda^{(MM)}, \Lambda^{(MD)}$ and $\Lambda^{(DD)}$ by scalars $\Lambda^{(MM)}, \Lambda^{(MD)}$ and $\Lambda^{(DD)}$)

$$ \begin{aligned} \frac{\partial u_i}{\partial t} &= -\partial_j(u_i \Phi_{uj}) - \partial_i \pi - \partial_j \sigma_{ij} \\ \frac{\partial \psi_1}{\partial t} &= \frac{\partial}{\partial p_i} (\psi \mathbf{p} \times \mathbf{\Phi}_M)_i \\ \frac{\partial M_i}{\partial t} &= -\psi_1 \left( \mathbf{p} \times \frac{\partial}{\partial \mathbf{p}} (\Phi_{\psi_1}) \right)_i + \frac{\partial}{\partial p_j} (M_i \mathbf{p} \times \mathbf{\Phi}_M)_j - (\mathbf{M} \times \mathbf{\Phi}_M)_i \qquad (41) \\ & \quad - M_j \left( \mathbf{p} \times \frac{\partial}{\partial \mathbf{p}} (\Phi_{M_j}) \right)_i - \Lambda^{(MM)} \psi_1 \left( \Phi_{M_i} - \frac{1}{2} (\nabla \times \mathbf{\Phi}_u)_i \right) \\ & \quad - \Lambda^{(MD)} \psi_1 (\mathbf{p} \times \mathbf{Dp})_i \end{aligned} $$

where $\pi$ is the hydrostatic pressure given by $\pi(\mathbf{r}) = -\varphi(\mathbf{u}, \mathbf{M}; \mathbf{r}) + u_i \Phi_{ui}$; $\varphi(\mathbf{u}, \mathbf{M}; \mathbf{r})$, and the extra stress tensor $\sigma$ (we leave out its antisymmetric part) is given by

$$ \sigma_{ij} = - \int d\mathbf{p} \, \psi_1 \left[ \frac{1}{2} \Lambda^{(DD)} \left( (p_i p_k D_{jk} + p_j p_k D_{ki}) - 2p_i p_j p_k p_l D_{kl} \right) \right. \\ \left. + \Lambda^{(DM)} (p_i p_k \Omega_{kj} + p_j p_k \Omega_{ki}) + \Lambda^{(DM)} \frac{1}{2} (\epsilon_{ijk} p_i p_k \Phi_{M_i} + \epsilon_{ijk} p_j p_k \Phi_{M_j}) \right] \qquad (42) $$

where $\Omega_{ij} = \frac{1}{2}(\partial_i \Phi_{uj} - \partial_j \Phi_{ui})$.

### 2.3.3. Elimination of the field of the angular momentum M(r)

Jeffery's equation involves only the field $\psi_1(\mathbf{p})$. If we want to compare the kinetic Eq. (41) with Jeffery's equation, we have to therefore eliminate the angular momentum field $\mathbf{M}(\mathbf{p})$. We recall

that we absolutely need $\mathbf{M}(\mathbf{p})$ for formulating the nondissipative Hamiltonian dynamics. But after introducing dissipation, we can eliminate it. The elimination is based on the assumption that the evolution of $(\mathbf{M}, \psi_1)$ proceeds in two stages. In the first (fast) stage the field $\mathbf{M}$ rapidly reaches a stationary state at which it becomes enslaved to the field $\psi_1$ (i.e. becomes a function of $\psi_1$). In the second (slow) stage only the filed $\psi$ continues to evolve together with $\mathbf{M}$ that, since it does not have anymore an autonomous existence, follows passively the evolution of $\psi_1$. It is then the time evolution in the second stage that is expected to be described by the Jeffery equation.

In order to find $\mathbf{M}(\psi_1)$ that is reached after completing the fast stage in the time evolution, we turn to the third equation in (41). This equation, if we omit in it the time derivative and the terms quadratic and higher order in $\mathbf{M}$ and $\mathbf{\Phi}_M$, yields

$$ \mathbf{\Phi}_M = -\frac{1}{2} \nabla \times \mathbf{\Phi}_u - (\Lambda^{(MM)})^{-1} \left( \mathbf{\Lambda}^{(MD)} \mathbf{X}^{(D)} + \mathbf{p} \times \frac{\partial}{\partial \mathbf{p}} (\Phi_{\psi}) \right) \qquad (43) $$

This, if inserted into the second equation in (41), yields

$$ \begin{aligned} \frac{\partial \psi_1}{\partial t} &= -\frac{\partial}{\partial p_i} (\psi_1 \Omega_{ij} p_j) - \frac{\partial}{\partial p_i} (\psi_1 \mathbf{p} \times (\mathbf{\Lambda}^{(MM)})^{-1} \mathbf{\Lambda}^{(MD)} \mathbf{X}^{(D)})_i \\ & \quad - \frac{\partial}{\partial p_i} \left( \psi_1 \mathbf{p} \times (\mathbf{\Lambda}^{(MM)})^{-1} (\mathbf{p} \times \frac{\partial}{\partial \mathbf{p}} (\Phi_{\psi_1})) \right)_i \\ &= -\frac{\partial}{\partial p_i} (\psi_1 \Omega_{ij} p_j) - \frac{\partial}{\partial p_i} (\psi_1 \lambda_{il} (D_{lj} p_j - D_{jk} p_l p_j p_k)) \\ & \quad - \frac{\partial}{\partial p_i} \left( \psi_1 \Gamma_{il} (p_i p_j - \delta_{ij}) \frac{\partial}{\partial p_j} (\Phi_{\psi_1}) \right) \end{aligned} \qquad (44) $$

where $\lambda = (\Lambda^{(MM)})^{-1} \Lambda^{(MD)}$ and $\Gamma = (\Lambda^{(MM)})^{-1}$. We shall hereafter choose, for the sake of simplicity, both $\Gamma$ and $\lambda$ to be diagonal matrices with the entries $\Gamma$ and $\lambda$ respectively. Eq. (44) is now an equation that we can compare with Jeffery's equation.

We make a few observations. First, we note that (44) still involves an unspecified free energy $\Phi(\mathbf{u}, \psi_1)$. We have not yet committed ourselves to a specific fiber suspension.

Next, we note that $\lambda$ appearing in (44) is in Jeffery's equation the structure factor $\lambda = \frac{a^2-1}{a^2+1}$, where $\lambda = \lambda \delta$ and $a$ is the fiber aspect ratio (the length of the fiber divided by its diameter). This is indeed consistent with the way $\lambda$ has arisen in (44). It does not depend directly on the dissipative kinetic coefficients entering the matrix $\mathbf{\Lambda}$ (since $\lambda = (\mathbf{\Lambda}^{(MM)})^{-1} \mathbf{\Lambda}^{(MD)}$), it only waits the influence of the force $\mathbf{X}^{(D)}$ on the relaxation of the angular momentum driven by $\mathbf{X}^{(M)}$ toward equilibrium. In the case of the sphere (i.e. $a=1$) the force $\mathbf{X}^{(D)}$ does not contribute to the relaxation at all and in the case of a slender fiber (i.e. $a \to \infty$) both forces $\mathbf{X}^{(M)}$ and $\mathbf{X}^{(D)}$ are of the same importance for the relaxation of the angular momentum.

Finally, we note that the kinetic Eq. (44) has arisen already equipped with the dissipative term (the last term on the right hand side of (44)). Such term does not arise in Jeffery's microhydrodynamical investigation of the fluid-fiber interactions. Also the extra stress tensor $\sigma$ (i.e. the quantity in which the coupling of (44) to the time evolution of the overall momentum equation is expressed – the first equation in (41)) is not a part of the original Jeffery's analysis. We shall derive below an expression for $\sigma$ from the top-down arguments.

### 2.3.4. Extra stress tensor

The kinetic Eq. (44) has arisen as an intermediate step on the way leading from (41) to its solutions. We expect that solutions to (41) and to (44) are very close to each other for the time that is larger than the time needed to complete the first (fast) stage of the time evolution governed by (41). However, the Hamiltonian structure of the time reversible part of (41) is not preserved in the passage from (41)–(44). The latter does not possess the Hamiltonian structure and consequently we cannot follow the method that led us to the expression (42) for the extra stress tensor $\sigma$. We can still however use its thermodynamic content, namely: (i)
$\Phi(\mathbf{u}, \psi_1)$ remains unchanged during the reversible part of the time evolution, and (ii) $\Phi(\mathbf{u}, \psi_1)$ does not increase during the complete time evolution. From the requirement (i) we obtain (we make a simplifying assumption $\lambda = \lambda \delta$)

$$ \sigma_{ij} = -\lambda \frac{1}{2} \left[ \int d\mathbf{p} \, \psi_1 p_i \frac{\partial}{\partial p_j} (\Phi_{\psi_1}) + \int d\mathbf{p} \, \psi_1 p_j \frac{\partial}{\partial p_i} (\Phi_{\psi_1}) - 2 \int d\mathbf{p} \, \psi_1 p_i p_j p_k \frac{\partial}{\partial p_k} (\Phi_{\psi_1}) \right] \quad (45) $$

This is an expression for the extra stress tensor that corresponds to the kinetic Eq. (44). We easily verify that with (45) the following equality and inequality hold: $(\dot{\Phi})_{\text{rev}} = \int d\mathbf{r} \Phi_{\mathbf{u}} \frac{\partial \psi}{\partial t} \Big|_{\text{rev}} + \int d\mathbf{p} \Phi_{\psi_1} \frac{\partial \psi_1}{\partial t} \Big|_{\text{rev}} = 0$ and $\dot{\Phi} = \int d\mathbf{r} \Phi_{\mathbf{u}} \frac{\partial \psi}{\partial t} + \int d\mathbf{p} \Phi_{\psi_1} \frac{\partial \psi_1}{\partial t} \le 0$. By the symbol ()$_{\text{rev}}$ we denote the time reversible part of the time evolution. We refer the interested reader to [30] for details about this thermodynamic method of calculating the stress tensor and about its comparison with other methods.

### 2.3.5. Free energy

It remains to discuss the last modulus (i.e. the potentials $E$ and $S$). Since we are restricting ourselves in this paper to isothermal suspensions, it suffices to specify only one potential, namely the free energy that we shall denote by the symbol $\Phi$. We recall that the free energy

$$ \Phi(x) = E(x) - k_B T S(x) \quad (46) $$

where $x$ denotes the state variables, $E(x)$ is the energy, $k_B$ the Boltzmann constant, $T$ is the temperature (a constant due to our limitation to isothermal suspensions) and $S(x)$ the entropy. As for the energy, we take the overall kinetic energy

$$ E(x) = \int d\mathbf{r} \frac{\mathbf{u}^2}{\rho} \quad (47) $$

where $\rho$ is the overall mass density (a constant due to the assumption of the overall incompressibility), and

$$ S(x) = -n_f \int d\mathbf{p} \psi_1(\mathbf{p}) \ln \psi_1(\mathbf{p}) - B_0(n_f)^2 \int d\mathbf{p} \int d\mathbf{p}' |\mathbf{p} \times \mathbf{p}'| \psi_1(\mathbf{p}) \psi_1(\mathbf{p}') \quad (48) $$

is the Onsager entropy [1]. By $B_0$ we denote a constant proportional to the excluded volume, $n_f$ is the number density of the fibers. The first term on the right hand side of (48) is the Boltzmann entropy and the second term expresses fiber-fiber topological interactions (note that this term is the smallest if the two fibers, one with the orientation $\mathbf{p}$ and the other with the orientation $\mathbf{p}'$ are parallel). The free energy thus becomes

$$ \Phi(\mathbf{u}, \psi_1) = \int d\mathbf{r} \frac{\mathbf{u}^2}{\rho} + k_B T \left[ n_f \int d\mathbf{p} \psi_1(\mathbf{p}) \ln \psi_1(\mathbf{p}) + B_0(n_f)^2 \int d\mathbf{p} \int d\mathbf{p}' |\mathbf{p} \times \mathbf{p}'| \psi_1(\mathbf{p}) \psi_1(\mathbf{p}') \right] \quad (49) $$

Inside the equilibrium theory, the free energy represents the complete information about the suspension. For example, Onsager has shown [1] that the experimentally observed isotropy-anisotropy phase transition at equilibrium is displayed as a pitchfork bifurcation in the solutions to the equation $\Phi_{\psi_1} = 0$ appearing at a critical value for $B_0 n_f$. It is interesting to note that while the mathematical demonstration of this Onsager result is rather complex (see [1,31]) it becomes elementary on the more macroscopic level on which the distribution function $\psi_1(\mathbf{p})$ is replaced by its second moment $a_{ij} = \int d\mathbf{p} p_i p_j \psi_1(\mathbf{p})$. We shall briefly sketch the demonstration.

First, we have to replace the Onsager free energy (46)-(48) with a corresponding to it free energy on the more macroscopic level on which the tensor $\mathbf{a}$ replaces the distribution function $\psi_1(\mathbf{p})$ as the

state variable. Such free energy (we write only its part involving $\mathbf{a}$) is well known, see [32]): $\Phi(\mathbf{a}) = -k_B T n_f (\frac{1}{2} \ln \det \mathbf{a} - B_0 n_f (1 - \text{tra} \mathbf{a}))$. We note that the second term in this expression is in fact exactly the second term on the right hand side of (48) in which $|\mathbf{p}_1 \times \mathbf{p}_2|$ is replaced by $(\mathbf{p}_1 \times \mathbf{p}_2)^2$. We look now for the minimum of this free energy subjected to the constraint $\text{tra} \mathbf{a} = 1$ (expressing, on the level that uses $\mathbf{a}$ as the state variable, the constraint $|\mathbf{p}| = 1$). Without loss of generality, we look for $\mathbf{a}$ in the form of a diagonal matrix with the diagonal $(a, a, (1-2a))$. Elementary calculations involving the elimination of the Lagrange multiplier used to take into account the constraint $\text{tra} \mathbf{a} = 1$, lead to the following equation: $(a - \frac{1}{3})(4B_0 n_f a^2 - 2B_0 n_f a + 1) = 0$. Consequently, we see that the solution $a = \frac{1}{3}$, corresponding to the isotropic state, bifurcates for $B_0 n_f = 4$ into another solution corresponding to an anisotropic state.

We now briefly compare the above top-down derivation of Jeffery's kinetic equation with the bottom-up derivation based on solving the Stokes problem for a fiber suspended in a Newtonian fluid. This microhydrodynamical result enters the Jeffery equation the third line in Eq. (44) with the material parameter $\lambda$ specifically expressed in terms of the microscopic material parameters entering the microhydrodynamic analysis of the Stokes problem. In the above bottom-up derivation the same terms have arisen as a consequence of the choices (38) and (39) for the thermodynamic forces and the dissipation potential respectively. The material parameter $\lambda$ is only required to be positive but otherwise remains undetermined. The symmetry type arguments on the basis of which the forces and the dissipation potential have been chosen (see the paragraph preceding Eq. (38)) replace the arguments used in the approximative solution of the Stokes problem in the bottom-up derivation. We note here that in addition to the in addition to deriving the Jeffery rotation rate our derivation also includes the dissipative part of the kinetic equation (the last line in Eq. (44)) and the formula (45) for the extra stress tensor, both involving the free energy (49) (that is absent in Jeffery's analysis).

Before leaving the one-fiber kinetic theory, we collect all the material parameters that are involved in the governing Eqs. (44)-(46). With the simplifying assumption $\Gamma = \Gamma \delta$ and $\lambda = \lambda \delta$ that we make about the kinetic coefficients entering the dissipation potential, the parameters are:

$$ \lambda, \Gamma, B_0, n_f \quad (50) $$

The parameter $\Gamma$ is required to be positive; both $\lambda$ and $\Gamma$ can depend on $\psi(\mathbf{p})$.

### 2.4. Two-fiber orientation kinetic theory

The only place in the one-fiber orientation kinetic theory discussed above where we can take into account fiber-fiber interactions directly is in the free energy $\Phi$. We have done it in the second term on the right hand side of (48). In addition, indirectly, the fiber-fiber interactions can be expressed in an appropriately chosen dependence of the material parameters (50) on the state variable fields $(\mathbf{u}, \psi)$. The setting of the two-fiber kinetic theory (with the state variables (17)) that we shall develop in this section will give us a more opportunity to involve fiber-fiber and more complex fiber-fluid interactions.

One of the advantages of the top-down viewpoint that we are following in this paper is that we can extend in a straightforward manner the derivation of one-fiber kinetic equations to the derivation of two-fiber kinetic equations. Below, we shall leave out details and concentrate only on the places in the derivation where some changes are made.

The modifications needed to adapt Section 2.3 to the setting of two-fiber kinetic theory are quite obvious. The unit vector $\mathbf{p}$ is replaced by a pair $(\mathbf{p}_1, \mathbf{p}_2)$ of unit vectors, the distribution function
$\psi_1(\mathbf{p})$ by $\psi_2(\mathbf{p}_1, \mathbf{p}_2)$, and the field $\mathbf{M}(\mathbf{p})$ by a pair of fields $(\mathbf{M}_1(\mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_2(\mathbf{p}_1, \mathbf{p}_2))$ that are related by $\mathbf{M}_1(\mathbf{p}_1, \mathbf{p}_2) = \mathbf{M}_2(\mathbf{p}_2, \mathbf{p}_1)$. The Poisson bracket expressing the kinematics of $(\mathbf{M}_1(\mathbf{p}_1, \mathbf{p}_2), \mathbf{M}_2(\mathbf{p}_1, \mathbf{p}_2))$ is clearly a sum of two brackets (27), one with the field $\mathbf{M}_1$ and the other with the field $\mathbf{M}_2$.

The fluid-fiber interactions are brought into the investigations in Section 2.3.2 in the dissipation. It is also at this point where we shall introduce the fluid-fiber but now also fiber-fiber interactions. In addition to the thermodynamic forces

$$ 
\begin{aligned}
\mathbf{X}^{(M1)} &= \Phi_{\mathbf{M}_1} + \frac{1}{2} \nabla \times \Phi_{\mathbf{u}} \\
\mathbf{X}^{(M2)} &= \Phi_{\mathbf{M}_2} + \frac{1}{2} \nabla \times \Phi_{\mathbf{u}} \\
\mathbf{X}^{(D1)} &= \mathbf{p}_1 \times \mathbf{D}\mathbf{p}_1 \\
\mathbf{X}^{(D2)} &= \mathbf{p}_2 \times \mathbf{D}\mathbf{p}_2
\end{aligned}
\qquad (51) $$

that extend in an obvious manner the one-fiber thermodynamic forces (38) to the two-fiber setting, we introduce new forces

$$ 
\begin{aligned}
\mathbf{Y}^{(D1)} &= \mathbf{p}_2 \times \mathbf{D}\mathbf{p}_1 \\
\mathbf{Y}^{(D2)} &= \mathbf{p}_1 \times \mathbf{D}\mathbf{p}_2
\end{aligned}
\qquad (52) $$

representing both fiber-fluid and fiber-fiber interactions.

Adopting the notation $\hat{\mathbf{X}}^{(M)} = \begin{pmatrix} \mathbf{X}^{(M1)} \\ \mathbf{X}^{(M2)} \end{pmatrix}$ and similarly for the other thermodynamic forces, the dissipation potential $\Xi$ corresponding to (39) becomes

$$ \Xi = \int d\mathbf{p}_1 \int d\mathbf{p}_2 \left( \hat{\mathbf{X}}^{(M)}, \hat{\mathbf{X}}^{(D)}, \hat{\mathbf{Y}}^{(D)} \right) \psi_2 \hat{\mathbf{\Lambda}} \begin{pmatrix} \hat{\mathbf{X}}^{(M)} \\ \hat{\mathbf{X}}^{(D)} \\ \hat{\mathbf{Y}}^{(D)} \end{pmatrix} \qquad (53) $$

with the matrix $\mathbf{\Lambda}$ given by

$$ \hat{\mathbf{\Lambda}} = \frac{1}{2} \begin{pmatrix} \hat{\mathbf{\Lambda}}^{(MM)} & \hat{\mathbf{\Lambda}}^{(MD)} & \hat{\mathbf{\Lambda}}^{(MY)} \\ \hat{\mathbf{\Lambda}}^{(MD)} & \hat{\mathbf{\Lambda}}^{(DD)} & \hat{\mathbf{\Lambda}}^{(DY)} \\ \hat{\mathbf{\Lambda}}^{(MY)} & \hat{\mathbf{\Lambda}}^{(DY)} & \hat{\mathbf{\Lambda}}^{(YY)} \end{pmatrix} \qquad (54) $$

where the entries are now 6 by 6 matrices. For the sake of simplicity, we choose them hereafter to have the form $\hat{\mathbf{\Lambda}}^{(MM)} = \begin{pmatrix} \mathbf{\Lambda}^{(MM)}, & 0 \\ 0, & \mathbf{\Lambda}^{(MM)} \end{pmatrix}$ and similarly for the other entries of the matrix $\hat{\mathbf{\Lambda}}$.

Following closely the analysis of the previous section, we arrive at

$$ 
\begin{aligned}
\Phi_{\mathbf{M}_1} &= -\frac{1}{2} \nabla \times \Phi_{\mathbf{u}} - (\mathbf{\Lambda}^{(MM)})^{-1} \left( \mathbf{\Lambda}^{(MD)} \mathbf{X}^{(D1)} + \mathbf{\Lambda}^{(MY)} \mathbf{Y}^{(D1)} + \mathbf{p}_1 \times \frac{\partial}{\partial \mathbf{p}_1} (\Phi_{\psi_2}) \right) \\
\Phi_{\mathbf{M}_2} &= -\frac{1}{2} \nabla \times \Phi_{\mathbf{u}} - (\mathbf{\Lambda}^{(MM)})^{-1} \left( \mathbf{\Lambda}^{(MD)} \mathbf{X}^{(D1)} + \mathbf{\Lambda}^{(MY)} \mathbf{Y}^{(D2)} + \mathbf{p}_2 \times \frac{\partial}{\partial \mathbf{p}_2} (\Phi_{\psi_2}) \right)
\end{aligned}
\qquad (55) $$

that replaces now (43). Adopting the notation $\lambda = (\mathbf{\Lambda}^{(MM)})^{-1} \mathbf{\Lambda}^{(MD)}$, $\xi = (\mathbf{\Lambda}^{(MM)})^{-1} \mathbf{\Lambda}^{(MY)}$, and $\Gamma = (\mathbf{\Lambda}^{(MM)})^{-1}$ we arrive finally at the following two-fiber kinetic equation

$$ 
\begin{aligned}
\frac{\partial \psi_2}{\partial t} = & -\frac{\partial}{\partial p_{1i}} (\psi_2 \Omega_{ij} p_{1j}) \\
& -\frac{\partial}{\partial p_{1i}} (\psi_2 \lambda_{il} (D_{ij} p_{1j} - D_{jk} p_{1l} p_{1j} p_{1k})) \\
& -\frac{\partial}{\partial p_{1i}} (\psi_2 \xi_{il} (D_{ij} p_{1j} - D_{jk} p_{2l} p_{2j} p_{1k})) \\
& -\frac{\partial}{\partial p_{2i}} (\psi_2 \Omega_{ij} p_{2j}) \\
& -\frac{\partial}{\partial p_{2i}} (\psi_2 \lambda_{il} (D_{ij} p_{2j} - D_{jk} p_{2l} p_{2j} p_{2k})) \\
& -\frac{\partial}{\partial p_{2i}} (\psi_2 \xi_{il} (D_{ij} p_{2j} - D_{jk} p_{1l} p_{1j} p_{2k})) \\
& -\frac{\partial}{\partial p_{1i}} (\psi \Gamma_{il} (p_{1l} p_{1j} - \delta_{lj}) \frac{\partial}{\partial p_{1j}} (\Phi_{\psi_2})) \\
& -\frac{\partial}{\partial p_{2i}} (\psi \Gamma_{il} (p_{2l} p_{2j} - \delta_{lj}) \frac{\partial}{\partial p_{2j}} (\Phi_{\psi_2}))
\end{aligned}
\qquad (56) $$

We shall hereafter choose: $\lambda = \lambda \delta$ and $\xi = \xi \delta$.

Before continuing to discuss the extra stress tensor and the free energy, we make two comments about the new thermodynamic forces $\mathbf{Y}^{(D)}$ (see (52)).

First, we recall what we know about the physical interpretation of the forces $\mathbf{X}^{(D)}$. From the way they have been introduced, they represent the friction forces originating in the fiber-fluid interactions. This is then also confirmed indirectly by comparing the top-down derivation developed in this paper with Jeffery's microhydrodynamic derivation. Now we turn to the forces $\mathbf{Y}^{(D)}$. From the way they have been introduced, they represent again the friction force originating in the fiber-fluid interactions but now in a way in which also fiber-fiber interactions participate. It is the simplest thermodynamic force involving both fibers "1" and "2". Another type of arguments supporting the choice (52) of the forces $\mathbf{Y}^{(D)}$ comes from comparing the above analysis with the analysis reported in [5]. We note that the first three lines in the kinetic Eq. (56) is the right hand side of the Liouville equation corresponding to Eq. (38) in [5] with $\mathbf{p}_1$ replacing $\mathbf{p}$ in [5] and $\mathbf{p}_2$ replacing $\mathbf{n}$ in [5]. Eq. (38) in [5] represents a modification of the Jeffery equation caused by replacing the isotropic medium in which the fiber characterized by the unit vector $\mathbf{p}$ with a linear anisotropic medium characterized by the unit vector $\mathbf{n}$. This equation has been derived in [5] in both bottom-up manner by solving (in the slender-body approximation) the corresponding Stokes problem and in the bottom-up manner using certain symmetries and restriction to the linear dependence on the velocity gradient. From the comparison between the kinetic Eqs. (56) and (38) in [5] we thus see that in our analysis the second fiber influences the rotation of the first fiber by making the medium in which the first fiber moves anisotropic. Another way in which both fibers influence their motion will appear the free energy (58) entering the dissipative part of the time evolution.

Second, we note in the context of our derivation of the kinetic Eq. (56) we could replace the force $\mathbf{Y}^{(D1)}$ with the force $\mathbf{Y}^{(D2)}$ and the force $\mathbf{Y}^{(D2)}$ with $\mathbf{Y}^{(D1)}$. There is indeed no reason in the context of our derivation to insists on the order that we have made in the definition (52). If such change is made, the kinetic Eq. (56) changes. The third line changes into $-\frac{\partial}{\partial p_{1i}} (\psi_2 \xi_{il} (D_{ij} p_{2j} - D_{jk} p_{1l} p_{1j} p_{2k}))$ and the sixth line into $-\frac{\partial}{\partial p_{1i}} (\psi_2 \xi_{il} (D_{ij} p_{2j} - D_{jk} p_{1l} p_{1j} p_{2k}))$. Here we have to use the argument developed in [5] to prefer (56) to the modified kinetic equation. The argument is that if $\mathbf{p}_1 \mathbf{D} \mathbf{p}_1 = 0$ (or alternatively $\mathbf{p}_2 \mathbf{D} \mathbf{p}_2 = 0$) then the second fiber does not influence the rotation of the first fiber. This observation has arisen in [5] as a result of the analysis (in the slender-body approximation) of the Stokes problem in a linear anisotropic medium.

We turn now to the extra stress tensor. The same argument as the one used in Section 2.3.4 leads now to the following expression:

$$ 
\begin{aligned}
\sigma_{ij} = & -\lambda \frac{1}{2} \int d\mathbf{p}_1 \int d\mathbf{p}_2 \psi_2 \left[ p_{1i} \frac{\partial}{\partial p_{1j}} (\Phi_{\psi_2}) + p_{1j} \frac{\partial}{\partial p_{1i}} (\Phi_{\psi_2}) \right. \\
& \left. + p_{2i} \frac{\partial}{\partial p_{2j}} (\Phi_{\psi_2}) + p_{2j} \frac{\partial}{\partial p_{2i}} (\Phi_{\psi_2}) \right. \\
& \left. - 2p_{1i} p_{1j} p_{1k} \frac{\partial}{\partial p_{1k}} (\Phi_{\psi_2}) - 2p_{2i} p_{2j} p_{2k} \frac{\partial}{\partial p_{2k}} (\Phi_{\psi_2}) \right] \\
& - \xi \frac{1}{2} \int d\mathbf{p}_1 \int d\mathbf{p}_2 \psi_2 \left[ p_{1i} \frac{\partial}{\partial p_{1j}} (\Phi_{\psi_2}) + p_{1j} \frac{\partial}{\partial p_{1i}} (\Phi_{\psi_2}) \right. \\
& \left. + p_{2i} \frac{\partial}{\partial p_{2j}} (\Phi_{\psi_2}) + p_{2j} \frac{\partial}{\partial p_{2i}} (\Phi_{\psi_2}) \right. \\
& \left. - (p_{1i} p_{2j} + p_{2i} p_{1j}) p_{2k} \frac{\partial}{\partial p_{1k}} (\Phi_{\psi_2}) - (p_{1i} p_{2j} + p_{2i} p_{1j}) p_{1k} \frac{\partial}{\partial p_{2k}} (\Phi_{\psi_2}) \right]
\end{aligned}
\qquad (57) $$

Finally, we specify the two-fiber free energy. The expression (49), if rewritten in terms of the two-fiber distribution function, becomes
$$
\Phi(\mathbf{u}, \psi_2) = \int d\mathbf{r} \frac{\mathbf{u}^2}{\rho} + k_B T n_f \left[ \int d\mathbf{p}_1 \int d\mathbf{p}_2 \psi_2(\mathbf{p}_1, \mathbf{p}_2) \ln \psi_2(\mathbf{p}_1, \mathbf{p}_2) \right. \\
\left. + B_0 n_f \int d\mathbf{p}_1 \int d\mathbf{p}_2 |\mathbf{p}_1 \times \mathbf{p}_2| \psi_2(\mathbf{p}_1, \mathbf{p}_2) \right] \quad (58)
$$

We note that $\frac{\partial}{\partial p_{1k}}(\Phi_{\psi_2}) = \frac{\partial}{\partial p_{2k}}(\Phi_{\psi_2})$.

We make now somewhat surprising but in fact very well known
observation. The free energy (58) is at the same time more detailed
than (49) (since (58) involves the two-fiber distribution function
and (49) the one-fiber distribution function) and less able to de-
scribe at least some macroscopic phenomena. In particular, the
equilibrium isotropy–anisotropy phase transition is not displayed
in (58) in the same way as in (49). Indeed, the equation $\Phi\psi = 0$ with
$\Phi$ given in (58) has only one solution for all $B_0 n_f$. This observation, a
sort of uncertainty principle between the amount of microscopic
details and the ability to see overall features, is well known in
Gibbs equilibrium statistical mechanics. For example, it is well
known that phase transitions can be seen in Gibbs’ equilibrium sta-
tistical mechanics in the same way as they are seen (and in fact de-
fined) in the equilibrium thermodynamics only after the Gibbs
theory has been modified. The modification consists for example
in carrying the Gibbs theory to the so called thermodynamic limit
[33] or in bringing into it some complex geometry (nonlinearities)
by making restrictions in the spaces of distribution functions (see
[34]).

Summing up, the kinetic Eq. (56), the expression for the extra
stress tensor (57), and the free energy (58) are the governing equa-
tions of the two-fiber kinetic theory of fiber suspensions. The
material parameters entering them (i.e. the parameters through
which the individual nature of the suspension is expressed in
them) are:

$$
\lambda, \xi, \Gamma, B_0, n_f \tag{59}
$$

The advection of the fibers is influenced by the shape of the
fibers (expressed in the parameter $\lambda$) and by the fiber–fiber
interactions (expressed in the parameter $\xi$). The fiber–fiber
interactions entering the free energy (in the Onsager term) is ex-
pressed in the parameter $B_0$. If $B_0 = 0$ and $\xi = 0$ then the two-fiber
kinetic theory developed in this paper reduces to the Jeffery one-
fiber kinetic theory. This is also demonstrated below in the
numerical solutions.

**3. Numerical solution of the two-fiber orientation kinetic equations**

For numerical resolution, the time has been made dimension-
less by multiplying it by the norm of the velocity gradient tensor.
Then the velocity gradient tensor has been normalized. Finally
the stresses has been made dimensionless by dividing them by
$k_B T$. The value of $n_f$ has been put equal to one in all numerical
calculations.

The two-fiber distribution function $\psi$ depends in 2D only on two orientation angles $(\phi_1, \phi_2)$ and the time. The numerical solutions in 2D can be thus found in the framework of the classical finite element discretization.

On the other hand, in 3D, the two-fiber distribution function de-
pends on four angles and the time. The classical finite element dis-
cretization becomes in this case practically unfeasible. We turn
therefore to the Proper Generalized Decomposition (PGD) method
that is briefly described in the Appendix. In the framework of the
PGD method, the numerical solution of the two-fiber kinetic equa-
tion is searched in the form of a sum of tensor products of a pair of
functions each of them depending on two angles.

3.1. 2D kinetic equations

In 2D the unit vectors **p**₁ and **p**₂ become

$$
p_i = (\cos \phi_1, \sin \phi_i)^T; \quad i = 1, 2 \tag{60}
$$

The kinetic Eq. (56) can be rewritten into the following form:

$$
\frac{\partial \psi}{\partial t} + T_0(\phi_1, \phi_2)\psi + T_1(\phi_1, \phi_2)\frac{\partial \psi}{\partial \phi_1} + T_2(\phi_2, \phi_1)\frac{\partial \psi}{\partial \phi_2} - \Gamma \frac{\partial^2 \psi}{\partial \phi_1^2} + \frac{\partial^2 \psi}{\partial \phi_2^2} = 0 \quad (61)
$$

with

$$
\begin{align}
T_0(x, y) &= E_0(x) + E_0(y) + G_0(x, y) + G_0(y, x) + H_0(x, y) \nonumber \\
T_1(x, y) &= E_1(x) + G_1(x, y) + H_1(x, y) \tag{62}
\end{align}
$$

and

$$
\begin{align*}
E_0(x) &= \lambda((\kappa_{11} - \kappa_{22}) \cos(2x) + (\kappa_{21} + \kappa_{12}) \sin(2x)) \\
E_1(x) &= (\kappa_{21} - \kappa_{12} + \lambda((\kappa_{21} + \kappa_{12}) \cos(2x) + (\kappa_{22} - \kappa_{11}) \sin(2x))) / 2
\end{align*}
$$

$$
G_0(x, y) = -\frac{1}{2}\xi(\cos x(2\kappa_{11} \cos y + (\kappa_{21} + \kappa_{12}) \sin y) \\
\qquad \qquad + \sin x((\kappa_{21} + \kappa_{12}) \cos y + 2\kappa_{22} \sin y))
$$

$$
G_1(x, y) = -\frac{1}{2}\xi(\sin x(2\kappa_{11} \cos y + (\kappa_{21} + \kappa_{12}) \sin y) \\
\qquad - \cos x((\kappa_{21} + \kappa_{12}) \cos y + 2\kappa_{22} \sin y))
$$

$$
\begin{align}
H_0(x, y) &= 2B_0 \Gamma |\sin(x - y)| \nonumber \\
H_1(x, y) &= -B_0 \Gamma |\sin(x - y)| \cot(x - y) \tag{63}
\end{align}
$$

where $\kappa$ is the velocity gradient. We see that the functions $G$ and $H$
disappear if $\xi = 0$ and $B_0 = 0$ (i.e. when the fiber–fiber interactions
are absent).

The domain $\Omega$ of the two-fiber distribution function is given by:
$\Omega = [0,2\pi] \times [0,2\pi]$. The problem is formulated in the finite element framework using a weighting function $\psi^*$:

$$
\int_{\Omega} \psi^* \left[ \frac{\partial \psi}{\partial t} + T_0(\phi_1, \phi_2)\psi + T_1(\phi_1, \phi_2) \frac{\partial \psi}{\partial \phi_1} + T_2(\phi_2, \phi_1) \frac{\partial \psi}{\partial \phi_2} - \Gamma \frac{\partial^2 \psi}{\partial \phi_1^2} + \frac{\partial^2 \psi}{\partial \phi_2^2} \right] d\Omega = 0 \quad (64)
$$

The domain $\Omega$ is partitioned into a collection of non-overlapping finite elements (rectangular facets). A bilinear, continuous interpolation of the distribution function is then built in each finite element:

$$
\psi^e = \sum_{i=1}^{4} N_i(\phi_1, \phi_2) \psi_i^e \quad (65)
$$

where $\psi_i^e$ are the values at node $i$ and $N_i(\phi_1, \phi_2)$ is the associated
shape function which takes a unit value at the node $i$ and vanishes
at the other nodal positions.

By parts integration of the resulting equation (taking into ac-
count that the domain is unbounded) yields

$$
\int_{\Omega} \psi^* \left[ \frac{\partial \psi}{\partial t} \right] d\Omega + \int_{\Omega} \psi^* \left[ T_0(\phi_1, \phi_2)\psi + T_1(\phi_1, \phi_2) \frac{\partial \psi}{\partial \phi_1} + T_2(\phi_2, \phi_1) \frac{\partial \psi}{\partial \phi_2} \right] \\
+ \Gamma \left[ \frac{\partial \psi^*}{\partial \phi_1} \frac{\partial \psi}{\partial \phi_1} + \frac{\partial \psi^*}{\partial \phi_2} \frac{\partial \psi}{\partial \phi_2} \right] d\Omega = 0 \quad (66)
$$

Due to the advection–diffusion character of the equation, an
appropriate stabilization of the finite element scheme is needed
in order to avoid numerical instabilities induced by the convection
term. Stabilization is achieved by the usual upwinding formulation
(SUPG) which modifies the weighting functions $\psi^*$.
The discrete system

$$N_{ij} \frac{\partial \psi_j}{\partial t} + A_{ij} \psi_j = 0 \qquad (67)$$

is then solved by using the implicit time integration scheme constrained by the normality condition (that is verified at every time step).

### 3.2. 3D kinetic equations

The domain $\Omega$ is in 3D the tensorial product of two unit spheres $S_{p_1}$ and $S_{p_2}$. Due to the large number of the variables in $\Omega$ we turn to the PGD method (see Appendix).

We are looking for a solution of (56) in the form

$$\Psi(p_1, p_2) = \sum_{i=1}^{n_F} F_1^i(p_1) \cdot F_2^i(p_2) \qquad (68)$$

For numerical calculations we need to write (68) in a discrete form that uses the nodal values of each function:

$$\Psi = \sum_{j=1}^{n_F} F_1^j \otimes F_2^j \qquad (69)$$

The weak form of Eq. (56) transformed into a discrete form becomes

$$\Psi^T \mathcal{N} \dot{\Psi} + \Psi^T A \Psi = 0 \qquad (70)$$

with

$$A = \sum_{j=1}^{n_A} \underline{A}_1^j \otimes \underline{A}_2^j \qquad (71)$$

$$\mathcal{N} = \underline{N}_1 \otimes \underline{N}_2 \qquad (72)$$

and

$$\underline{N}_1 = \int_{S_{p_1}} \mathbf{N} \mathbf{N}^T dS_{p_1} \qquad (73)$$

$$\underline{N}_2 = \int_{S_{p_2}} \mathbf{N} \mathbf{N}^T dS_{p_2} \qquad (74)$$

where $\mathbf{N}$ are the vectors containing the shape functions associated to the unit sphere. As for the operator $A$, we only write the first term:

<figure><img src="image_1.png" /></figure>

$$\begin{aligned} \underline{A}_1^1 &= \int_{S_{p_1}} \mathbf{N} \frac{\partial}{\partial p_{1i}} (\Omega_{ij} p_{1j} \mathbf{N}^T) dS_{p_1} \\ \underline{A}_2^1 &= \int_{S_{p_2}} \mathbf{N} \mathbf{N}^T dS_{p_2} \end{aligned} \qquad (74)$$

The other terms are calculated in the same way.

When the system is completely specified we apply the resolution strategy of the PGD formalism. For more details about the PGD method we refer to Appendix.

### 3.3. Numerical solutions

The imposed velocity gradient in 2D is

$$\mathbf{K} = \begin{pmatrix} 0 & G \\ 0 & 0 \end{pmatrix} \qquad (75)$$

and in 3D

$$\mathbf{K} = \begin{pmatrix} -E & 0 & G \\ 0 & E & 0 \\ 0 & 0 & 0 \end{pmatrix} \qquad (76)$$

By using the implicit time integration scheme we calculate at each time step the fiber distribution function. The initial state is given by the isotropic orientation (i. e. the constant value $1/2\pi$ in 2D and $1/4\pi$ in 3D). At each time step we calculate the stress by using the expressions derived in the text.

Results in 2D (both the distribution function and the stresses) are depicted in Figs. 1–3. In Fig. 1 we show results for two-fiber kinetic theory with $B_0 = 0$ and $\xi = 0$. The distribution function shown in the figure corresponds to the stationary state. In order to see the fiber–fiber correlations we calculate from the distribution function $\psi(\mathbf{p}_1, \mathbf{p}_2)$ a correlation function $g(\mathbf{p}_1, \mathbf{p}_2)$ defined by

$$\psi(\mathbf{p}_1, \mathbf{p}_2) = \Psi(\mathbf{p}_1) \Psi(\mathbf{p}_2) g(\mathbf{p}_1, \mathbf{p}_2) \qquad (77)$$

where

$$\Psi(\mathbf{p}_1) = \int d\mathbf{p}_2 \psi(\mathbf{p}_1, \mathbf{p}_2), \qquad (78)$$

The calculated correlation function turns out to be independent of $(\mathbf{p}_1, \mathbf{p}_2)$ and equals to one. This means that in the absence of the coupling terms in the two-fiber kinetic equation (i.e. the terms, proportional to $B_0$ and $\xi$, that represent in the kinetic equation the fiber–fiber interactions) the two-fiber kinetic theory reduces

<figure><img src="image_2.png" /><figcaption>Fig. 1. Two-fiber kinetic theory in 2D with $\xi = B_0 = 0$. There are no fiber–fiber correlations in the orientation and the stresses are the same as in the one-fiber kinetic theory.</figcaption></figure>
to the one-fiber kinetic theory. The two-fiber distribution function is separable, i.e. it is a product of two one-fiber distribution functions.

In Figs. 2 and 3 we show results for the kinetic equation that takes into account the fiber-fiber interactions in the advection (the case $\xi \neq 0$) and in the free energy (the case $B_0 \neq 0$). In the depicted two-fiber orientation distribution function (corresponding to the stationary state) we see clearly a loss of the separability. We see emergence of fiber-fiber correlations in the orientation. The correlation function $g(p_1, p_2)$ is shown on Fig. 3. We also note

that the fiber-fiber interactions bring about more complex rheology with overshoots.

In the remaining two Figs. 4 and 5 we show results in 3D. In Fig. 4 we again demonstrate that for $\xi = 0$ and $B_0 = 0$ the two-fiber distribution function is separable (i.e. the orientation fiber-fiber correlations are absent) and the rheological results are the same in the one-fiber kinetic theory (calculated and displayed in [30]).

The case of $\xi \neq 0$, $B_0 \neq 0$ is represented in Fig. 5. From the depicted distribution function we see emergence of correlations in the orientation. We also note that the fiber-fiber interactions

<figure><img src="image_3.png" /><figcaption>Fig. 2. Two-fiber kinetic theory in 2D with (i) $\xi \neq 0, B_0 = 0$, (ii) $\xi \neq 0, B_0 \neq 0$, (iii) $\xi = 0, B_0 \neq 0$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Correlation function in the two-fiber kinetic theory in 2D with (i) $\xi \neq 0, B_0 = 0$, (ii) $\xi \neq 0, B_0 \neq 0$, (iii) $\xi = 0, B_0 \neq 0$.</figcaption></figure>

expressed in the value $\xi = 0.5$ increase the stress tensor by a factor of, approximatively, 2.

**4. Concluding remarks**

Two types of interactions in fiber suspensions are shown to give rise to flow induced fine-scale orientation ordering. One such interaction is the fiber–fiber interaction introduced by Onsager in his

equilibrium analysis of liquid crystals [1]. The other is a combined fiber–fluid and fiber–fiber interaction influencing fiber rotations. The limitation of the investigation reported in this paper on the orientation order (leaving aside a possible fine-scale spatial order and fiber migration) prevents us from clearly identifying the types of suspensions in which the two interactions become important and consequently in which the fine-scale order in the orientation emerges. Physical interpretation of the two interactions suggest that these are the suspensions in which the average length of the fibers is of the same order as the average distance between the centers of mass of neighbor fibers.

We have demonstrated that the flow induced fine-scale order in the orientation joins then the flow induced overall orientation of fibers to constitute the morphology that is strongly coupled to the rheology of fiber suspensions. The main results of this paper are the formulation of the two-fiber orientation kinetic Eq. (56), formulation of the corresponding to it expression (57) for the extra stress tensor, and the following results concerning solutions of the kinetic equations:

(i) As $t \to \infty$, solutions to (56) and (31) describing externally unforced fiber suspensions approach equilibrium states at which their behavior is well described by equilibrium thermodynamics. The fundamental thermodynamic relation of the equilibrium theory is determined by the free energy entering (56) evaluated at the equilibrium states.

(ii) In the case of the imposed simple shear flow, numerical solutions of the two-fiber kinetic equation in 2D and 3D show that in the absence of the fiber–fluid interactions influencing the relaxations of fiber rotation (i.e. if $\xi = 0$) and the fiber–fiber interactions in the free energy (i.e. if $B_0 = 0$) the two-fiber orientation distribution function is separable (i.e. it is a product of two one-fiber orientation distribution functions) and the implied rheology is the same as in the one-fiber orientation kinetic theory. The separability of the orientation distribution function means the absence of the fine-scale order in the orientation of the fibers. In the case when $\xi \neq 0$ and/or $B_0 \neq 0$ the two-fiber orientation distribution function ceases to be separable (i.e. the fine scale order in the orientation emerges) and the implied rheology is more complex then in the absence of the fiber–fiber interactions. The main challenge in the investigation of numerical solutions of the kinetic equations introduced in this paper is the large number of independent variables (five in the case of Eq. (56) in 3D) of the unknown functions. This type of challenge was exactly the motivation for developing the Proper Generalized Decomposition (PGD) method [35–38] – see also Appendix). The method combines in an appropriate way the standard methods of discretization with a numerical viewpoint of the method of separation of variables.

What are other types of investigations of fiber suspensions with which these results could be compared? We suggest three: direct numerical simulations, experimental observations, and other types of modeling of fiber suspensions.

As for direct simulations, the point of departure could be the same physics the one we used as the starting point of our derivation of kinetic equations (expressed however on a more microscopic level suitable for formulating governing equations of direct simulations). Appropriate statistical averages of calculated trajectories of suspended fibers would then provide an information about pair correlations that could be compared with the results obtained in this paper. Unfortunately, we cannot do such comparison since we are not aware of any simulation in which orientation correlations are calculated.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Two-fiber kinetic theory in 3D with $\xi = B_0 = 0$. There are no fiber-fiber correlations in the orientation and the stresses are the same as in the one-fiber kinetic theory.</figcaption></figure>

Next, we turn to experimental observations. As we have already noted, we do not know about any and, at least at this stage, we cannot even suggest any experimental observations that would provide a direct information about orientation correlations in fiber suspensions. Results of rheological observations do reflect, as it is shown in this paper, the correlations but they also reflect many other aspects of fiber suspensions. Comparison of the observed rheological behavior with the behavior predicted by the two-fiber kinetic equation derived in this paper does not therefore constitute a confirmation of the presence or absence of orientation correlations.

The investigation reported in this paper has a closest relation to the investigation (reported in [5–7]) of extensions of Jeffery’s single-fiber dynamics to more concentrated suspensions. We have shown that their modification of Jeffery’s dynamics, if put into the context of two-fiber kinetic theory, represents one of the sources (in addition to the Onsager entropy) of orientation correlations.

The one and two-fiber kinetic equations are derived in this paper by using the top-down method. The point of departure is an overall framework for the governing equations, called GENERIC, guaranteeing agreement of their solutions with certain experimental observations. The derivation consists of filling the framework

(or in other words, constructing its particular realizations) with the particular physics involved in the system under consideration. As for the framework itself, its physical content is, on the one hand, compatibility with the microscopic mechanics of the particles composing the macroscopic system, and, on the other hand, compatibility with the observations made in equilibrium thermodynamics. Specifically, the principal elements (modules) of the structure are: (i) state variables, (ii) their kinematics expressed in a Poisson bracket, (iii) thermodynamic forces and dissipation potential, (iv) energy and entropy potentials. Two particular realizations of the structure have been developed in this paper. In the one-fiber orientation kinetic theory the state variables are (15), the Poisson bracket is (27), the thermodynamic forces are (38) and the thermodynamic potential is (39). In the two-fiber orientation kinetic theory the state variables are (17), the Poisson bracket is a sum of two brackets (27), the thermodynamic forces are (51) and (52) and the thermodynamic potential is (53). After writing the time evolution equations (with still unspecified energy and entropy potentials), we consider the momentum fields $M$ as fast variables, assume that they have already reached their stationary states and become thus functions of the correlation functions $\psi$. The physical arguments on which the choice of the thermodynamic
<figure><img src="image_6.png" /><figcaption>Fig. 5. Two-fiber kinetic theory in 3D with $\xi \neq 0, B_0 = 0$.</figcaption></figure>

forces is made also supported by both microhydrodynamic and phenomenological arguments developed in [5].

Another way to approach the derivation of kinetic equations would be the bottom-up method in which the starting point is a detailed microhydrodynamic viewpoint of suspensions that is then gradually reduced to the mesoscopic level of kinetic theory. In both

top-down and bottom-up methods the physical basis is microscopical, but in the top-down method we are deriving separately the moduli (that are subsequently put together, according to (19) to form the complete time evolution equation together with the corresponding to it expression for the stress tensor) while in the bottom-up method the complete time evolution equation and
the expression for the stress tensor are derived directly. The mod-
ular nature of the top-down method can be seen as a guidance to
make reductions and approximations that are compatible among
themselves and with mechanics and thermodynamics.

Finally, we indicate one possible continuation of the analysis
presented in this paper. In order to be able to investigate theoret-
ically the relative importance of correlations in the position and in
the orientation of fibers in fiber suspensions a complete position
and orientation kinetic theory have to be developed. We intend
to do it in a future paper The same method that we have used in
this paper to introduce one-fiber and two-fiber orientation kinetic
theories can be used to develop the complete theory. In the reali-
zation of GENERIC corresponding to the complete position and ori-
entation kinetic theory we choose the state variables (16). Their
kinematics (i.e. the Poisson bracket) can be derived from mechan-
ics of one and two fibers (but now involving also the motion of the
centers of mass of the fibers) in the same way as we have derived
in this paper the bracket (27). In addition to the thermodynamic
forces (51) and (52), new forces involving F and D have to be
introduced. Details of the physics taking place in the suspension
will enter into the specification of the matrix A in the dissipation
potential and finally in the specification of the entropy and the en-
ergy potentials. In the case of suspensions of spherical particles
this program has been realized recently in [39-41].

Appendix A

Let

$$
\mathcal{A} \Psi = \mathcal{B} \tag{79}
$$

where

$$
A = \sum_{j=1}^{n_A} \underline{A}_1^j \otimes \underline{A}_2^j \otimes \dots \otimes \underline{A}_N^j \quad (80)
$$

$$
\mathcal{B} = \sum_{j=1}^{n_B} \underline{B}_1^j \otimes \underline{B}_2^j \otimes \dots \otimes \underline{B}_N^j \quad (81)
$$

be the equation to be solved for the unknown Ψ.

We look for solutions of (79) in the form of a tensor product

$$
\Psi = \sum_{j=1}^{\infty} \underline{F}_1^j \otimes \underline{F}_2^j \otimes \dots \otimes \underline{F}_N^j \qquad (82)
$$

Next, we use the optimization strategy consisting of minimizing
the error functional

$$
\mathcal{J}(\Psi) = \|\mathcal{A}\Psi - \mathcal{B}\|_2 \quad (83)
$$

This error function can be rewritten into the form

$$
\mathcal{J}(\Psi) = (\Psi - \mathcal{A}^{-1}\mathcal{B})^T \mathcal{A}^T \mathcal{A} (\Psi - \mathcal{A}^{-1}\mathcal{B}) = \|\Psi - \mathcal{A}^{-1}\mathcal{B}\|_{\mathcal{A}^T \mathcal{A}} \quad (84)
$$

The equation $\partial_t \Psi = 0$ then provides a new form

$$
\partial_t \mathcal{A}^T \mathcal{A} \Psi = \partial_t \mathcal{A}^T \mathcal{B} \quad (85)
$$

for the equation to be solved.

The adaptive strategy that we use for its resolution consists of
finding at every $(n_F + 1)$th iteration the best new tensor product $\Psi_R$

$$
\Psi = \underbrace{\sum_{j=1}^{n_F} \underline{F}_1^j \otimes \underline{F}_2^j \otimes \dots \otimes \underline{F}_N^j}_{\Psi_F} + \underbrace{\underline{R}_1 \otimes \underline{R}_2 \otimes \dots \otimes \underline{R}_N}_{\Psi_R} \quad (86)
$$

We proceed by making sub-iterations (i.e. iterations made in-
side a single iteration from $n_F$ to $(n_F + 1)$). Let all functions $R_1 = \underline{R}_1$, $-1 = \underline{R}_1, \dots, \underline{R}_{j-1} = \underline{R}_{j-1}, \underline{R}_{j+1} = \underline{R}_{j+1}, \dots, \underline{R}_N = \underline{R}_N$ in (87) be already
known. Our problem is to find the single unknown $R_j$. For this pur-
pose we choose the test function

$$
\Psi^* = \underline{R}_1 \otimes \dots \otimes \underline{R}_{j-1} \otimes \underline{R}_j^T \otimes \underline{R}_{j+1} \otimes \dots \otimes \underline{R}_N \quad (87)
$$

Consequently, the system to be solved writes

$$
\Psi^* A^T A \Psi_R + \Psi^* A^T A \Psi_F = \Psi^* A^T B \quad (88)
$$

which with (88) becomes

$$
\sum_{i=1}^{n_F} \sum_{k'=1}^{n_A} \sum_{k=1}^{n_A} \left( \underline{A}_j^{k'} \underline{A}_j^k \underline{F}_j^i \prod_{\substack{h=1 \\ h \neq j}}^{N} \underline{R}_h^T \underline{A}_h^{k'} \underline{A}_h^k \underline{F}_h^i \right) + \\
\sum_{k'=1}^{n_A} \sum_{k=1}^{n_A} \left( \underline{A}_j^{k'} \underline{A}_j^k \underline{R}_j \prod_{\substack{h=1 \\ h \neq j}}^{N} \underline{R}_h^T \underline{A}_h^{k'} \underline{A}_h^k \underline{R}_h \right) = \\
\sum_{k'=1}^{n_A} \sum_{k=1}^{n_B} \left( \underline{A}_j^{k'} \underline{B}_j^k \prod_{\substack{h=1 \\ h \neq j}}^{N} \underline{R}_h^T \underline{A}_h^{k'} \underline{A}_h^k \underline{B}_h^k \right) \tag{89}
$$

This is then the equation determining $R_j$. More details about the
Proper Generalized Decomposition (PGD) method sketched above
can be found in [35-38].

References

[1] L. Onsager, The effect of shape on the interaction of colloidal particles, Ann. N.Y. Acad. Sci. 51 (1949) 627-659.

[2] J.F. Brady, R.J. Phillips, R.J. Lester, G. Bossis, Dynamic simulation of hydrodynamically interacting suspensions, J. Fluid Mech. 195 (1988) 257-280.

[3] M.B. Mackaplow, E.S.G. Shaqfeh, Numerical study of rigid non-Brownian fibers, J. Fluid Mech. 329 (1969) 155-186.

[4] G.B. Jeffery, The motion of ellipsoidal particles immersed in a viscous fluid, Proc. Roy. Soc. Lond. Ser. A 102 (1922) 161-179.

[5] D.L. Koch, h E.S.G. Shaqfe, The average rotation rate of a fiber in the linear flow of a semidilute suspension, Phys. Fluids A 2 (1990) 2093-2102.

[6] M. Rahnama, D.L. Koch, Y. Iso, C. Cohen, Hydrodynamic, translational diffusion in fiber suspensions subject to simple shear flow, Phys. Fluids A 5 (1992) 849-862.

[7] C.A. Stover, The Dynamics of Fibers Suspended in Shear Flows, PhD Thesis, School of Chemical Engineering, Cornell University, 1991.

[8] F. Du, R.C. Scogna, W. Zhou, S. Brand, J.E. Fischer, K. Winey, Nanotube networks in polymer nanocomposites: rheology and electrical conductivity, Macromolecules 37 (2004) 9048-9055.

[9] J. Xu, S. Chaterjee, K.W. Koelling, Y. Wang, S.E. Bechtel, Shear and extensional rheology of carbon nanofiber suspensions, Rheol. Acta 44 (2005) 537-562.

[10] Q. Zhang, I.A. Archer, Polyethylene oxide/silica nanocomposites: structure and rheology, Langmuir 18 (2002) 10435-10442.

[11] A. Ma, M. Mackley, F. Chinesta, The microstructure and rheology of carbon nanotube suspensions, Int. J. Mater. Form. 2 (2008) 75-81.

[12] A. Ma, F. Chinesta, A. Ammar, M. Mackley, Rheological modeling of carbon nanotube aggregate suspensions, J. Rheol. 52 (2008) 1311-1330.

[13] A. Clebsch, Über die Integration der hydrodynamische Gleichungen, J. Reine Angew. Math. 56 (1895) 1-10.

[14] V.I. Arnold, Sur la géométrie différentielle des groupes de Lie de dimension infini et ses applications dans l'hydrodynamique des fluides parfaits, Ann. Inst. Fourier 16 (1966) 319.

[15] J.E. Marsden, A. Weinstein, Coadjoint orbits, vortices and Clebsch variables for incompressible fluids, Physica D 7 (1983) 305-323.

[16] M. Grmela, H.C. Öttinger, Dynamics and thermodynamics of complex fluids: general formulation, Phys. Rev. E 56 (1997) 6620-6633.

[17] M. Grmela, Öttinger, Dynamics and thermodynamics of complex fluids: illustration of the general formalism, Phys. Rev. E 56 (1997) 6633-6650.

[18] I.E. Dzyaloshinskii, G.E. Volovick, Poisson brackets in condense matter physics, Ann. Phys. (NY) 125 (1980) 67-97.

[19] M. Grmela, Particle and bracket formulations of kinetic equations, Contemporary Math. 28 (1984) 125-132. Physics Letters A 102, 355 (1984).

[20] A.N. Kaufman, Dissipative Hamiltonian systems: a unifying principle, Phys. Lett. A 100 (1984) 419.

[21] P.J. Morrison, Bracket formulation for irreversible classical fields, Phys. Lett. A 100 (1984) 423.

[22] M. Grmela, Bracket formulation of diffusion-convection equations, Physica D 21 (1986) 179-212.

[23] A.N. Beris, B.J. Edwards, Thermodynamics of Flowing Systems, Oxford Univ. Press, Oxford, 1994.

[24] H.C. Öttinger, Beyond Equilibrium Thermodynamics, Wiley Interscience, 2005.

[25] M. Grmela, Multiscale equilibrium and nonequilibrium thermodynamics in chemical engineering, Adv. Chem. Eng. 39 (2010) 76-128.

[26] M. Grmela, Why GENERIC, J. NonNewton. Fluid Mech. 165 (2010) 980-986.

[27] S.R. de Groot, P. Mazur, Non-Equilibrium Thermodynamics, Dover, New York, 1984.
[28] D.D. Holm, J.E. Marsden, T. Ratiu, A. Weinstein, Nonlinear stability of fluid and plasma equilibria, Phys. Rep. 123 (1985) 1.

[29] J.E. Marsden, T.S. Ratiu, Introduction to Mechanics and Symmetry, Texts in Applied Mathematics, vol. 17, Springer, New York, 1999.

[30] M. Grmela, A. Ammar, F. Chinesta, Extra stress tensor in fiber suspensions: mechanics and thermodynamics, J. Rheol. 55 (2011) 17–42.

[31] R.F. Kayser, H.J. Raveché, Bifurcation in Onsager's model of the isotropic-nematic transition, Phys. Rev. A 17 (1978) 2067–2071.

[32] W. Maier, A.Z. Saupe, Naturforschung 15a (1960) 287.

[33] D. Ruelle, Thermodynamic Formalism: The Mathematical Structure of Equilibrium Statistical Mechanics, Cambridge Mathematical Library, 2004.

[34] M.S. Green, Generalized Ornstein–Zernike approach to critical phenomena, J. Math. Phys. 9 (1968) 875.

[35] A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, J. Non-Newton. Fluid Mech. 139 (2006) 153–176.

[36] A. Ammar, M. Normandin, F. Chinesta, Solving parametric complex fluids models in rheometric flows, J. Non-Newton. Fluid Mech. 165 (2010) 1588–1601.

[37] F. Chinesta, A. Ammar, E. Cueto, Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models, Arch. Comp. Methods Eng. 17 (2010) 327–350.

[38] F. Chinesta, A. Ammar, A. Leygue, R. Keunings, An overview of the proper generalized decomposition with applications in computational rheology, J. Non-Newton. Fluid Mech. 166 (2011) 578–592.

[39] V. Zmievski, M. Grmela, M. Bousmina, S. Dagréou, Nonlinear microstructure and rheology of semidilute colloidal suspensions of structurless partiales, Phys. Rev. E 71 (2005) 051503.

[40] [40] G. Maitrejean, A. Ammar, F. Chinesta, M. Grmela, Deterministic solution of the kinetic theory model of colloidal suspensions of structurless particles, Rheol. Acta 51 (2012) 527–543.

[41] M. Grmela, G. Maitrejean, A. Ammar, F. Chinesta, Kinetic theory of colloidal suspensions: rheology and migration, Rheol. Acta, submitted for publication.