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

HAL
open science

# Extra stress tensor in fiber suspensions: Mechanics and thermodynamics

Miroslav Grmela, Amine Ammar, Francisco Chinesta

► To cite this version:

Miroslav Grmela, Amine Ammar, Francisco Chinesta. Extra stress tensor in fiber suspensions: Mechanics and thermodynamics. Journal of Rheology, 2011, 55, pp.17-42. 10.1122/1.3523538. hal-01516400

HAL Id: hal-01516400

https://hal.science/hal-01516400

Submitted on 30 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
# Extra stress tensor in fiber suspensions: mechanics and thermodynamics

Miroslav Grmela¹*, Amine Ammar², Francisco Chinesta³

¹ Ecole Polytechnique de Montreal, C.P.6079 suc. Centre-ville, Montreal, H3C 3A7, Quebec, Canada

² Laboratoire de Rheologie, INPG, UJF, CNRS (UMR 5520), 1301 rue de la piscine, BP 53 Domaine Universitaire, F-38041 Grenoble Cedex 9, France

³ EADS Corporate International Chair GeM, Ecole Centrale de Nantes BP 92101, 44321 Nantes, cedex 3, France

September 8, 2010

## Abstract

Formulas expressing the extra stress tensor, $\sigma$, in fiber suspensions in terms of microstructural state variables are derived by using two types of arguments: mechanical and thermodynamical. Results are compared for the distribution function ($\psi$) and the orientation tensor, $A$, playing the role of state variables. The main results are the following: (i) In the thermodynamical analysis the formulas arise as compatibility conditions between the time evolution of the fluid velocity and the time evolution of the internal structure. (ii) A complete agreement among the formulas arising in mechanics and thermodynamics is seen only in kinetic theory (i.e. with $\psi$ as the state variable) and only with the Chan-Terentjev mechanical formula. (iii) Theoretical arguments as well as numerical illustrations indicate that the larger is the role of the reversible part of the time evolution of the microstructure the larger is the difference in predicted stresses (i.e. the formulas for $\sigma$ evaluated at solutions of the microstructural equations) calculated with the thermodynamic and the Dinh-Armstrong mechanical formulas.

## 1 Introduction

Let our interest in complex fluids be mainly directed towards their macroscopic flow behavior. The setting for their theoretical investigations has to include

<footnote>*corresponding author: e-mail: miroslav.grmela@polymtl.ca</footnote>
therefore classical hydrodynamics as its one (macroscopic) component. The second (microscopic) component is needed because the complex fluids involve a microstructure (e.g. various suspended particles, macromolecules or membranes) that changes in time on the scale that is comparable to the time scale of changes of hydrodynamic fields. The time evolution of the microstructure cannot be thus ignored even if our main interest is focused on the macroscopic flow.

In this paper we limit ourselves to isothermal and incompressible fluids. The overall fluid behavior is thus described only by the overall velocity fields $v(r)$. By $r$ we denote the position vector. The fluid density $\rho$ and the temperature $T$ are constants. Because of the convenience in both physical interpretations and mathematical formulations, we shall use systematically the momentum field $u(r)$ instead of the velocity field. If the only term in the energy that depends on the velocity is the kinetic energy $\int dr \frac{1}{2} \rho v^2$ then the relation between $u$ and $v$ is simply $u(r) = \rho v(r)$. The state variable characterizing the microstructure will be denoted by the symbol $\xi(r)$. The field $\xi$ can depend also on other variables (as for example $\xi$ chosen in Section 2). The complete set of state variables used in this paper (denoted formally by the symbol $x$) is thus

$$x = (u, \xi) \qquad (1)$$

Now we turn to the equations governing the time evolution of (1). As for the momentum field $u(r)$, the equation governing its time evolution is the local conservation law

$$\frac{\partial u}{\partial t} = -\nabla \cdot \left( \frac{uu}{\rho} \right) - \nabla p - \nabla \cdot \boldsymbol{\sigma} \qquad (2)$$

where $p$ is the hydrostatic scalar pressure (determined by the incompressibility requirement), and $\boldsymbol{\sigma}$ is the extra stress tensor. We shall consider in this paper only symmetric extra stress tensors.

The complex fluids under investigation in this paper will be suspensions of fibers. The microstructure characterized by $\xi(r)$ is thus a characteristics of the distribution of suspended fibers. Let the microscopic component of the setting of theoretical rheology be formally represented by the equation

$$\frac{\partial \xi}{\partial t} = \chi(u, \xi) + \omega(u, \xi) \qquad (3)$$

Advection of $\xi$ with the flow is expressed in the term $\chi(u, \xi)$. The physics behind the advection is the Stokes problem describing interaction of a fluid with fibers suspended in it. The term $\chi$ can also be seen as a term in which the external forcing on the microstructure (by imposing a flow) is expressed. The dependence on $u$ is such that $\chi(-u, \xi) = -\chi(u, \xi)$. The second term $\omega(u, \xi)$ represents dissipation. Most often $\omega$ is independent of $u$. If it depends on $u$, the dependence is such that $\omega(-u, \xi) = \omega(u, \xi)$. The physics that is behind the term $\omega(\xi)$ is best seen in the setting of thermodynamics. We shall discuss it below in Section 2.2.
Equations (2) and (3) are coupled in $\sigma(\xi)$ and $\chi(\mathbf{u}, \xi)$. The main objective of this paper is to compare investigations, made by following two routes (mechanical and thermodynamical), of the relation between $\sigma(\xi)$ and $\chi(\mathbf{u}, \xi)$ in the context of fiber suspensions with two choices of $\xi$ ($\xi$ = one fiber distribution function, and $\xi$ = an orientation tensor). The thermodynamical route is moreover divided into two: classical nonequilibrium thermodynamics and GENERIC.

# 2 Kinetic theory

The internal structure (i.e. the distribution of suspended fibers) is chosen in this section to be characterized by one fiber configuration space distribution function $\psi(\mathbf{r}, \mathbf{p})$, where $\mathbf{p}$ is the unit vector along the fiber.

## 2.1 Mechanics

Equation (2) has two mechanical interpretations. First, it is a local conservation law for the momentum $\mathbf{u}(\mathbf{r})$ (i.e. the right hand side of (2) is divergence of a flux), and second, it is a continuum version of Newton's law. From the second interpretation we then see that $\sigma$ is a force acting on surface. We can thus find $\sigma$ by identifying the local surface forces acting on fibers. This is indeed the way the expression for $\sigma$ has originally been derived. The mechanical approach to calculating $\sigma$ in fiber suspensions has been initiated by Batchelor (1970), (1971) and later continued in Evans (1975), Dinh and Armstrong (1984) and Lipscomb et al. (1988). The Dinh-Armstrong formula is a special case of the Lipscomb et al. formula. The most recent and the most complete analysis of the forces inside fiber suspensions has been made in Chan and Terentjev (2007). The expressions derived by Chan and Terentjev will be denoted $\sigma^{(mechCT)}(\psi)$ and the expression derived by Dinh and Armstrong $\sigma^{(mechDA)}(\psi)$. The upper indices $(mechCT)$ and $(mechDA)$ denote the provenance of the formulas.

The advection term $\chi(\psi)$ used in Chan and Terentjev (2007) is

$$ \chi(\mathbf{u}, \psi) = -\frac{\partial}{\partial r_j} \left( \psi \frac{u_j}{\rho} \right) - \frac{\partial}{\partial p_j} (\dot{p}_j \psi) \quad (4) $$

where

$$ \dot{p}_i = \Omega_{jk} p_k + \lambda D_{jk} p_k - \lambda D_{jm} p_j p_l p_m, \quad (5) $$

The second term on the right hand side (3) takes the form

$$ \omega(\psi) = \frac{\partial}{\partial p_j} \left( \Lambda (\delta_{jk} - p_j p_k) \frac{\partial \psi}{\partial p_k} \right) \quad (6) $$

We use hereafter the summation convention $i, j, k = 1, 2, 3$, the tensors $D$ and $\Omega$ are respectively the symmetric and antisymmetric parts of the velocity gradient $\nabla \mathbf{u}$, $\lambda = \frac{a^2-1}{a^2+1}$; $a$ is the fiber aspect ratio (fiber length to fiber diameter ratio).
The formula for the stress tensor $\sigma^{(mechCT)}$ derived in Chan and Terentjev (2007) is the following:

$$ \sigma_{ij}^{(mechCT)} = -\lambda \frac{1}{2} \left( \int d\mathbf{p} \psi p_i \frac{\partial}{\partial p_j} (\Phi_\psi) + \int d\mathbf{p} \psi p_j \frac{\partial}{\partial p_i} (\Phi_\psi) - 2 \int d\mathbf{p} \psi p_i p_j p_k \frac{\partial}{\partial p_k} (\Phi_\psi) \right) \quad (7) $$

where

$$ \Phi = \int dr \frac{u^2}{2\rho} + k_B T n_{fib} \int dr \int d\mathbf{p} \, \psi \ln \psi, \quad (8) $$

$n_{fib}$ is the number density of the fibers, $k_B$ is the Boltzmann constant, $T$ is the temperature. We use hereafter the shorthand notation: $\frac{\partial\Phi}{\partial\mathbf{u}(\mathbf{r})} = \Phi\mathbf{u}(\mathbf{r})$ and $\frac{\partial\Phi}{\partial\xi} = \Phi_\xi$. Moreover, we also use the symbol $\partial$ for both the ordinary partial derivative and the Volterra functional derivative. For example, $\frac{\partial\Phi}{\partial\mathbf{u}(\mathbf{r})}$ is the Volterra functional derivative since $\Phi$ is a real valued function of a function $\mathbf{u}(\mathbf{r})$ and $\frac{\partial\Phi}{\partial\mathbf{u}_i}$ is an ordinary partial derivative of $\Phi$ that is a real valued function of a finite number $N$ of independent variables ($\mathbf{u}_1, ..., \mathbf{u}_N$) where $\mathbf{u}_i = \mathbf{u}(\mathbf{r}_i)$; $i = 1, 2, ..., N$ and ($\mathbf{r}_1, ..., \mathbf{r}_N$) is a discretization of $\mathbf{r} \in \mathbb{R}^3$.

If we replace in (7) $\Phi_\psi$ with $k_B T n_{fib}(\ln \psi + 1)$ then indeed (7) is the symmetric part of the stress tensor appearing in Eq.(64) in Chan and Terentjev (2007). We are choosing to write the formula in the form (7) because it appears in this way in thermodynamics (in the following Section 2.2). We shall see that the quantity $\Phi$ has the physical interpretation of free energy.

The formula for the fiber contribution to the extra stress tensor derived by Dinh and Armstrong (1984) is the following:

$$ \sigma_{ij}^{(mechDA)} = 2\eta_m N_p A_{ijkl} D_{kl} \quad (9) $$

where $A_{ijkl} = \int d\mathbf{p} \, p_i p_j p_k p_l \psi(\mathbf{r}, \mathbf{p})$, $\eta_m$ is the matrix viscosity, and $N_p$ is a dimensionless parameter called the particle number that represents the relative importance of the fibers, $D$ is the symmetric velocity gradient.

We shall compare the formulas (7), (9) and also other formulas derived below in Sections 2.2 and 2.3 in Section 2.4.

## 2.2 Thermodynamics

In this subsection we consider externally unforced fiber suspensions. The experimental observation on which we shall concentrate is the approach to equilibrium states at which the behavior of suspensions is seen to be well described by classical equilibrium thermodynamics. We look for a structure of the equations governing the time evolution of the state variables (1) that will guarantee that solutions of the equations agree with the above experimental observation.

We begin with a short preparation. We say that Eqs.(2), (3) are time reversible if the inversion of the sign of the time $t$ is compensated by the transformation $(\mathbf{u}, \xi) \rightarrow (-\mathbf{u}, \xi)$ (i.e. a simultaneous application of $t \rightarrow -t$ and
($u, \xi) \to (-u, \xi)$ leaves (2), (3) invariant). If the simultaneous application of $t \to -t$ and $(u, \xi) \to (-u, \xi)$ does not leave (2), (3) invariant, the two equations are called time irreversible.

First, we turn to Eq.(2). We split $\sigma$ into two parts:

$$ \sigma = \sigma^{(+)} + \sigma^{(-)} \quad (10) $$

where $\sigma^{(+)}$ is invariant with respect to $(u, \xi) \to (-u, \xi)$, and $\sigma^{(-)}$ changes its sign if $(u, \xi) \to (-u, \xi)$ is applied. We see now that Eq.(2) with $\sigma = \sigma^{(+)}$ is time reversible.

Next, we turn to Eq.(3). Since $\omega$ is assumed to depend on $u$ in such a way that $\omega$ remains invariant with respect to the transformation $u \to -u$, Eq.(3) without $\chi$ is time irreversible and Eq.(3) without $\omega$ is time reversible (since we assume that $\chi(-u, \xi) = -\chi(u, \xi)$).

Now, we are in position to present the thermodynamic argument that will lead us to an expression for the extra stress tensor $\sigma$. The argument is the same for both choices of the internal state variables (i.e. for both $\psi(r, p)$ and $A(r)$). We present it therefore below with $\xi(r, p)$ representing both choices. Let $\Phi(u, \xi)$ be the free energy. We choose it in such a way that $\Phi(u, \xi) = \Phi(-u, \xi)$. We recall its role in the time evolution. As a consequence of the observed approach of externally unforced fluids to equilibrium states (denoted $(u, \xi)_{eq}$), the following inequality holds

$$ \frac{d\Phi}{dt} = \int dr \left[ \Phi_{u_j} \frac{\partial u_j}{\partial t} + \int dp \, \Phi_{\xi} \frac{\partial \xi}{\partial t} \right] \le 0 \quad (11) $$

The free energy $\Phi$ plays thus the role of the Lyapunov function for the approach. The equilibrium states $(u, \xi)_{eq}$, reached as $t \to \infty$, are the state at which $\Phi$ reaches its minimum and are therefore solutions to $\Phi u = 0, \Phi \xi = 0$) [we use hereafter the shorthand notation: $\frac{\partial \Phi}{\partial u} = \Phi u$ and $\frac{\partial \Phi}{\partial \xi} = \Phi \xi$] Moreover, $\Phi$ evaluated at the equilibrium state $(u, \xi)_{eq}$ becomes the equilibrium thermodynamic free energy determining the equilibrium thermodynamic properties of the fluid under consideration.

We shall replace (11) by a somewhat stronger statement:

$$ \left(\frac{d\Phi}{dt}\right)^{(-)} = \int dr \left[ \Phi_{u_i} \frac{\partial(u_i v_j + p \delta_{ij} + \sigma_{ij}^{(+)})}{\partial r_j} + \int dp \, \Phi_{\xi} \chi \right] = 0 \quad (12) $$

$$ \left(\frac{d\Phi}{dt}\right)^{(+)} = \int dr \left[ -\Phi_{u_i} \frac{\partial \sigma_{ij}^{(-)}}{\partial r_j} + \int dp \, \Phi_{\xi} \omega \right] < 0 \quad (13) $$

It is clear that Eqs.(12) and (13) imply (11) (since $\frac{d\Phi}{dt} = (\frac{d\Phi}{dt})^{(+)} + (\frac{d\Phi}{dt})^{(-)}$) but (11) does not, in general, imply (12),(13).

The stress tensors arising from Eqs.(12) and (13) will be hereafter denoted $\sigma^{(th)}$ in order to indicate their provenance.
It is easy to verify (see Grmela (1985)) that Eq.(12) is satisfied if

$$ \sigma_{ij}^{(th+)} = - \int dp \, \Phi_{\xi} \frac{\partial \chi}{\partial \left( \frac{\partial \Phi_{u_i}}{\partial r_j} \right)} \quad (14) $$

Since $\omega$ is independent of $u$, the inequality (13) implies:

$$ \int dr \, \Phi_{u_i} \frac{\partial \sigma_{ij}^{(-)}}{\partial r_j} = - \int dr \, \frac{\partial \Phi_{u_i}}{\partial r_j} \sigma_{ij}^{(-)} > 0 \quad (15) $$

and

$$ \int dr \int dp \, \Phi_{\xi} \omega < 0 \quad (16) $$

The first equality in (15) is a result of by part integration in which the term involving integral over boundary equals zero due to boundary conditions. The inequality (15) restricts the freedom of choice of $\sigma^{(-)}$ but, unlike the equality (12) that determines $\sigma^{(+)}$, it does not determine it. We note that, in particular,

$$ \sigma_{ij}^{(th-)} = -C_{ik}(\xi)D_{kj}; \quad C(\xi) \text{ is a positive definite tensor,} \quad (17) $$

($D$ is the symmetric velocity gradient) will always satisfy (15) independently of the particular choice of the tensor $C(\xi)$.

In conclusion, the thermodynamic argument that we have just recalled provides an expression (14) for $\sigma^{(th+)}$ but leaves $\sigma^{(th-)}$ undetermined (only constrained by the inequality (15)). In addition, we have also arrived at the inequality (16) that we shall exploit later.

Now, we apply (14)-(16) to Eq.(3) with (4) and (6). From (14) we get

$$ \sigma_{ij}^{(th+)} = - \int dp \, \psi \frac{\partial \Phi_{\psi}}{\partial p_k} \frac{\partial \dot{p}_k}{\partial D_{ij}} \quad (18) $$

which equals to the right hand side of (7).

So far, we have not yet identified the free energy $\Phi$. We have only used $\Phi u = \frac{u}{\rho} = v$ which is the velocity field. We note that if $\Phi = \frac{u^2}{2\rho}$, i.e. the kinetic energy then indeed $\Phi u$ is the velocity field. But what is the remaining part of the free energy? It seems that this question does not need to be answered in the mechanical argument. However, it does have to be answered if we want to be certain that the time evolution equation (3) with (4) and (6) are physically sound. The role of (6) is to bring (in the absence of external influence) solutions to (3) to equilibrium states. It is exactly the inequality (16) that gives $\omega(\psi)$ its physical interpretation. We shall now use it to determine the free energy $\Phi(u, \psi)$.

We proceed as follows. First we propose the free energy (8) and then we verify that the inequality (16) holds for $\omega$ given in (6).

We note that if we succeed to cast $\omega$ into the form $-\frac{\partial \Xi}{\partial \Phi_{\psi}}$ then indeed the inequality (16) holds provided $\Xi$, a real valued function of $\Phi_{\psi}$ that we shall call
a dissipation potential, is such that $\Xi(0) = 0$, $\Xi$ reaches its minimum at 0, and $\Xi$ is convex in a neighborhood of 0 [recall that if $\Xi$ is a real valued function of $x \in \mathbb{R}$ satisfying the properties listed above than indeed $x \frac{d\Xi}{dx} > 0$].

Now we identify the dissipation potential $\Xi$ (see Grmela (2008)):

$$ \Xi(\Phi_{\psi}) = \int dr \int dp (R\Phi_{\psi})_j \frac{1}{2} \psi \Lambda(R\Phi_{\psi})_j, \quad (19) $$

where $\Lambda > 0$ is a coefficient and $(R\Phi_{\psi})_j = (p \times \frac{\partial}{\partial p} \Phi_{\psi})_j; \times$ denotes the vector product.

We have thus proven that (8) is the free energy (implicitly) present in the kinetic equation (3), (4), (6), and (7). The first term in (8) is the kinetic energy and the second term is the entropy multiplied by $-1$.

We end this subsection with three observations.

**Observation 1** According to the virtual work principle (Doi (1983), Doi and Edwards (1986) p.70, see also Feng et al. (2000), Wang (2002), Sircar and Wang (2009)) the variation of the microstructural contribution to the free energy density (i.e. the variation of the second term on the right hand side of (8)) equals the work done to the material by the elastic stress (i.e. in our notation $\sigma^{(+)}$). We note now that this is in fact an alternative physical interpretation of the equality (12). Indeed, the first term on the right hand side of (12) (after making by parts integration and using the expression $p = \phi - u_k \Phi_{u_k} - \int dp \psi \Phi_{\psi}$, where $\phi$ defined by $\Phi = \int dr \phi$ for the hydrostatic pressure $p$) becomes the work done to the material by the elastic stress (i.e. $\sigma^{(+)}$) multiplied by $-1$. The second term is clearly the variation of the microstructure contribution to the free energy density. We thus conclude that the virtual work principle is equivalent to the requirement (expressed in (12)) that the total free energy remains unchanged during the reversible time evolution.

**Observation 2** We have seen that the thermodynamic argument (i.e. the requirement that solutions of the time evolution equations agree with the experimental observation constituting the basis of classical equilibrium thermodynamics) does leave the irreversible part $\sigma^{(-)}(\psi, D)$ of the extra stress tensor undetermined. We have seen that, for example, any expression of the form (17) is admissible. We shall now argue that $\sigma^{(-)}(\psi, D) = 0$. The physics expressed in $\sigma^{(-)}(\psi, D)$ is the physics of dissipative processes taking place in suspensions. But we have already taken such processes into account in the term $\omega(\psi)$ in Eq.(3) governing the time evolution of the microstructural state variable $\psi(r, p)$. Assuming that the dissipative processes are expressed completely in $\omega(\psi)$, we have no physical basis (inside the thermodynamic argument) for constructing $\sigma^{(-)}(\psi, D)$. We thus conclude that (inside the thermodynamic argument) $\sigma^{(-)}(\psi, D) = 0$. We shall discuss this point further in Section 2.4.

**Observation 3** The split of the time evolution into the reversible and irreversible parts depends on the choice of the morphological state variable $\xi$ in (1). In other words, the split of the time evolution into the reversible and irreversible parts is level dependent. By a level of description we mean the way we regard the system under consideration. For example, the choice of the morphological
state variable that we have made in this section corresponds to a more microscopic (i.e. involving more details) level than the level corresponding the choice that we shall make in Section 3. We shall illustrate the level dependence of the reversible-irreversible splitting of the time evolution on two examples. Let the complex fluid under investigation be composed, for the sake of simplicity of the illustration, of only point particles. In the first example we choose $\xi$ to consist of the position coordinates and velocities of all particles composing the fluid. The time evolution in this case is completely time reversible, there is no time irreversible part. In particular, all the interactions among the particles enter the time reversible part. In the second example we choose the morphological state variable to be the one particle distribution function, either in the configuration space (i.e. $\xi = \psi(\mathbf{r})$) or in the phase space (i.e. $\xi = \psi(\mathbf{r}, \mathbf{v})$, where $\mathbf{v}$ is the velocity). Let the particles interact among each other and let the interaction be expressed in a mean-field type potential $V^{(mf)}(\psi)$. In the case of $\xi = \psi(\mathbf{r}, \mathbf{v})$, this potential enters into the time reversible part (in the form of $\frac{\partial}{\partial \mathbf{v}} \left( \psi \frac{\partial V^{(mf)}(\psi)}{\partial \mathbf{r}} \right)$) as well as into the time irreversible part. On the other hand, in the case of $\xi = \psi(\mathbf{r})$, the mean-field potential is absent in the time reversible part.

## 2.3 GENERIC

As in the previous subsection, we consider externally unforced fiber suspensions. Our objective is to identify the structure of the time evolution equations guaranteeing that their solutions agree with the experimentally observed approach to equilibrium. The structure identified in this subsection will be richer that the one identified in the previous section. In addition to requiring the compatibility with thermodynamics we shall require that the mechanics (Newton's law) does not enter only Eq.(2) but also Eq.(3). In this subsection we shall combine mechanics and thermodynamics into a single formalism addressing both Eqs.(2) and (3).

To begin with, we have to choose a formulation of mechanics that is appropriate for making the combination. Among several possible mathematical formulations of Newton's law, we shall choose the Hamiltonian formulation. It has been first used in the context of hydrodynamics by Clebsch (1895). It appears to be the most convenient for expressing the advection term $\chi(\mathbf{u}, \xi)$ in (3) as a result of the mechanical interaction of the fluid with an obstacle and for combining mechanics with thermodynamics.

Clebsch (1895) and Arnold (1966) have shown that the Euler equation, i.e. Eq.(2) with $\sigma = 0$, can be cast into the form

$$ \frac{dA}{dt} = \{A, E\}^{(u)} \quad \text{holds for all } A \qquad (20) $$

where $A$ is a sufficiently regular real valued function of $\mathbf{u}$ and $E$ is the energy
$\int dr \frac{u^2}{2\rho}$ (i.e. (8) with the second term missing) and

$$ \{A, B\}^{(u)} = \int dr u_i \left( \frac{\partial A_{u_i}}{\partial r_j} B_{u_j} - \frac{\partial B_{u_i}}{\partial r_j} A_{u_j} \right) \quad (21) $$

is a Poisson bracket in which the kinematics of continuum (i.e. the Lie group of transformations $\mathbb{R}^3 \to \mathbb{R}^3$) is expressed.

Now, we supplement the momentum field $u(r)$ with the fiber distribution function $\psi(r, p)$ and ask the question of what is the Poisson bracket expressing kinematics of $(u(r), \psi(r, p))$. If we assume that $\psi$ is passively advected (Lie dragged) by $u$ then the bracket is given by (see Grmela (1988))

$$ \{A, B\}^{(u, \psi)} = \{A, B\}^{(u)} + \int dr \int dp \, \psi \left( \frac{\partial A_{\psi}}{\partial r_j} B_{u_j} - \frac{\partial B_{\psi}}{\partial r_j} A_{u_j} \right) \\ + \lambda \int dr \int dp \, \psi n_j (\delta_{kl} - p_k p_l) \left( \frac{\partial A_{\psi}}{\partial p_l} \frac{\partial B_{u_k}}{\partial r_j} - \frac{\partial B_{\psi}}{\partial p_l} \frac{\partial A_{u_k}}{\partial r_j} \right) \quad (22) $$

The mechanical content of the time reversible part of the evolution of $(u(r), \psi(r, p))$ is now expressed by requiring that the governing equations have still the form (20) (expressing the mechanical content of the Euler equation) but with $\{A, B\}^{(u)}$ replaced by $\{A, B\}^{(u, \psi)}$ and $E$ replaced by the total energy. The explicit form of the governing equations (including the explicit formula for the extra stress tensor) emerges now in the following calculations. The left hand side of (20) is written as $\int dr A_{u_i(r)} \frac{\partial u_i(r)}{\partial t} + \int dr \int dp A_{\psi(r,p)} \frac{\partial \psi(r, p)}{\partial t}$. The right hand side of (20) is written (with the use of by parts integration) in the form $\int dr A_{u_i(r)}(\bullet)_i + \int dr \int dp A_{\psi(r,p)}(\bullet, \bullet)$, where $(\bullet)_i$ and $(\bullet, \bullet)$ represent the expressions obtained in the calculations. Since Eq.(20) is required to hold for all $A$, we obtain $\frac{\partial u_i(r)}{\partial t} = (\bullet)_i$ and $\frac{\partial \psi(r, p)}{\partial t} = (\bullet, \bullet)$. These equations are the same as Eqs.(2),(3) with $\chi(u, \psi)$ given in (4), the term $\omega(\psi)$ missing, and $\sigma$ given in (7). The extra stress tensor $\sigma$ arising in this analysis will be hereafter denoted by the symbol $\sigma^{(GENERIC)}$ to denote its provenance. We have thus shown that $\sigma^{(th+)} = \sigma^{(GENERIC+)}$. This result is not surprising if we realize that (20) implies $\frac{d\Phi}{dt} = \{\Phi, \Phi\} = 0$ and that this equality served as a basis (see (12)) on which (7) was obtained in Section 2.2.

Having addressed the time reversible part of the evolution, we now proceed to the time irreversible part. Here we join the previous subsection. The only requirement is the compatibility with thermodynamics expressed by the inequality $d\Phi/dt \le 0$ required to hold during the time evolution. As in the previous subsection, we note that if

$$ \left( \frac{\partial \phi}{\partial t} \right)_{irr} = - \frac{\partial \Xi}{\partial \Phi_{\psi}} \quad (23) $$

$$ \left( \frac{\partial \mathbf{u}}{\partial t} \right)_{irr} = - \frac{\partial \Xi}{\partial \Phi_{\mathbf{u}}} \quad (24) $$
are the equations governing the time irreversible evolution then indeed $d\Phi/dt \le 0$ holds provided $\Xi(\Phi_\psi, \Phi_u)$ is a potential, called a dissipation potential, satisfying: (i) $\Xi(0, 0) = 0$, (ii) $\Xi$ reaches its minimum at $(0, 0)$, and (iii) $\Xi$ is convex in a neighborhood od $(0, 0)$. As we have seen at the end of Section 2.2, the choice (19) of the dissipation potential corresponds to the choice (6) of $\omega$ (see Eq.(3)).

The potential (19) is not however the only dissipation potential for which $-\frac{\partial \Xi}{\partial \Phi_\psi} = \frac{\partial}{\partial p_j} \left( \Lambda(\delta_{jk} - p_j p_k) \frac{\partial \psi}{\partial p_k} \right)$ and $d\Phi/dt \le 0$. For example, if we choose

$$ \Xi(\Phi_\psi, \Phi_u) = \int dr \int dp (R\Phi_\psi)_j \frac{1}{2} \psi \Lambda(R\Phi_\psi)_j + \frac{1}{2} \int dr D_{ij} C_{jk} D_{ki} \quad (25) $$

$(D_{ij} = 1/2 (\partial_i \Phi_{u_j} + \partial_j \Phi_{u_i})$, $\partial_i = \partial/\partial r_i$, and $C$ is a positive definite tensor) then we still satisfy the inequality $d\Phi/dt \le 0$ and the time evolution of $(u(r), \psi(r, p)$ is governed by Eq.(2) (with $\sigma = \sigma^{(+)} + \sigma^{(-)}$, where $\sigma^{(+)}$ is given in (7) and $\sigma^{(-)}$ in (17)) and Eq.(3) with $\chi$ and $\omega$ given in (4),(6). We thus conclude that the stronger requirement that we have used in this subsection still leaves the part $\sigma^{(-)}$ undetermined.

We end this subsection with three observations.

**Observation 1** Exactly the same argument as the one introduced in *Observation 2* in the previous section leads us again to $\sigma^{(-)} = 0$. Section 2.4 below provides additional arguments in favor of this conclusion.

**Observation 2** If we combine the Hamiltonian reversible time evolution presented above with the dissipative time evolution generated by $\omega(\psi)$ (expressed with the help of the dissipation potential discussed in Section 2.2) in such a way that the inequality (11) is guaranteed, we arrive at a time evolution generated by equations that have started to appear in Dzyaloshinski and Volovick (1980) and later in Grmela (1984), Kaufman (1984), Morrison (1984), Beris and Edwards (1994). An abstract equation representing an appropriately combination of the reversible Hamiltonian and the irreversible dissipative vector field has been called GENERIC in Grmela and Ottinger (1997) and Ottinger and Grmela (1997). GENERIC has been then further developed in Grmela (2001), (2002), (2004), (2010) and in a different direction in Ottinger (1998),(2005). The compatibility of the reversible and irreversible time evolution is essentially a requirement guaranteeing satisfaction of (12) and (13). Beside this constraint, the irreversible part of the time evolution remains undetermined. This then means that, like in the preceding section, the irreversible part $\sigma^{(GENERIC-)}$ of the extra stress tensor remains undetermined (except that (13) is required to hold). The geometrical interpretation in which the time evolution generated by GENERIC becomes a continuous sequence of Legendre transformations (see Grmela (2010)) provides an additional insight into the coupling of the reversible and irreversible time evolutions.

**Observation 3** The third observation is about the bracket (22). It can easily be verified that (22) is a Poisson bracket only for $\lambda = 1$ and $\lambda = -1$ with the tensor $\partial B_u / \partial r$ transposed. Otherwise, the Jacobi identity is not verified. We can understand this observation as follows: From the mathematical point of view, it is well known (see the concept of the semi direct product in
e.g. Marsden and Ratiu (1999)) that the scalar $\psi(\mathbf{r}, \mathbf{p})$ has to be passively advected (in other words, Lie dragged) in order that the extended space with coordinates $(\mathbf{u}(\mathbf{r}), \psi(\mathbf{r}, \mathbf{p}))$ be equipped with a Poisson bracket extending the Poisson bracket (21). The extended Poisson bracket is then the bracket (22). From the physical point of view, the advection (in our case Eq. (5)) is a result, as we have already mentioned in the text following Eq. (3), of solving the Stokes problem. The passive advection that is expressed in (22) corresponds to a special solution to the Stokes problem, namely to the solution in which the flow of the fluid remains completely undisturbed by the presence of suspended fibers. This special solution can be considered to be an acceptable approximation of the actual solution only for fibers with special shapes like an infinitely thin fiber and an infinitely thin plate. Indeed, if we recall that $\lambda = \frac{a^2-1}{a^2+1}$ where $a$ is the fiber aspect ratio (fiber length to fiber diameter ratio) then $\lambda = 1$ corresponds to the infinitely thin fiber and $\lambda = -1$ to an infinitely thin plate.

The question arises as to whether it is possible to consider the case $\lambda \neq \pm 1$ inside GENERIC. Below, we shall briefly indicate how the non-passive character of the advection can be incorporated into the dissipative part of the GENERIC structure. Instead of the kinematics (22), we shall regard a fiber as a rigid rod. The state variables $(\mathbf{u}(\mathbf{r}), \psi(\mathbf{p}))$ (for the sake of simplicity we shall limit ourselves in this observation to homogeneous suspensions, i.e. $\psi$ is independent of $\mathbf{r}$) are replaced by $(\mathbf{u}(\mathbf{r}), \psi(\mathbf{p}), \mathbf{M}(\mathbf{p}))$, where $\mathbf{M}(\mathbf{p})$ is the angular momentum. The kinematics of these state variables is expressed (see Grmela and Lafleur (1998)) by the Poisson bracket

$$ \{A, B\}^{(\text{rod})} = \int d\mathbf{r} \int d\mathbf{p} \psi \left[ \left( \frac{\partial A_{\psi}}{\partial r_j} B_{u_j} - \frac{\partial B_{\psi}}{\partial r_j} A_{u_j} \right) \right. \\ \left. + M_k \left( \frac{\partial A_{M_k}}{\partial r_j} B_{u_j} - \frac{\partial B_{M_k}}{\partial r_j} A_{u_j} \right) + M_j (A_{\mathbf{M}} \times B_{\mathbf{M}})_j \right. \\ \left. + M_j p_i \left( (A_{\mathbf{M}} \times \frac{\partial B_{M_j}}{\partial \mathbf{p}})_i - (B_{\mathbf{M}} \times \frac{\partial A_{M_j}}{\partial \mathbf{p}})_i \right) \right] \quad (26) $$

where $(\mathbf{a} \times \mathbf{b})_i = \epsilon_{ijk} a_j b_k$ is the vector product of two vectors $\mathbf{a}$, and $\mathbf{b}$; $\epsilon$ is the alternating tensor; $\zeta_i = \epsilon_{ijk} \frac{\partial v_k}{r_j}$ is the vorticity; and $\vartheta_i = \epsilon_{ijk} p_l p_k D_{jl}$. The same calculations as those that we sketched in the text following (22) lead us to

$$ \frac{\partial \psi}{\partial t} = \frac{\partial (\psi \mathbf{p} \times \Phi_{\mathbf{M}})_k}{\partial p_k} \quad (27) $$

We note now that if we let to dissipate the angular momentum $\mathbf{M}(\mathbf{p})$ in such a way that $\Phi_{\mathbf{M}}$ approaches rapidly $-\frac{1}{2}\lambda\vartheta$ then, as an approximation, we can replace in the first equation of (27) $\Phi_{\mathbf{M}}$ by $-\frac{1}{2}\lambda\vartheta$. The resulting equation becomes equivalent to the kinetic equation Eq.(3) with the term $\omega$ missing and the term $\chi(\mathbf{u}, \psi)$ given by (4) and (5).
## 2.4 Comparison of different formulas for the extra stress tensor

The results obtained above in this section imply

$$ \sigma^{(\text{mechCT})}(\psi) = \sigma^{(\text{th})}(\psi) = \sigma^{(\text{GENERIC})}(\psi) \quad (28) $$

$$ \sigma^{(\text{mechDA})}(\psi) \neq \sigma^{(\text{mechCT})}(\psi) \quad (29) $$

and

$$ \sigma^{(\text{mechDA})}(\psi) \neq \sigma^{(\text{th})}(\psi) \quad (30) $$

The first equality in (28) provides an additional argument (in addition to the one presented in Observation 2 in Section 2.2) for $\sigma^{(\text{th}-)} = \sigma^{(\text{GENERIC}-)} = 0$. The inequalities (29) and (30) arise due to the following two reasons: forces and morphology of the microstructure are investigated on different levels of description, and the stresses are investigated in driven systems and not in externally unforced systems considered in Sections 2.2 and 2.3. We now discuss these two reasons in some detail.

**Stresses and the microstructure investigated on different levels of description** We note that the extra stress tensor that is actually observed in experiments (we shall call it hereafter a predicted stress) is $\sigma(\psi|_{\text{Eq.}(3)})$, where by $\psi|_{\text{Eq.}(3)}$ we denote solutions to Eq.(3) with a given $u$. The predicted stress is thus independent of $\psi$ and depends only on the symmetric velocity gradient $D$. Let now the process of finding the solution to Eq.(3) be regarded as a gradual process representing, from the physical point of view, a gradual descend to more macroscopic levels of description, and from the mathematical point of view a gradual restriction of the functions $\psi(r, p)$ to submanifolds of the original state space. The well known example of the restricted submanifold is the local Maxwellian distribution function and its Chapman-Enskog deformations (see more in Grmela (2010)) providing a passage from the phase space kinetic theory to hydrodynamics. Another example is the van Wiechen-Booij (1971) configuration space distribution function $\psi(p) \sim \exp(B_{ij}p_i p_j)$, where the tensor $B$ is proportional to the inverse of the orientation tensor $\tilde{A}$, providing a passage from the configuration space kinetic theory to the orientation tensor theory.

Let the distribution function $\psi$ restricted to such submanifold (we can also see it as a partial solution) be denoted by $\psi|_{\widetilde{\text{Eq.}(3)}}$. If $\sigma(\psi)$ is replaced by $\sigma(\psi|_{\widetilde{\text{Eq.}(3)}})$ then clearly both formulas will give exactly the same predicted stresses. In general, we can think of an infinite number of $\psi|_{\widetilde{\text{Eq.}(3)}}$ representing an exact partial solutions to Eq.(3) and thus an infinite number of formulas for the extra stress tensor, all looking very differently but all implying exactly the same predicted stresses.

The formula $\sigma(\psi|_{\widetilde{\text{Eq.}(3)}})$ can also be seen as being directly derived on a more macroscopic (i.e. less detailed) level than the level on which the distribution function $\psi$ serves as the microstructural state variable. The state variable used on the more macroscopic level is $\psi$ restricted to the manifold on which $\psi =$
$\psi|_{\widetilde{\text{Eq.}(3)}}$). We can thus regard the formula $\boldsymbol{\sigma}^{(\text{mechDA})}$ as obtained in this way. The microstructure equation (3) can indeed be, in principle, investigated on a level of description and the formulas for the extra stress tensor on another different level of description. If however such a route is taken, then the main problem is the compatibility, i.e. the problem of guaranteeing that the physics put into the derivation of the microstructure equation (3) is the same as the physics put (on another level of description) into the derivation of the formula for the extra stress tensor. From the mathematical point of view, this problem is then the problem of proving that $\tilde{\psi}$ for which $\boldsymbol{\sigma}^{(\text{th})}(\tilde{\psi}) = \boldsymbol{\sigma}^{(\text{mechDA})}(\phi)$ is indeed $\psi|_{\widetilde{\text{Eq.}(3)}}$, i.e. a partial solution of Eq.(3).

Can we identify $\tilde{\psi}$ for which $\boldsymbol{\sigma}^{(\text{mechCT})}(\tilde{\psi}) = \boldsymbol{\sigma}^{(\text{mechDA})}(\psi)$? It has been noted in Grmela (2008) that $\boldsymbol{\sigma}^{(\text{mechCT})}(\psi)$ turns into $\boldsymbol{\sigma}^{(\text{mechDA})}(\psi)$ if $\psi$ in $\boldsymbol{\sigma}^{(\text{mechCT})}$ is replaced by

$$ \tilde{\psi}(\boldsymbol{p}) = \text{constant} \exp(p_i p_j D_{ij}) \qquad (31) $$

which is indeed an approximate solution of (3) corresponding to small symmetric velocity gradients $\boldsymbol{D}$. The better (31) approximates solutions to (3) the closer are the predicted stresses calculated from both formulas.

**Externally unforced versus driven suspensions** What is observed in rheological investigations are responses of fluids to external forces. The external forces are typically imposed overall flows and the responses are observed in both stresses and the microstructural morphology. The advection term $\chi(\boldsymbol{u}, \psi)$ in Eq.(3) represent the external influence. Given $\boldsymbol{u}$, solution to (3) provides the response of the internal structure. In order to obtain the stresses induced by the externally imposed flow, the internal structure obtained as a solution to (3) has to be inserted into the formula $\boldsymbol{\sigma}(\xi)$ for the extra stress tensor. Two questions arise: (i) different external forces require different levels of description to investigate the responses, (ii) can the formulas $\boldsymbol{\sigma}(\xi)$ derived in Section 2.2 and 2.3 for externally unforced fluids (from the requirement that such fluids approach equilibrium states) be used in rheological investigations? The first question is clearly out of the scope of this paper. Our starting point is the microstructural equation (3), we do not investigate its domain of applicability. We only note that there are always external forces for which any given level of description (i.e. any given choice of $\xi$ and Eq.(3)) will appear to be inadequate for discussing the fluid response. As for the second question, we see that the formulas for the extra stress tensor derived on the basis of the compatibility with thermodynamics do not have indeed a universal applicability in rheology. Ideally, the formulas for $\boldsymbol{\sigma}$ should be obtained from a mechanical analysis, i.e. from an analysis of forces inside the suspension. Such analysis will, in general, differ for different external forces. Moreover, we have to require that the physics involved in the analysis of the forces should always be the same as the one involved in the derivation of the microstructural equation. The fact that the Chan-Terentjev mechanical formula for the extra stress tensor is the same as the formula derived from thermodynamics is a strong indication of a rather broad applicability of thermodynamic considerations in rheology.
**Comparison with results of rheological measurements** Predicted stresses (formulas for the stress tensor evaluated at solutions to the microstructural equations) that can be compared with results of experimental observations depend on both the microstructural equation (3) and the expression for the extra stress tensor $\sigma(\xi)$. Comparison with experiments does not therefore contribute to answering the question asked in this paper (i.e. the question: given microstructural equations, what is the corresponding to it expression for the extra stress tensor).

**Numerical illustration** We shall now solve the microstructural equation (3) numerically and compare the predicted stresses implied by different formulas $\sigma(\psi)$. The discussion presented above in this subsection indicates the complexity of the issues involved in the comparison. Consequently, the numerical investigation, however detailed, will not really represent a substantial contribution to the discussion. We shall therefore limit ourselves to showing, for selected values of the material parameters involved in (3), that the predicted stresses are indeed different for different formulas but the difference is not large.

In the numerical illustration presented below we assume that the distribution function $\psi(\mathbf{p})$, where $|\mathbf{p}| = 1$, is independent of the position vector $\mathbf{r}$, the velocity gradient of the imposed flow is

$$ \nabla v = \begin{pmatrix} -E & 0 & G \\ 0 & E & 0 \\ 0 & 0 & 0 \end{pmatrix} \qquad (32) $$

where $E$ and $G$ are constants.

We begin the numerical solution by passing to a week form of the kinetic theory equation. We use as degrees of freedom the nodal values of the distribution function over a triangular tetellation of the unit sphere (which is here the configurations space of all possible 3D orientations). The model used in the calculations involves about 2500 degree of freedom. Then we consider a particular velocity gradient given by By using an implicit time integration scheme we calculate at each time step the distribution function. The initial state is given by the isotropic orientation i. e. $\psi(\mathbf{p}) = 1/4\pi$. Simulations are done for $E = 0.1$, $G = 1$ and three values 0.005, 0.05, 0.5 of the parameter $\Lambda$. Figures 1, 2 and 3 depict the distribution function $\psi(\mathbf{p})$, the orientation tensor $\mathbf{A}$ calculated as the second moment of the distribution function (i.e. $A_{ij} = \int d\mathbf{p} p_i p_j \psi(\mathbf{p})$), and predicted stress calculated by using the formulas $\sigma^{(mechDA)}$ and $\sigma^{(mechCT)}$.
PDF at time = 100 with Λ = 0.5, E = 0.1, G = 1

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

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

Kinetic Theory with Λ = 0.5, E = 0.1, G = 1

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

Kinetic Theory with Λ = 0.5, E = 0.1, G = 1

<figure><img src="image_5.png" /><figcaption>Figure 1: Distribution function $\psi(p)$, corresponding to its orientation tensor $A = \int d\mathbf{p}\mathbf{p}\psi(p)$, and predicted stresses $\sigma$ for $E = 0.1$, $G = 1$ and $\Lambda = 0.5$</figcaption></figure>
PDF at time = 100 with $\Lambda = 0.05$, $E = 0.1$, $G = 1$

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

Kinetic Theory with $\Lambda = 0.05$, $E = 0.1$, $G = 1$

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

Kinetic Theory with $\Lambda = 0.05$, $E = 0.1$, $G = 1$

<figure><img src="image_8.png" /><figcaption>Figure 2: Distribution function $\psi(p)$, corresponding to its orientation tensor $A = \int dp \, p \psi(p)$, and predicted stresses $\sigma$ for $E = 0.1$, $G = 1$ and $\Lambda = 0.05$</figcaption></figure>
PDF at time = 100 with Λ = 0.005, E = 0.1, G = 1

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

Kinetic Theory with Λ = 0.005, E = 0.1, G = 1

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

Kinetic Theory with Λ = 0.005, E = 0.1, G = 1

<figure><img src="image_11.png" /><figcaption>Figure 3: Distribution function $\psi(p)$, corresponding to its orientation tensor $A = \int dp pp \psi(p)$, and predicted stresses $\sigma$ for $E = 0.1, G = 1$ and $\Lambda = 0.005$</figcaption></figure>
The smaller is the value of the parameters $\Lambda$ the more important is the role the reversible part $\chi(u, \psi)$ plays in the time evolution. In the solutions to the microstructural equation presented on Figures 1, 2,3 this then manifests itself in more intense elastic-type character of the responses. From the way we have derived in this section the thermodynamic formula $\sigma^{(th)}$, we can anticipate that the difference between $\sigma^{(mechDA)}$ and $\sigma^{(th)}$ will be larger in the situations in which the advection (i.e. the term $\chi(u, \psi)$) plays in the time evolution plays a larger role (i.e. in the situation in which the suspension appears as a more elastic fluid). This is indeed seen on the Figures 1,2,3.

## 3 Orientation tensor theory

In order to gain an additional insight into the relation between the microstructural equation and the formula for the extra stress tensor, we shall repeat the analysis made in Section 2 with another microstructural state variable $\xi$. We choose in this section $\xi = A$, where $A$ is an orientation tensor. From the mathematical point of view, $A$ is a symmetric and positive definite matrix. From the physical point of view, the orientation of the principal axis of the ellipsoid given by the graph of $\langle R, AR \rangle = 1$; $R \in \mathbb{R}^3$ and $\langle \cdot, \cdot \rangle$ denotes the inner product, represents the average orientation of the fibers and the thickness of the ellipsoid the dispersion in the fiber distribution. In order to express the fact that the length of fibers is fixed, we shall constrain $A$ by $tr A = 1$. This constraint plays now the same role as the constraint $|p| = 1$ that we used in kinetic theory in the preceding section. For the sake of simplicity we shall limit ourselves in this section to homogeneous suspensions, i.e. $A$ is independent of the position vector $r$.

As an example of Eq.(3) with $\xi = A$ we take the equation introduced recently in Wang et al. (2008). We shall call it WOT equation. The advection term $\chi$ in this equation is particularly interesting. This makes then also the problem of the correct formula for the extra stress tensor particularly pertinent. The advection term in the WOT equation is given by

$$ \chi(u, A) = \Omega \cdot A - A \cdot \Omega + \lambda D \cdot A + \lambda A \cdot D - 2\lambda[A + (1 - \kappa)(L - M : A)] : D \quad (33) $$

where $\lambda$ is the same parameter as in (5), $\kappa$ is a phenomenological scalar parameter, and $L$, $A$, and $M$ are certain functions (specified in Wang et al. (2008)) of $A$. The part $\omega$ of Eq.(3) in the WOT equation is given by

$$ \omega(A) = 2\kappa C_I \dot{\gamma}(I - 3A) \quad (34) $$

where $\dot{\gamma} = \sqrt{2D : D}$, $I$ is the unit matrix, and $C_I > 0$ is a phenomenological parameter. Regarding the dependence of $L$, $A$, and $M$ on $A$, there are several possibilities to choose from.

In the particular case of the WOT equation corresponding to $\kappa = 1$ and $\lambda = 1$, the advection is passive. The the case $\kappa = 1$ and $\lambda \neq 1$, corresponds to the
well known and often considered active advection. The new WOT advection corresponds to $\kappa \neq 1$. The physics behind this new advection is the following. The influence of the imposed flow on the orientation of the fibers is separated from the influence on the dispersity in the distribution. In terms of the orientation tensor, the separation into orientation and dispersity is made as follows. The tensor $\mathbf{A}$ can always be represented as $(\mathbf{Q}, \mathbf{E})$ where $\mathbf{A} = \mathbf{Q} \cdot \mathbf{E} \cdot \mathbf{Q}^T$; $\mathbf{Q} \cdot \mathbf{Q}^T = \mathbf{I}$, $\mathbf{I}$ is the unit matrix, $()^T$ denotes the transpose, and $\mathbf{E}$ is a diagonal matrix. If we see $\mathbf{A}$ as an ellipsoid (i.e. graph of $\langle \mathbf{R}, AR \rangle = 1$; $\mathbf{R} \in \mathbb{R}^3$ and $\langle \cdot, \cdot \rangle$ denotes the inner product) then $\mathbf{Q}$ characterizes the rotation and $\mathbf{E}$ the shape of the ellipsoid (i.e. a measure of the dispersion). Given an equation governing the time evolution of $\mathbf{A}$, a coupled system of equations governing the time evolution of the two coordinates is implied. Wang et al. (2008) have suggested to modify this system by modifying separately, and in a different way, equations governing the time evolution of $\mathbf{Q}$ and $\mathbf{E}$.

We have now specified the state variables (1) as well as the microstructural equation (3) and we can, following Section 2, proceed to the mechanical, thermodynamical, and GENERIC investigations of the extra stress tensor. Before doing it we shall note that there is also another way to regard the orientation tensor $\mathbf{A}$. The alternative to regarding $\mathbf{A}$ as a state variable with its own autonomous physical interpretation is to regard it as a reduction of the distribution function $\psi(\mathbf{p})$, namely as $A_{ij} = \int d\mathbf{p} \, p_i p_j \psi(\mathbf{r}, \mathbf{p})$. Both the equation governing the time evolution of $\mathbf{A}$ and the expressions for the extra stress tensor could be then seen as reductions of the microstructural time evolution equations and of the expressions for the extra stress tensor discussed in Section 2. This is indeed an interesting avenue to follow but we shall not do it in this paper. Wang et al. (2008) have derived (33) and (34) independently of kinetic equations governing the time evolution of $\psi(\mathbf{p})$. The only reference to the level of kinetic theory appear in their analysis in the dependence of $\mathcal{L}$, $\mathcal{A}$, and $\mathcal{M}$ on $\mathbf{A}$. These functions arise in applying various closures (i.e. mappings $\mathbf{A} \hookrightarrow \psi(\mathbf{p})$) identifying the submanifold of $\psi(\mathbf{p})$ on which the kinetic level reduces to the level on which the orientation tensor serves as the state variable (see e.g. van Wieche-Booij (1971) distribution). In our discussion below, we shall simply consider $\mathcal{L}$, $\mathcal{A}$, and $\mathcal{M}$ as given functions of $\mathbf{A}$. The only time where we shall need their specific form will be in the numerical illustration. All the results regarding the formulas for the extra stress tensor will be obtained below with unspecified functions $\mathcal{L}$, $\mathcal{A}$, and $\mathcal{M}$.

## 3.1 Mechanics

The formula that Wang et al. (2008) used to calculate the fiber contribution to the extra stress tensor is the formula (9) in which $\mathbf{A}$ is not the fourth moment of the distribution function $\psi$ as it in the context of the kinetic theory discussed in the preceding section but it is $\mathcal{A}$ appearing in (33).
## 3.2 Thermodynamics

We now completely follow Section 2.2. with the microstructural state variable $\psi$ replaced by the orientation tensor $\mathbf{A}$. The formula (14) takes in this context the form

$$ \sigma_{ij}^{(th+)} = -\Phi_{A_{kl}} \frac{\partial \chi_{kl}}{\partial \left( \frac{\partial u_i}{\partial r_j} \right)} \quad (35) $$

With $\chi$ given in (33), we obtain

$$ \sigma_{ij}^{(th+)} = \lambda[\Phi_{A_{ik}}A_{jk} + \Phi_{A_{jk}}A_{ik} - 2\Phi_{A_{kl}}\mathbf{A}_{klij} - 2(1-\kappa)\Phi_{A_{kl}}(\mathbb{L}_{klij} - \mathbb{M}_{klmn}\mathbf{A}_{mnij})] \quad (36) $$

Next, we turn to the specification of the free energy $\Phi$. As in Section 2.2, the free energy $\Phi(\mathbf{u}, \mathbf{A})$ is a sum of the kinetic energy $\int dr \frac{\mathbf{u}^2}{2\rho}$ and a term that is independent of $\mathbf{u}$ and that represents the contribution of the internal structure (i.e. the contribution of the fibers). The inequality (16) takes now the form

$$ \Phi_{A_{ij}}\omega_{ij} < 0 \quad (37) $$

To identify the free energy, we proceed in the same way as in Section 2.2. We note that we can cast $\omega(\mathbf{A})$ into the form

$$ \omega_{ij}(\mathbf{A}) = -\Xi_{\Phi_{A_{ij}}} \quad (38) $$

for the dissipation potential

$$ \Xi = \frac{1}{2} \Lambda \Phi_{A_{ij}} A_{jk} \Phi_{A_{ki}} \quad (39) $$

and the free energy

$$ \Phi = \int dr \frac{\mathbf{u}^2}{2\rho} + \frac{nk_B T}{2} (3\text{tr}(\mathbf{A}) - \ln \det(\mathbf{A})), \quad (40) $$

where

$$ \Lambda = \frac{4C_I \kappa \dot{\gamma}}{nk_B T} > 0 \quad (41) $$

To verify the algebra involved we recall that $\frac{\partial(\text{tr}\mathbf{A})}{\partial A_{ij}} = \delta_{ij}$ and $\frac{\partial(\ln \det \mathbf{A})}{\partial A_{ij}} = A_{ij}^{-1}$, where $\mathbf{A}^{-1}$ denote the inverse of $\mathbf{A}$ and $A_{ij}^{-1}$ are its elements.

From the physical point of view, the second term in (40) represents the entropy (first identified in Sarti and Marrucci (1973)), the first term guarantees the constraint $\text{tr}\mathbf{A} = \text{const}$. The equilibrium states $(\mathbf{u}, \mathbf{A})_{\text{eq}}$, i.e. solutions of $\Phi_{\mathbf{u}} = 0$ and $\Phi_{\mathbf{A}} = 0$, are $(\mathbf{u})_{\text{eq}} = 0$ and $(A_{ij})_{\text{eq}} = \frac{1}{3}\delta_{ij}$.

The simplest dissipation potential (by definition a function of the derivative of the free energy with respect to the state variable) satisfying the three properties listed at the end of Section 2.2. is a positive definite quadratic function of $\Phi_{\mathbf{A}}$. The dissipation potential (19) in Section 2.2. has been chosen in this way (in this example $\Phi_{\mathbf{A}}$ is, of course, replaced by $\Phi_{\psi}$). The dissipation function (39) is again a positive definite quadratic function. Its positive definiteness follows from (41) and the positive definiteness of the orientation tensor $\mathbf{A}$.
## 3.3 GENERIC

The advection represented by (33) is passive only for $\lambda = 1$ and $\kappa = 1$. As in kinetic theory, the reversible part of the time evolution governed by (2), (3) is Hamiltonian only in this case. The Poisson bracket expressing this passive advection is given by (see Grmela (1988), Edwards et al. (2003))

$$ \{C, B\}^{(u,A)} = \{C, B\}^{(u)} + \int dr \left[ A_{ij} \left( \frac{\partial C_{A_{ij}}}{\partial r_k} B_{u_k} - \frac{\partial B_{A_{ij}}}{\partial r_k} B_{u_k} \right) \right. \\ \left. + A_{kl} \left( C_{A_{lm}} \frac{\partial B_{u_m}}{\partial r_k} - B_{A_{lm}} \frac{\partial C_{u_m}}{\partial r_k} \right) \right. \\ \left. + A_{km} \left( C_{A_{lm}} \frac{\partial B_{u_l}}{\partial r_k} - B_{A_{lm}} \frac{\partial C_{u_l}}{\partial r_k} \right) \right. \\ \left. - 2A_{ij}A_{kl} \left( C_{A_{kl}} \frac{\partial B_{u_j}}{\partial r_i} - B_{A_{kl}} \frac{\partial C_{u_j}}{\partial r_i} \right) \right] \quad (42) $$

It can easily be verified that $\chi(\mathbf{u}, \mathbf{A})$ implied by (20) and (42) is indeed (33) with $\lambda = 1$ and $\kappa = 1$.

Again, as in kinetic theory, the nonpassive advection (33) with $\lambda \neq 1$ and/or $\kappa \neq 1$ can be seen as an approximation of a time evolution taking place in an extended state space and involving both a passive advection and an irreversible part. In addition to the state variables $(\mathbf{u}, \mathbf{A})$, we adopt also $\mathbf{w}$, having the physical meaning of a conjugate to the gradient of the velocity disturbed by the presence of fibers (see more in Gu and Grmela (2008)). The passive advection of $(\mathbf{u}, \mathbf{A}, \mathbf{w})$ is expressed in the Poisson bracket $\{C, B\}^{u,A,w} = \{C, B\}^{(u,A)} + \{C, B\}^{(w)}$, where $\{C, B\}^{(u,A)}$ is given in (42) and $\{C, B\}^{(w)}$ is given by (see more in Gu and Grmela (2008))

$$ \{C, B\}^{(w)} = A_{jk} (C_{A_{mj}} B_{w_{mk}} - B_{A_{mj}} C_{w_{mk}}) \\ + A_{mk} (C_{A_{mj}} B_{w_{jk}} - B_{A_{mj}} C_{w_{jk}}) \\ + w_{ij} (C_{w_{im}} B_{w_{mj}} - B_{w_{im}} C_{w_{mj}}) \\ + A_{ij} (\partial_m (C_{A_{ij}}) B_{u_m} - \partial_m (B_{A_{ij}}) C_{u_m}) \\ + A_{ki} (C_{A_{im}} \partial_k (B_{u_m}) - B_{A_{im}} \partial_k (C_{u_m})) \\ + A_{km} (C_{A_{im}} \partial_k (B_{u_i}) - B_{A_{im}} \partial_k (C_{u_i})) \\ + w_{ik} (C_{w_{im}} \partial_k (B_{u_m}) - B_{w_{im}} \partial_k (C_{u_m})) \\ - w_{kj} (C_{w_{mj}} \partial_m (B_{u_k}) - B_{w_{mj}} \partial_m (C_{u_k})) \\ + w_{ij} (\partial_m (C_{w_{ij}}) B_{u_m} - \partial_m (B_{w_{ij}}) C_{u_m}) \quad (43) $$

The property $\{C, B\} = -\{B, C\}$ is clearly visible in (43), the Jacobi identity follows from the way the bracket (43) is derived in Gu and Grmela (2008).

By inserting $\{C, B\}^{u,A,w} = \{C, B\}^{(u,A)} + \{C, B\}^{(w)}$ into (20) we obtain equations governing the time evolution of $\mathbf{A}$ and $\mathbf{w}$. The former is

$$ \frac{dA_{ij}}{dt} = A_{ki}(\partial_k(\Phi_{u_j}) + \Phi_{w_{jk}}) + A_{kj}(\partial_k(\Phi_{u_i}) + \Phi_{w_{ik}}) \quad (44) $$
If we add to the right hand side of that latter (see Gu and Grmela (2008) for its explicit form) an appropriate dissipation term that brings $\Phi_w$ rapidly to the stationary state for which $\frac{d\Phi_w}{dt} = 0$ and insert $\Phi_w$ at the stationary state into (44) we arrive at $\chi$ given in (33).

## 3.4 Comparison of different formulas for the extra stress tensor

We have arrived at two formulas for the extra stress tensor. One, $\sigma^{(mechDA)}$, given by (9), derived from mechanical arguments and the other, $\sigma^{(th)}$, given by (36), derived from thermodynamic arguments or from the GENERIC structure. The formulas are different. All what we have said in Section 2.4 can now be repeated (except the arguments that uses the first equality in (28) since we do not have in this section any analog of the mechanical formula $\sigma^{(mechCT)}$ that is the same as the thermodynamic formula). Our conclusion is thus again that the thermodynamic formula is physically more justified than the mechanical formula. An additional indirect argument supporting the expression (36) over the expression (9) is that the parameters $\lambda$ and $\kappa$ quantifying the nonpassivity of the advection do not enter the formula (9) but do enter the formula (36). Details of the advection clearly influence the forces in the suspension and we can therefore expect that the parameters appearing in the advection should show up also in the formula for the extra stress tensor.

As in Section 2., we can suggest that the Dinh-Armstrong formula relates to the thermodynamic formula by

$$ \sigma^{(mechDA)}(\mathbf{A}, \mathbf{u}) \approx \sigma^{(th+)}(\mathbf{A}|_{\widetilde{\text{Eq.}(3)}}) \quad (45) $$

where $\mathbf{A}|_{\widetilde{\text{Eq.}(3)}}$ denotes a partial solution of Eq.(3). In this setting we are however unable to identify $\mathbf{A}|_{\widetilde{\text{Eq.}(3)}}$ for which (45) would hold.

Now we comment about the particular choice (41) of the coefficient $\Lambda$ made in Wang et al. (2008) (i.e. $\Lambda \sim \dot{\gamma}$). While this choice is admissible if both state variables ($\mathbf{u}, \mathbf{A}$) in Eqs.(2),(3) remain unconstrained (i.e. the setting used in thermodynamic and GENERIC arguments), it becomes inadmissible if $\mathbf{u}$ is seen as imposed from outside of the system. The former statement is true because solutions to Eqs.(2),(3) converge to the equilibrium state (i.e. to solutions of $\Phi_{\mathbf{u}} = 0$ and $\Phi_{\mathbf{A}} = 0$) as $t \to \infty$. This then means that the stress converges to zero as $t \to \infty$. If however the overall momentum of the fluid $\mathbf{u}$ is controlled from outside, the stress does not always approaches zero in the absence of flow. This happens in particular if the flow is suddenly stopped. In that case $\dot{\gamma} = 0$, the term $\omega(\mathbf{A}) = 0$, and the orientation of the fibers (i.e. $\mathbf{A}$) ceases to evolve towards equilibrium. Consequently, $\Phi_{\mathbf{A}}$ remains different from zero and thus also the extra stress (36) remains different from zero. This is, of course, a contradiction since the stress tensor appearing in Eq.(2) can be different from zero only in the presence of a flow.

How can be this contradiction (and thus a thermodynamics inadmissibility of the choice (41) of the coefficient $\Lambda$) resolved?
The modification that will solve the problem is to replace $\dot{\gamma}$ in (41) (and thus in the last term on the right hand side of (33)) by $(K + \dot{\gamma})$, where $K > 0$ can be arbitrarily small but different from zero. With this modification, Eqs.(2),(3), (33), (34), (41) become intrinsically compatible. In addition, the cessation of the imposed flow will be followed by relaxation of the fiber orientation and of the stress. We note that, in the case of lamellae replacing the fibers as suspended particles, a good agreement with results of experimental observations has been achieved in Eslami et al. (2007) with a model that is essentially the same as (33) with $\mathbf{A}, \mathbf{L}, \mathbf{M}$ absent and $\Lambda$ that does not involve $\dot{\gamma}$ at all.

### Numerical illustration

Finally, we turn, as in Section 2, to a numerical illustration. Again, intention is not to make a systematic numerical investigation of the differences in the implied predicted stresses but to show, for a selected model parameters, that the predicted stresses are different but the difference is not large.

The WOT equation is considered within the Euler time integration scheme (in the first order). As for the function $\mathbf{A}$, we choose the one that is in Wang et al. (2008) referred to as a hybrid closure. In figure 4 we depict the orientation tensor solving Eq.(3) with (33) and (34) and the predicted stresses calculated with mechanical and thermodynamic formula. The imposed flow is again given by (32). The initial state is given by the isotropic tensor. The values of the parameters are $C_I = 0.01, \kappa = 0.1, G = 1, E = 0.15, 0.2, 0.25$.
<figure><img src="image_12.png" /></figure>

$C_I = 0.01, \kappa = 0.1, E = 0.2, G = 1$

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

$C_I = 0.01, \kappa = 0.1, E = 0.25, G = 1$

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

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

$C_I = 0.01, \kappa = 0.1, E = 0.2, G = 1$

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

$C_I = 0.01, \kappa = 0.1, E = 0.25, G = 1$

<figure><img src="image_17.png" /><figcaption>Figure 4: Orientation tensor and predicted stress for $C_I = 0.01, \kappa = 0.1, G = 1, E = 0.15, 0.2, 0.25$.</figcaption></figure>

Solutions to the microstructural equation and the corresponding to it predicted stresses calculated from $\sigma^{(mechDA)}$ and $\sigma^{(th)}$ are presented on Figure 4 for three types of imposed flows. The flows differ by the strength of their elongational components (the larger is $E$ the larger is the elongational component). We see, for example, that the first normal stress difference calculated from the thermodynamic formula displays more elastic-like behavior than the one calculated from the Dinh-Armstrong mechanical formula. This is indeed what we can anticipate from the way the formulas were derived (the thermodynamic formula takes into account more explicitly the advection) and also from what we have
seen on Figures 1,2,3 in the setting of kinetic theory.

# 4 Concluding remarks

The complex fluid under consideration in this paper is a suspension of fibers. The distribution function $\psi(r, p)$ and the orientation tensor $A$ have been chosen as microstructural state variables. We have followed two routes leading to formulas expressing the stress tensor as a function of the microstructural state variables: 1. mechanical and 2. thermodynamical.

## 1. Mechanics

This is the most common and most frequently used approach to calculate stresses. The physics on which this approach is based is the interpretation of the equation governing the time evolution equation of the fluid velocity as a continuum version of Newton's law. The stress tensor is then calculated by identifying the local surface forces acting on the microstructure. The directness of the relation between what is calculated and what is measured that is inherent to this method is certainly the main reason for its popularity. The main problem that we see in the mechanical approach is that the microstructural time evolution equation on the one hand and the formula for the extra stress tensor on the other hand arise in two parallel and largely independent of each other considerations. Their compatibility is thus in question. The intrinsic compatibility of the microstructural dynamics and macroscopic stresses does not come naturally in the mechanical approach and is never completely guaranteed.

## 2. Thermodynamics

The thermodynamic approach is based on the requirement that solutions to the equations governing the time evolution of the overall fluid velocity and the microstructure agree with one particular experimental observation. The observation is that externally unforced suspensions reach states, called equilibrium states, at which their behavior is well described by classical thermodynamics. The formula for the stress tensor arises as a compatibility condition relating the microstructural and the overall velocity equations. The physics entering the analysis of the microstructural dynamics and the physics entering the analysis of the extra stress tensor is thus guaranteed to be identical. The disadvantage of the thermodynamic approach is that the formula for the extra stress tensor is derived for externally unforced systems and is used then in rheological investigations of driven systems.

We have also explored a stronger version (called GENERIC) of the thermodynamic approach. In addition to the requirement of the agreement of solutions to the governing equations with the experimental observations constituting the basis of equilibrium thermodynamics, GENERIC requires that the reversible parts of the time evolution of the velocity and the microstructure represent together a Hamiltonian time evolution for which the free energy is a constant of motion. We have shown that the GENERIC approach leads to the same formulas for the extra stress tensor as the thermodynamic approach that does not use the Hamiltonian structure. Moreover, we have noted that the Hamiltonian
structure, that is present only when the microstructure is passively advected, can be recovered also in the case of a physically more realistic nonpassive advection. What is needed is an extended setting in which the microstructural state variables involve an additional velocity type variable (as e.g. the angular momentum of the fibers). The nonpassive advection is manifested in such extended setting in the dissipation of the added velocity type state variable, in particular then with its coupling with the dissipation of the fluid velocity.

**Results for fiber suspensions and two choices of microstructural state variables.**

Given a microstructural equation there are, in general, infinitely many corresponding to it formulas for the extra stress tensor that look very differently but are equivalent in the sense that they all imply the same predicted stresses. By a predicted stress we mean the stress obtained by evaluating the expression for the stress tensor at solutions to the microstructural equation. Predicted stresses are the stresses measured in rheological observations. The reason why there are, in general, infinitely many equivalent formulas for the extra stress tensor is that there are, in general, infinitely many partial solutions to the microstructural equation. Let us see the process of getting solution to the microstructural equation as a gradual process in which the space in which the solutions are searched is being gradually restricted. Partial solutions are, in our terminology, the microstructural state variables restricted to such submanifolds. From the physical point of view, the gradual process of restrictions leading to partial solutions can be interpreted as the process of a descend to more macroscopic (i.e. less detailed) levels of description. The submanifolds are also called closures (see more in Section 2.4 and in Grmela (2010)).

In this paper we have seen that the mechanical and the thermodynamic approaches lead to identical formulas for the extra stress tensor only in kinetic theory, i.e. with the distribution function $\psi$ playing the role of microstructural state variable, and only with the mechanical formula derived in Chan and Terentjev (2007). The Dinh-Armstrong mechanical formula, that can be applied also for the orientation tensor $A$ playing the role of microstructural state variable, is different from the formula arising in thermodynamics. Is this formula equivalent (in the sense of the previous paragraph) to the thermodynamic formula? In the context of the kinetic theory we identify a restricted distribution function that, if inserted into the thermodynamic formula, makes it look similar to the Dinh-Armstrong mechanical formula. The identified restricted distribution function is however only an approximative solution (valid for small velocity gradients of the imposed flow) to the kinetic equation. The absence in the Dinh-Armstrong formula of details of the advection of the microstructure is an indication of the non-equivalence of the Dinh-Armstrong and the thermodynamic formulas. The difference is also seen in numerical illustrations developed in both the kinetic theory and the orientation tensor theory. In general, we see that the larger is the role of the time reversible part of the time evolution of the microstructure (i.e. the role of the advection) the larger is the difference in predicted stresses.
Acknowledgement

One of the authors, M.G., would like to thank Natural Sciences and Engi-
neering Research Council of Canada for financial support and Charles Tucker
for stimulating discussion and for suggestions to improve the presentation of
this paper.

References

Arnold, V.I., "Sur la geometrie differentielle des groupes de Lie de dimension infini et ses applications a l'hydrodynamique des fluides parfaits", Ann. Inst. Fourier **16**, 319-361 (1966)

Batchelor, G.K., "The stress system in suspensions of force-free particles", J. Fluid Mech. **41**, 545 (1970)

Batchelor, G.K., "The stress generated in a non-dilute suspension of elon-
gated particles by pure straining motion", J. Fluid Mech. **46**, 593 (1971)

Beris, A.N.and Edwards, B.J., "Thermodynamics of flowing systems", Oxford Univ. Press, Oxford (1994)

Chan, C.J. and Terentjev, E.M., "Non-equilibrium statistical mechanics of liquid crystals: relaxation, viscosity and elasticity", J. Phys. A: Math. Theor. **40**, R103 (2007)

Clebsch, A., "Ueber die Integration der hydrodynamische Gleichungen", J. Reine Angew. Math, **56** 1-10 (1895)

Dinh, S. M., and R. C. Armstrong, A rheological equation of state for semi-
concentrated fiber suspensions, J. Rheol. **28**, 207227 (1984).

Doi, M. "Variational principle for the Kirkwood theory for the dynamics of polymer solutions and suspensions", J. Chem. Phys. **79**, 5080 (1983)

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

Edwards, B.J., Dressler, M., Grmela, M., Ait-Kadi A., "Rheological models with microstructural constraints", Rheol. Acta **42**, 64-72 (2003)

Eslami, H. Grmela, M. Bousmina, M. "A mesoscopic rheological model of polymer/layered silicate nanocomposites", J. Rheol. **51**, 1189-1222 (2007)

Evans, J.G., Ph.D Dissertation, Cambridge University (1975)

Feng, J., Sgalari, G., Leal, L.G., "A theory for flowing polymers with orienta-
tional distribution", J. Rheol. **44**, 1085 (2000)

Grmela, M., "Particle and Bracket Formulations of Kinetic Equations", Con-
temporary Math. **28**, 125-132 (1984), Physics Letters A **102**, 355 (1984)

Grmela, M., "Stress tensor in generalized hydrodynamics", Phys.Letters A, **111A**, 41-44 (1985)

Grmela, M., "Hamiltonian dynamics of incompressible elastic fluids", Phys. Lett. A **130**, 81-86 (1988)
Grmela M. and Ottinger, H.C., "Dynamics and Thermodynamics of Complex Fluids:General Formulation", Phys.Rev.E **56**, 6620-6633 (1997)

Grmela, M., Lafleur, P.G., "Kinetic theory and hydrodynamics of rigid body fluids" J. Chem. Phys. **100**, 6956-6972 (1998)

Grmela, M. "Complex fluids subjected to external influences", J. Non-Newtonian Fluid Mech. **96**, 221-254 (2001)

Grmela, M. "Reciprocity relations in thermodynamics", Physica A **309**. 304-328 (2002)

Grmela, M. "Geometry of mesoscopic dynamics and thermodynamics", J. Non-Newtonian Fluid Mech. **120**, 137-147 (2004)

Grmela, M. "Stress tensor in fiber suspensions", Phys.Lett.A, **372**, 4267-4270 (2008)

Grmela, M. "Multiscale equilibrium and nonequilibrium thermodynamics in chemical engineering", Advances in Chemical Engineering, Vol. **39**, 75-128 (2010) edited by D.H. West and G. Yablonsky, Elsevier Inc.

Gu, J.F. and Grmela, M. "GENERIC model of active advection", J. Non-Newtonian Fluid Mech. **152**, 12-26 (2008)

Kaufman, A.N. "Dissipative Hamiltonian systems: A unifying principle", Phys. Letters A **100**, 419 (1984)

Lipscomb, G. G., M. M. Denn, D. U. Hur, and D. V. Boger, The flow of fiber suspensions in complex geometries, J. Non-Newtonian Fluid Mech. **26**, 297-325 (1988)

Marsden, J.E., Ratiu, T.S. "Introduction to Mechanics and Symmetry", Texts in Applied mathematics., Volume 17, Springer, New York, ((1999)

Morrison, P.J., Bracket formulation for irreversible classical fields, Phys. Letters A **100**, 423 (1984)

Ottinger, H.C. "General projection operator formalism for the dynamics and thermodynamics of complex fluids", Phys. Rev. E **47**, 1416 (1998)

Ottinger, H.C. "Beyond equilibrium thermodynamics", Wiley (2005)

Ottinger, H.C. and Grmela, M. "Dynamics and Thermodynamics of Complex Fluids: Illustration of the General Formalism", Phys.Rev.E **56**, 6633-6650 (1997)

Sarti, G.C., G. Marrucci, "Thermomechanics of dilute polymer solutions: multiple bead-spring model", Chem.Eng.Sci. **28**, 1053 (1973) (1973)

Sircar, S., Wang, Q., "Dynamics and rheology of biaxial liquid crystal polymers in shear flows" J. Rheol. **53**, 819 (2009)

van Wiechen, P.H., Booij, H.C., J. Eng. Math. "A General Solution to the Necklace Model Problem in the Rheology of macromolecules" **5**, 89-98 (1971)

Wang, J., J.F. O'Gara, C.L. Tucker III, "An objective model for slow orientation kinetics in concentrated fiber suspensions: Theory and rheological evidence", J.Rheol. **52** 1179-1200 (2008)

Wang, Q. "A hydrodynamic theory for solutions of nonhomogeneous nematic liquid crystalline polymers of different configurations", J. Chem. Phys. **116**, 9120 (2002)