**1. Introduction**

There is a class of complex fluids which can be considered as suspensions of rod-like particles. Examples include linear polymer solutions, worm-like micelles, FD-virus, liquid crystals, etc. Such a class enjoys interesting properties like anisotropy [1], gradient and vorticity banding [2–4], shear banding instabilities [5], transition between isotropic and nematic phases [6], and cluster formation [7,8].

Here, we formulate a new mathematical model which is good for concentrated suspensions and show that it predicts anisotropy and some other properties of suspensions of rod-like particles. To this end, we study Poiseuille-like shear flows. The practice of pumping oil in pipelines shows that the total oil flux can depend not only on the pressure gradient, but on the history of pumping as well [9]. We establish that the developed model captures such an effect and show its relationship with the hysteresis phenomenon.

Studies of rodlike particles flow in fluids go back to Jeffrey's work on interactions of a floating isolated ellipsoid with unbounded linear shear fluid flow [10]. It turns out that such a particle periodically rotates in Jeffrey's orbits, which depend on the geometry of the particle and its initial orientation. Jeffery's approach was developed further in a number of kinematic models [11,12], which include equations both for particle mass centre and for the direction vector with the help of a third rank shape tensor. Available experiments [13,14] confirmed applicability of the generalized Jeffery equations. Such an approach formed basis for extensions accounting for rod–rod interactions [15,16] and for the prediction fibre alignment distributions in moulded parts [17]. Equations proposed in [18] also allow for governing particles motion in a simplified situation where the rod orientation is restricted to the plane spanned by the direction of shear and the direction of gravity.

In a number of studies, the search for the rheology of suspensions of rodlike particles is reduced to establishing the relationship between stress and rate of strain in shear flows. In [19], starting from experiments with FD-viruses, it was studied how viscosity depends

**Citation:** Shelukhin, V. Flows of Linear Polymer Solutions and Other Suspensions of Rod-like Particles: Anisotropic Micropolar-Fluid Theory Approach. *Polymers* **2021**, *13*, 3679. https://doi.org/10.3390/polym 13213679

Academic Editors: Célio Bruno Pinto Fernandes, Salah Aldin Faroughi, Luís L. Ferrás and Alexandre M. Afonso

Received: 14 September 2021 Accepted: 16 October 2021 Published: 25 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

on concentration, shear rate and ionic strength. An expression for viscosity was derived in [20] with the use of friction coefficients parallel and perpendicular to the rod axis. We refer the reader to detailed description of viscosity representation formulas to [1,19,21]; the viscosity dependence on shear rate is also discussed there.

Our approach is different. We use methods of mechanics of continua by applying conservation laws only and not involving the concept of particle direction. To take into account particle rotation and form, we apply the theory of micropolar fluids, which allows for particle microinertia [22]. According to this theory, which is a part of rational mechanics, any infinitesimal volume contains sufficiently many particles. This is why such an approach is applicable for suspensions with a high concentration of particles. As is proved within the micropolar fluid theory in [23], it is due to particle rotation that the Segre–Silberberg effect occurs. Such an effect is known as a tubular pinch phenomenon, stating that particles tend to migrate towards a concentric annular region for the laminar flow of neutrally buoyant dilute suspension of rigid spheres through a circular tube [24]. There is one more effect caused by particle rotation and rotational diffusion. This is the separation of particles when flowing between two concentric rotating cylinders [25].

The micropolar fluid theory allows for intrinsic rotations and micro-inertia thanks to the concept of the Cosserat continuum where each material point is treated as a rigid body [26]. We formulate anisotropic constitutive law by including the micro-inertia tensor into stress/rate of strain relationships. Such an idea of anisotropy was first formulated in [27]. In a great number of papers, rotation of the particles is neglected and the anisotropy is taken into account by using the differences of normal stresses [28].

In the micro-polar fluid theory allowing for internal spins, stress tensor loses symmetry, couple stress appears, and the angular momentum equation should be included into conservation laws. Formulation of rheological constitutive laws in the present paper involves introduction of new viscosities both relative to the Cauchy stress tensor and to the couple stress tensor. Skew-symmetric and anisotropic viscosities are introduced in addition to the common shear viscosity, which we call here symmetric viscosity. While there are experiments [29] and theories [25] to determine the skew-symmetric viscosity, the question of measuring the anisotropic viscosity remains open. We cannot quantitatively confirm our equations by three-dimensional experiments, since the calculations were carried out on the basis of one-dimensional flows. However, we prove that it is precisely due to the anisotropic viscosity that these equations capture such effects as hysteresis, shear gradient banding instability and phase transition.

The flow of short polymer chains in the channels can be regarded as an example of the applicability of the method outlined in this work. The flow of such polymers between graphite surfaces is studied in [30] by the molecular dynamic simulation technique. It is established there that the polymer chains exhibit preferential alignment of oligomers parallel to the surfaces with increasing shear rate. Though in the present paper it is assumed that rods lie in the plane orthogonal to the bounding planes, we predict like in [30] that the apparent viscosity shows an oscillatory behaviour and its variation versus the shear rate corresponds to the shear thinning phenomenon. We perform calculations of the simple flow depending on one variable only; nevertheless, we capture appearance and instability of the nematic phase. More complicated phase transition was addressed in [31] for colloidal suspensions in water; the nucleation of a kagom lattice from solution was detected.

The goal of Section 2 is to remind foundations of the micropolar fluid theory and formulate conservation and constitutive laws obeying the basic principles of thermodynamics. In Section 3, we derive equations for one-dimensional Poiseuille-like shear flows. Finally, in Section 4, we perform calculations explaining different phenomena. In addition to anisotropy, the calculations predict gradient banding instability, phase transition between isotropic and nematic phases, sustained temporal oscillations of macroscopic viscosity, shear thinning and hysteresis. For the flow between two planes, we also establish that the total flow rate depends not only on the pressure gradient, but on the history of its variation as well.

#### **2. Anisotropic Micropolar Fluids**

We remind basic notions of the micro-polar fluid theory. Given a material point labelled by the Lagrangian coordinate-vector *ξ*, the position vector **x**(*t*, *ξ*) at the time instant *t* in the three dimensional Euclidean space jointly with orthogonal director-vectors **d***i*(*t*, *ξ*), *i* = 1, 2, 3, are assigned to such a point to treat it as a rigid body. Orientation of **d***<sup>i</sup>* is controlled by an orthogonal tensor Q(*t*, *ξ*):

$$\mathbf{d}\_{i}(t,\mathfrak{F}) = \mathbf{Q}(t,\mathfrak{F})\mathbf{d}\_{i}(0,\mathfrak{F}), \quad \mathbf{Q}^{\*}\mathbf{Q} = \mathbf{Q}\mathbf{Q}^{\*} = I.$$

Here, *I* is the identity matrix with the elements *δ i j* , Q∗ is the adjoint matrix, (Q∗ )*ij* = Q*ji*. The skew-symmetric tensor Ω(*t*, **x**) = Q*t*Q∗ defines the particle's rotation with the angular velocity

$$\omega(t, \mathbf{x}) = \mathbf{e}\_i \times (\Omega \mathbf{e}\_i) / 2 = \mathbf{e} : \Omega / 2, \quad (\Omega \mathbf{a})\_i = \Omega\_{ij} a\_j \,\forall \, \mathbf{a} \in \mathbb{R}^3 \,\forall i$$

where {**e***i*} 3 1 is an orthonormal basis in R<sup>3</sup> and *e* is the Levi–Civita third order tensor,

$$\varepsilon(\mathbf{a}, \mathbf{b}, \mathbf{c}) = \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}), \ \mathbf{e}\_{\mathrm{i}} \times \mathbf{e}\_{\mathrm{j}} = \varepsilon\_{\mathrm{sij}} \mathbf{e}\_{\mathrm{s}\prime} \ \varepsilon\_{\mathrm{sij}} = \varepsilon(\mathbf{e}\_{\mathrm{s}\prime} \mathbf{e}\_{\mathrm{i}} \ \mathbf{e}\_{\mathrm{j}}), \ (\mathbf{e} : \Omega)\_{\mathrm{i}} = \varepsilon\_{\mathrm{i}\prime \mathrm{k}} \Omega\_{\mathrm{j}\mathrm{k}}.$$

Given the velocity field **v**(*t*, **x**) = **x***t*(*t*, *ξ*), we introduce the rate of strain tensors [22].

$$B = \nabla \mathbf{v} - \Omega, \quad A = \nabla \omega\_{\nu}$$

where (∇**v**)*ij* = *∂vi*/*∂x<sup>j</sup>* . Observe that both *B* and *A* are objective relative to the change of frame of references.

Let *ρ*, **s**, *T* and *N* stand for the mass density, the specific internal spin, the Cauchy stress tensor and the angular moment tensor, respectively. We introduce the material derivative *ρ*˙ (or *dρ*/*dt*) related to the velocity field **v** as follows

$$
\dot{\rho} = \frac{\partial \rho}{\partial t} + v\_i \frac{\partial \rho}{\partial \mathbf{x}\_i} \quad \text{or} \quad \dot{\rho} = \frac{\partial \rho}{\partial t} + (\mathbf{v} \cdot \nabla)\rho. \tag{1}
$$

Conservation laws of mass, momentum and angular momentum are given by the equations

$$
\rho + \rho \text{div}\, \mathbf{v} = \mathbf{0}, \tag{2}
$$

$$
\rho \,\,\dot{\mathbf{v}} = \text{div}\,T + \rho \,\mathbf{f},\tag{3}
$$

$$
\rho \,\dot{\mathbf{s}} = \text{div}\,\mathbf{N} - \mathbf{e} : T + \rho \,\mathbf{l}, \tag{4}
$$

where **f** is the mass force density, *l* is the mechanical couple density, and

$$(\text{div}\,T)\_{\dot{i}} = \partial T\_{\dot{i}\dot{j}}/\partial x\_{\dot{j}}.$$

Observe that the stress tensor *T* is not symmetric. Given an orthonormal basis {**e***i*} 3 1 , the vector

$$\mathbf{t} = \mathbf{e}\_l \times (\mathbf{T} \cdot \mathbf{e}\_l) = \boldsymbol{\varepsilon} : \mathbf{T}$$

does not depend on the choice of the basis and it is a stress symmetry defect measure in the sense that the equality **t** = 0 implies T ∗ = T and vice versa. By the definition of **t**, we have the formula

$$\mathbf{t} \cdot \boldsymbol{\omega} = \mathbf{T} : \Omega. \tag{5}$$

The internal specific spin is defined by the formula **s** = *Jω*, where the symmetric inertia tensor *J*[cm<sup>2</sup> ] obeys the identity [22]

$$
\mathbf{J} - \boldsymbol{\Omega}\mathbf{J} + \mathbf{J}\boldsymbol{\Omega} = \mathbf{0}.\tag{6}
$$

Before proceeding to constitutive laws, we address the thermodynamics issue. Given a specific internal energy *<sup>e</sup>*, the total energy *<sup>E</sup>* <sup>=</sup> e <sup>+</sup> **<sup>v</sup>** · **<sup>v</sup>**/2 <sup>+</sup> **<sup>s</sup>** · *<sup>ω</sup>*/2 satisfies the equation [32]

$$\rho \, \dot{E} = \text{div}\,(T^\* \mathbf{v} + N^\* \omega - \mathbf{q}) + \rho \, \mathbf{f} \cdot \mathbf{v} + \rho \, \mathbf{l} \cdot \omega,\tag{7}$$

where **q** is the heat flux obeying the Fourier law **q** = − κ∇*θ*, with κ standing for the heat conductivity. Generally, internal energy *e* depends on *ρ*, *η* and *J*, *e* = *e*(*ρ*, *η*, *J*), where *η* is the specific entropy. It is common knowledge that absolute temperature and pressure are defined by the derivatives *θ* = *∂e*/*∂η*, *p* = *ρ* <sup>2</sup>*∂e*/*∂ρ* respectively [33]. We calculate that

$$\mathbf{e} = \mathbf{e}\_{\rho}\boldsymbol{\phi} + \mathbf{e}\_{\eta}\boldsymbol{\eta} + \nabla\_{\mathsf{J}}\mathbf{e} : \boldsymbol{\mathfrak{f}}, \quad \text{where} \quad \nabla\_{\mathsf{J}}\mathbf{e} : \boldsymbol{\mathfrak{f}} = \left(\mathfrak{f}\_{\mathrm{ij}}\frac{\partial}{\partial I\_{\mathrm{ij}}}\right)\mathbf{e}.$$

From the rheological point of view, the internal energy e(*ρ*, *η*, *J*) should be an isotropic function of *<sup>J</sup>*. Hence, <sup>∇</sup>Je is also an isotropic function of *<sup>J</sup>*; it implies that [34]

$$
\nabla\_\mathbf{J} \mathbf{e} = \mathfrak{a}\_0 \mathbf{I} + \mathfrak{a}\_1 \mathbf{J} + \mathfrak{a}\_2 \mathbf{J}^2,\tag{8}
$$

where the scalar functions *α<sup>i</sup>* depend on invariants of *J*. Now, it follows from (6) and (8) that <sup>∇</sup>Je : ˙*J* = 0.

We use Equation (6) to calculate that

$$\frac{d}{dt}(\mathbf{s}\cdot\boldsymbol{\omega}) = \dot{\mathbf{s}}\cdot\boldsymbol{\omega} + \mathbf{s}\cdot\dot{\boldsymbol{\omega}} = 2\dot{\mathbf{s}}\cdot\boldsymbol{\omega}.$$

Hence,

$$
\rho \dddot{\mathbf{E}} = \frac{p}{\rho} \rho + \rho \theta \dot{\eta} + \rho \dot{\mathbf{v}} \cdot \mathbf{v} + \rho \dot{\mathbf{s}} \cdot \boldsymbol{\omega} = -p \text{div} \, \mathbf{v} + \rho \dot{\mathbf{v}} \cdot \mathbf{v} + \rho \dot{\mathbf{s}} \cdot \boldsymbol{\omega}.\tag{9}
$$

Multiplying Equations (3) and (4) by **v** and *ω*, respectively, we arrive at the energetic equality

$$\rho \cdot \mathbf{v} \cdot \dot{\mathbf{v}} + \rho \,\omega \cdot \dot{\mathbf{s}} = \text{div}(T^\* \mathbf{v} + N^\* \omega) - T : B - N : A + \rho \,\mathbf{f} \cdot \mathbf{v} + \rho \,\mathbf{l} \cdot \omega \tag{10}$$

With *S* standing for the viscous part of the stress tensor *T*, we write the representation formula *T* =−*pI* + *S*. Hence, *T* : *B* = −*p*div **v**+ *S* : *B*, where *T* : *B* = *TijBij*. Now, it follows from (7), (9) and (10) that the entropy equation

$$
\rho \eta \dot{\eta} + \text{div} \left( \frac{\mathbf{q}}{\theta} \right) = \frac{R}{\theta} \,\text{'}\tag{11}
$$

holds, with the function

$$R = S:B + N:A + \frac{\varkappa |\nabla \theta|^2}{\theta}$$

standing for the entropy production. The second law of the thermodynamics *R* ≥ 0 can be formulated as

$$S: \mathcal{B} + \mathcal{N}: A \ge 0. \tag{12}$$

In what follows, we use the notations

$$B\_s = \frac{B + B^\*}{2}, \quad B\_a = \frac{B - B^\*}{2}$$

for the symmetric and skew-symmetric parts of *B*. We formulate anisotropic constitutive laws as follows:

$$S = 2\mu\_s B\_s + 2\mu\_a B\_a + 2\mu\_{an}\sigma^2 I B\_\prime \quad \text{N} = \frac{2\nu}{\sigma^2} A + 2\nu\_{an} A \,\text{J}\_\prime \tag{13}$$

where *µ<sup>s</sup>* , *<sup>µ</sup>a*, *<sup>µ</sup>an*, *<sup>ν</sup>*, *<sup>ν</sup>an*[g/(cm · s)] are the viscosities and *<sup>σ</sup>*[cm−<sup>1</sup> ] is the specific particles surface area. The first rheological equation in (13) suggests that the contributions of the symmetric part *B<sup>s</sup>* and skew-symmetric part *B<sup>a</sup>* of the rate of strain tensor *B* into local stress state are different. The fact that both *S* and *N* depend on *J* implies anisotropy. Such an approach was first formulated in [22]. Observe that the objectivity of the *S* and *N* results form the objectivity of *B*, *A* and *J* [22].

Due to the symmetry of *J*, one can verify that

$$JB:B = \sum\_{1}^{3} \lambda\_{j} |B^\* \mathbf{e}\_{j}|^2 \rho$$

where **e***<sup>j</sup>* and *λ<sup>j</sup>* are the eigenvectors and the eigenvalues of *J*. Observe that *λ<sup>j</sup>* ≥ 0 provided each suspension particle enjoys an axis of rotational symmetry. For such suspensions, we find that

$$S: B = 2\mu\_s B\_s: B\_s + 2\mu\_a B\_a: B\_a + 2\mu\_{an}\sigma^2 \sum\_{1}^{3} \lambda\_j |B^\* \mathbf{e}\_j|^2 \ge 0.1$$

Similarly, one can verify that

$$Af:J = \sum \lambda\_j |A\mathbf{e}\_j|^2.$$

Thus, the constitutive laws (13) satisfy the thermodynamic restriction (12), provided the suspension particles are axially symmetric.

#### **3. Poiseuille Flows**

We consider one-dimensional flows along the horizontal *x*-axis in the vertical layer |*y*| < *H* between two parallel planes under the prescribed pressure gradient ∇*p* = (*px*, 0, 0), *px*(*t*) < 0, Figure 1a. In this case, *v*<sup>2</sup> = *v*<sup>3</sup> = 0, *v*<sup>1</sup> = *v*(*y*, *t*). We assume that the suspension particles are the rods of the same size; they lie in the plane *z* = 0 and rotate around the *z*-axis. Hence, *ω* = (0, 0, *ω*).

**Figure 1.** (**a**) Schema of particle's position in one-dimensional flows. (**b**) The cylinder approximation of the rod-like particle.

Let us describe the micro-inertia tensor *J*. First, we consider a cylinder *V*<sup>0</sup> stretched along the *y*-axis with the height *h* and the radius *r*, Figure 1b. By definition, the inertia tensor *J*(*V*0) of *V*<sup>0</sup> is given by the formula

$$J(V\_0) = \int\_{V\_0} |\mathfrak{f}|^2 \cdot I - \mathfrak{f} \otimes \mathfrak{f} \, d\mathfrak{f} \quad \text{or} \quad J\_{ij}(V\_0) = \int\_{V\_0} |\mathfrak{f}\_k|^2 \delta\_{ij} - \mathfrak{f}\_i \mathfrak{f}\_j \, d\mathfrak{f}\_i$$

where *I* is the identity matrix and **a** ⊗ **b** stands for the tensor product of two vectors **a** and **b**, (**a** ⊗ **b**)*ij* = *aib<sup>j</sup>* . Calculations reveal that

$$f(V\_0) = \begin{pmatrix} r^2/4 + h^2/3 & 0 & 0 \\ 0 & r^2/2 & 0 \\ 0 & 0 & r^2/4 + h^2/3 \end{pmatrix}.$$

In the limit as *r* → 0, we obtain the inertia tensor of the rod particle stretched along the *y*-axis:

$$J^0 = \lim\_{r \to 0} J(V\_0) = j\_0 \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}. \quad j\_0 = h^2/3.$$

Let *J*(*V*) be the inertia tensor of the cylinder *V*, which is produced by rotation of *V*<sup>0</sup> around the *z*-axis by the angle *ϕ* counted from the axis *y* counter-clockwise, see Figure 1a. By definition of the spin **s**, we find that

$$\mathbf{s} = J(V) \cdot \boldsymbol{\omega} = \int\_{V} \mathbf{x} \times (\boldsymbol{\omega} \times \mathbf{x}) \, d\mathbf{x} = Q\_{\boldsymbol{\varphi}} l^{0} Q\_{\boldsymbol{\varphi}}^{\*} \cdot \boldsymbol{\omega} \, d\mathbf{x}$$

where *Q<sup>ϕ</sup>* is the orthogonal matrix such that

$$
\Omega = \dot{Q}\_{\phi} \mathcal{Q}\_{\varphi \prime}^{\*} \quad \Omega \cdot \mathbf{h} = \omega \times \mathbf{h} \,\forall \, \mathbf{h} \,, \quad (\mathcal{Q}^{\*})\_{\mathrm{ij}} = \mathcal{Q}\_{\mathrm{ji}} \,\tag{14}
$$

$$Q\_{\varphi} = \begin{pmatrix} \cos \varphi & -\sin \varphi & 0 \\ \sin \varphi & \cos \varphi & 0 \\ 0 & 0 & 1 \end{pmatrix}, \quad \Omega = \begin{pmatrix} 0 & -\omega & 0 \\ \omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \tag{15}$$

with "dot" standing for the material derivative (1) related to the velocity vector **v**. Thus, *J*(*V*) = *Q<sup>ϕ</sup> J*(*V*0)*Q*<sup>∗</sup> *ϕ* . In the limit as *r* → 0, we obtain the inertia tensor *J*(*ϕ*) of the rod particle with the position angle *ϕ*:

$$J(\boldsymbol{\varrho}) = \mathbb{Q}\_{\boldsymbol{\varrho}} \boldsymbol{\varrho}^{0} \mathbb{Q}\_{\boldsymbol{\varrho}}^{\*} = \boldsymbol{j}\_{0} \begin{pmatrix} \cos^{2} \boldsymbol{\varrho} & \sin \boldsymbol{\varrho} \cos \boldsymbol{\varrho} & 0 \\ \sin \boldsymbol{\varrho} \cos \boldsymbol{\varrho} & \sin^{2} \boldsymbol{\varrho} & 0 \\ 0 & 0 & 1 \end{pmatrix} \quad \frac{\partial \boldsymbol{\varrho}}{\partial t} = \boldsymbol{\varrho}\_{l} = \boldsymbol{\omega}. \tag{16}$$

Given an initial distribution of particle's angles *ϕ*0(*y*), we denote the initial microinertia tensor by *J*0(*y*) = *J*(*ϕ*0(*y*)):

$$J|\_{t=0} = J\_0. \tag{17}$$

For the described one-dimensional flows, the material derivative ˙*J* reduces to the time derivative *J<sup>t</sup>* . With the use of Equation (15), (6) can be written as follows:

$$\frac{\partial}{\partial t}I\_{11} = -2\omega J\_{12\prime} \quad \frac{\partial}{\partial t}I\_{12} = \omega (f\_{11} - f\_{22}), \quad \frac{\partial}{\partial t}I\_{22} = 2\omega J\_{12\prime} \tag{18}$$

and *Jij* = 0 otherwise. Observe that

$$\frac{\partial J\_{\mathrm{ij}}}{\partial t} = \frac{\partial J\_{\mathrm{ij}}}{\partial \varphi} \frac{\partial \varphi}{\partial t} = J'\_{\mathrm{ij}} \omega, \quad J'\_{\mathrm{ij}} = \frac{\partial J\_{\mathrm{ij}}}{\partial \varphi}.$$

Hence, the system (18) is equivalent to

$$\mathbf{J}'\_{11} = -2\mathbf{J}\_{12} \quad \mathbf{J}'\_{12} = \mathbf{J}\_{11} - \mathbf{J}\_{22} \quad \mathbf{J}'\_{22} = 2\mathbf{J}\_{12} \quad \mathbf{J}\_{ij}|\_{\boldsymbol{\phi} = \boldsymbol{\varphi}\_0} = \mathbf{J}\_0. \tag{19}$$

One can verify that the matrix *J* in (16) solves the system (19).

We calculate the rate of strain tensors and find that

$$B = \begin{pmatrix} 0 & v\_y + \omega & 0 \\ -\omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad B^\* = \begin{pmatrix} 0 & -\omega & 0 \\ v\_y + \omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \tag{20}$$

$$B\_s = \begin{pmatrix} 0 & v\_y/2 & 0 \\ v\_y/2 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad B\_d = \begin{pmatrix} 0 & v\_y/2 + \omega & 0 \\ -v\_y/2 - \omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix},$$

$$j\_0^{-1} J B = \begin{pmatrix} -\omega\cos\varphi\sin\varphi & \cos^2\varphi(v\_y + \omega) & 0 \\ -\omega\sin^2\varphi & \cos\varphi\sin\varphi(v\_y + \omega) & 0 \\ 0 & 0 & 0 \end{pmatrix},$$

$$A = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & \omega\_y & 0 \end{pmatrix}, \quad j\_0^{-1} A f = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ \omega\_y \cos\varphi\sin\varphi & \omega\_y \sin^2\varphi & 0 \end{pmatrix},\tag{21}$$

Let us denote

$$\varepsilon\_1 = \frac{\mu\_a}{\mu\_s}, \quad \varepsilon\_{20} = \frac{\mu\_{an} j\_0 \sigma^2}{\mu\_s}, \quad \varepsilon\_{30} = \frac{\nu\_{an} j\_0 \sigma^2}{\nu}, \quad B^0 = B\_s + \varepsilon\_1 B\_a + \varepsilon\_{20} I B\_a$$

Calculations reveal that matrix *B* 0 is equal to

$$
\begin{pmatrix}
\frac{v\_{y}(1-\varepsilon\_{1})}{2} - \omega(\varepsilon\_{1} + \varepsilon\_{20}\sin^{2}\varphi) & \varepsilon\_{20}\cos\varphi\sin\varphi(v\_{y} + \omega) & 0\\\
0 & 0 & 0
\end{pmatrix}.
$$

We consider incompressible fluids with the assumption *ρ* = const. For one-dimensional flows, the incompressibility condition div **v** = 0 is satisfied automatically. Other conservation laws (3) and (4) become

$$
\rho\_l = \omega, \quad \rho v\_l = -p\_x + \frac{\partial S\_{12}}{\partial y}, \quad \rho j\_0 \omega\_l = \frac{\partial N\_{32}}{\partial y} + S\_{21} - S\_{12}.\tag{22}
$$

For one-dimensional flows, the constitutive laws (13) reduce to

$$S\_{\rm ij} = 2\mu\_s B\_{\rm ij}^0 \quad \text{N}\_{\rm 32} = 2\frac{\nu}{\sigma^2} A\_{\rm 32} + 2\nu\_{\rm an} (Af)\_{\rm 32}. \tag{23}$$

Observe that

$$S\_{21} - S\_{12} = 2\mu\_s (B\_{21}^0 - B\_{12}^0), \\ B\_{21}^0 - B\_{12}^0 = -\upsilon\_y (\varepsilon\_1 + \varepsilon\_{20} \cos^2 \varphi) - \omega (2\varepsilon\_1 + \varepsilon\_{20}).$$

We formulate boundary conditions at |*y*| = *H* as follows:

$$\mathbf{v} = 0, \quad \omega = \frac{\alpha}{2} \nabla \times \mathbf{v}, \quad 0 \le \alpha \le 1. \tag{24}$$

The first condition in (24) states that velocity obeys the no-slip condition. The second condition in (24) has the meaning that the micro-rotation *ω* depends linearly on the macrorotation ∇ × **v**/2 [25].

Let *V* and *T* stand for the velocity and time reference values. We denote Ω = *V*/*H* and choose *T* = 1/Ω. We introduce dimensionless variables as follows:

$$S' = \frac{1}{2\mu\_s\Omega} S, \quad N' = \frac{H\sigma^2}{2\nu\Omega} N, \quad B'^0 = \frac{B^0}{\Omega}, \quad Re = \frac{H^2\rho\Omega}{\mu\_s} J$$

$$y' = \frac{y}{H'}, \\ v' = \frac{v}{V'}, \\ t' = \frac{t}{T}, \\ \omega' = \frac{\omega}{\Omega}, \\ \Pi = \frac{|p\_x|H^2}{2V\mu\_s}, \\ \gamma = \frac{\nu}{H^2\mu\_s\sigma^2}.$$

In new variables, Equation (22) become

$$
\varphi\_{l'} = \omega', \quad \frac{\text{Re}}{2} v'\_{l'} = \Pi + \frac{\partial S'\_{12}}{\partial y'}, \quad \frac{\text{Re} j\_0}{2H^2} \omega'\_{l'} = \gamma \frac{\partial N'\_{32}}{\partial y} + S'\_{21} - S'\_{12}.\tag{25}
$$

In what follows, we consider quasi-steady slow flows. Neglecting terms with small Reynolds numbers in (25), we arrive at the equations

$$
\varphi\_{l'} = \omega', \quad 0 = \Pi + \frac{\partial S\_{12}'}{\partial y'}, \tag{26}
$$

$$0 = \gamma \frac{\partial}{\partial y'} \left[ \omega\_{y'}^{\prime} (1 + \varepsilon\_{30} \sin^2 \varphi) \right] - \left[ v\_{y'}^{\prime} (\varepsilon\_1 + \varepsilon\_{20} \cos^2 \varphi) + \omega^{\prime} (2\varepsilon\_1 + \varepsilon\_{20}) \right] \tag{27}$$

where

$$S\_{12}' = \frac{v\_{y'}'(1 + \varepsilon\_1 + 2\varepsilon\_{20}\cos^2\varphi)}{2} + \omega'(\varepsilon\_1 + \varepsilon\_{20}\cos^2\varphi).$$

Because of the symmetry conditions

$$
\omega^{\prime}(-y^{\prime},t^{\prime}) = v^{\prime}(y^{\prime},t^{\prime}), \quad \omega^{\prime}(-y^{\prime},t^{\prime}) = -\omega^{\prime}(y^{\prime},t^{\prime}), \quad \varphi(-y^{\prime},t^{\prime}) = \varphi(y^{\prime},t^{\prime})
$$

we consider flows only in the upper half-layer 0 < *y* 0 < 1. In such a case the initial and boundary conditions take the form

$$\varphi|\_{t'=0} = \varphi\_0(y'), \ v'(1) = 0, \ v\_{y'}'(0) = 0, \ \omega'(1) = -0.5\omega v\_{y'}'(1), \ \omega'(0) = 0. \tag{28}$$

To perform numerical solution, one should fix the dimensionless parameters Π, *ε*1, *ε*20, *ε*30, *γ*, *α*.

#### **4. Results of Calculations**

Here, we address the system (26)–(28) by applying the Wolfram Mathematica solver.

It is well known in many complex fluids that a shear banding effect occurs when applied shear stress is above some critical value [4,35,36]. Such a phenomenon is characterized by coexisting bands of different shear rates and /or viscosities. Depending on the directional alignment of the banded structure, there are two types of shear banding for suspensions of rod-like particles: gradient banding and vorticity banding [2–4]. In the case of gradient banding, the flow separates into bands of different shear rates along the gradient direction. With reference to the coordinate system of Figure 1a, *x* is the flow direction along the velocity vector **v** = (*v*, 0, 0), *y* is the gradient direction along which the flow has non-zero derivative *∂v*/*∂y*. The *z*-axis is the vorticity direction along the non-zero macro-vorticity vector ∇ × **v**.

The system (26)–(28) cannot be applied for description of the vorticity banding since the corresponding one-dimensional flow does not depend on the *z*-variable. However, calculations reveal that the system (26)–(28) can really capture the gradient banding. Figure 2 depicts appearance of gradient banding when shear stress increases; calculations are performed at *t* = 10 for

$$
\varepsilon\_1 = 1, \quad \varepsilon\_{20} = \mathfrak{I}, \quad \varepsilon\_{30} = \mathfrak{I}, \quad \gamma = 1.3, \quad \mathfrak{a} = 0.3, \quad \varrho\_0 = 0. \tag{29}
$$

Intervals where *ϕ*(*y*) = const or *ω*(*y*) = const correspond to the nematic phase. The profiles of the intrinsic angular velocity *ω* at Figures 2b and 3 imply appearance and instability of the nematic phase. Figure 4b depicts the phase transition from the isotropic phase to the nematic phase.

**Figure 2.** From top to bottom, profiles of the dimensionless velocity *v*(*y*) and dimensionless microspin *ω*(*y*) on the upper half-layer 0 < *y* < 1 at dimensionless time *t* = 10 for dimensionless pressure gradient (**a**) Π = 0.85 and (**b**) Π = 2.85. Gradient banding development is observed at high pressure gradients (**b**).

**Figure 3.** Gradient banding instability with respect to time. From top to bottom, dimensionless velocity *v*(*y*) and dimensionless micro-spin *ω*(*y*) profiles at Π = 2.85 for different dimensionless times *t* = 15 (**a**) and *t* = 25 (**b**). Values of other parameters are as in the data list (29).

**Figure 4.** Gradient banding instability with respect to initial particles orientation. From top to bottom, profiles of dimensionless velocity *v*(*y*) and dimensionless micro-spin *ω*(*y*) at Π = 2.85 and at *t* = 15 for initial *ϕ*0(*y*) = 0 (**a**) and *ϕ*0(*y*) = 4*y* + 9*y* 2 (**b**). Values of other parameters are as in the data list (29).

Figure 3 shows gradient banding instability with respect to time. A treatment of time dependent phenomena for worm-like micelles can be found in [5].

It turns out that the gradient banding is also unstable with respect to initial particles orientation. When passing from spatially homogeneous initial orientation of particles *ϕ*0(*y*) = 0 to a spatially heterogeneous orientation (like *ϕ*0(*y*) = 4*y*+ 9*y* 2 ), the gradient banding effect becomes more pronounced, see Figure 4.

Many shear banding systems display oscillations or irregular fluctuations. Example systems include worm-like micelles [37]. Within the developed anisotropic model, one can observe a chaotic behaviour of the shear velocity even at a constant applied pressure gradient, see Figure 5. Basically, it is due to anisotropic viscosities in the rheological constitutive laws (13).

**Figure 5.** Time variation of the velocity in the middle of the channel at a constant pressure gradient in dimensionless variables (**a**) for homogeneous transversal initial particles orientation and (**b**) for non-homogeneous initial particles orientation along the channel.

Next, we consider questions motivated by oil transportation through pipelines. To optimize pumping, additives are used that change the microstructure of oil. As a result, it is discovered that friction factor can depend not only on oil discharge, but on its prehistory as well [38]. It turns out that the smallest friction losses are achieved by decreasing rather than increasing the flow rate to a predetermined level [9]. Let us show that the developed mathematical model in Section 1 captures such an a effect.

First, we establish that the system (26)–(28) does not provide one-to-one correspondence between the pressure gradient Π and the total fluid flux *Q* = 2 R 1 0 *v dy*. Given a time-dependent pressure gradient Π(*t*), one can calculate the corresponding total flux *Q*(*t*). Let us consider the parametric line

$$
\Pi = \Pi\_0 (1 + \sin \pi t), \quad \mathbb{Q} = \mathbb{Q}(t), \quad 0 < t < 1,\tag{30}
$$

which corresponds to the curve *O*\, *P*, *L* on the (Π, *Q*)-plane, Figure 6. The lower part *O*\, *A*, *P* of this curve corresponds to the time interval 0 < *t* < 1/2. Along this part, both *Q* and Π grow, Π<sup>0</sup> < Π < 2Π0. The top part *P*\, *B*, *L* of the curve corresponds to the time interval 1/2 < *t* < 1. Along this part, both *Q* and Π decrease.

For typical viscous fluids like a power law fluid, there is a one-to-one correspondence between Π and *Q*; as a consequence, the lines *O*\, *P*, *L* and *P*\, *B*, *L* coincide. It is not the case for the anisotropic fluid considered here. Given Π<sup>∗</sup> satisfying the inequalities Π<sup>0</sup> < Π<sup>∗</sup> < 2Π0, how can one determine a corresponding flux *Q*? It follows from Figure 6 that there are two values *Q<sup>A</sup>* and *Q<sup>B</sup>* corresponding to Π∗. Indeed, let us consider the intersection of the vertical line <sup>Π</sup> <sup>=</sup> <sup>Π</sup><sup>∗</sup> with the curve *<sup>O</sup>*\, *<sup>P</sup>*, *<sup>L</sup>*. On this way we arrive at the points *A* and *B*:

$$A = (\Pi\_{\ast \prime} Q\_A)\_{\prime} \quad B = (\Pi\_{\ast \prime} Q\_B)\_{\cdot}$$

Clearly, there are *t<sup>A</sup>* and *t<sup>B</sup>* such that

$$10 < t\_A < 1/2 < t\_B < 1,\ II(t\_A) = \Pi(t\_B) = \Pi\_{\ast \ast} Q\_A < Q\_{B \prime} \ Q\_i = Q(t\_i)\_{\prime \ast}$$

with *i* = *A*, *B*.

Let us choose the points *C* = (Π*C*, *QC*) and *D* = (Π*D*, *QD*) in such a way that

$$
\Pi\_{\mathbb{C}} < \Pi\_\* < \Pi\_{D^\*}.
$$

When the value of Π goes from the low value Π*<sup>C</sup>* to Π∗, the value of *Q* changes from *Q<sup>C</sup>* to

$$Q\_A = \lim\_{\Pi \nearrow \Pi\_\*} Q(\Pi).$$

When the value of Π goes from the upper value Π*<sup>D</sup>* to Π∗, the value of *Q* changes from *Q<sup>D</sup>* to

$$Q\_B = \lim\_{\Pi \nearrow \Pi\_\*} Q(\Pi).$$

Thus, total flux depends not only on pressure gradient, but on the evolution history of pressure gradient as well.

**Figure 6.** On the (Π, *Q*)-plane, the hysteresis loop corresponding to process (30) for rather small *ε*20,*ε*30, with Π<sup>0</sup> = 2.3.

Similarly, we consider determination of Π starting from values of *Q*. Again, one should know a prehistory of *Q*. Indeed, let us consider a total flux *Q*∗, which is between *Q*|*t*=<sup>0</sup> and *Q*|*t*=1/2. Let us consider the intersection of the horizontal line *Q* = *Q*<sup>∗</sup> with the hysteresis loop *O*\, *P*, *L*. In this way, we arrive at the points *N* and *M*:

$$N = (\Pi\_{N\prime} Q\_\*)\_{\prime} \quad M = (\Pi\_{M\prime} Q\_\*)\_{\prime}.$$

Clearly, there are *t<sup>N</sup>* and *t<sup>M</sup>* such that

$$0 < t\_N < 1/2 < t\_M < 1,\\ \mathcal{Q}(t\_N) = \mathcal{Q}(t\_M) = \mathcal{Q}\_{\*\prime} \, \Pi\_M < \Pi\_{N\prime} \, \Pi\_i = \Pi(t\_i),$$

with *i* = *N*, *M*. Let us choose points *R* = (Π*R*, *QR*) and *S* = (Π*S*, *QS*) in such a way that *Q<sup>R</sup>* < *Q*<sup>∗</sup> < *QS*.

If *Q* goes from the lower value *Q<sup>R</sup>* to *Q*<sup>∗</sup> then Π changes from Π*<sup>R</sup>* to

$$\Pi\_N = \lim\_{Q \nearrow Q\_\*} \Pi(Q).$$

If *Q* goes from the upper value *Q<sup>S</sup>* to *Q*<sup>∗</sup> then Π changes from Π*<sup>S</sup>* to

$$\Pi\_M = \lim\_{Q \searrow Q\_\*} \Pi(Q).$$

Thus, pressure gradient depends not only on total flux, but on the prehistory evolution of total flux as well.

As far as the oil pipelines are concerned, the designed oil flux can be obtained in two ways: by switching from a fast or slow flux. By the developed anisotropic model, the pressure drop to ensure the designed oil flux is less in the first case.

Now, we consider friction losses which play an important role in oil pumping through pipelines. Returning to dimension variables, we remind that the mean velocity *U* and the friction factor *λ* are defined as follows:

$$\mathcal{U} = \frac{1}{2H} \int\_{-H}^{H} v(y) \, dy, \quad |p\_x| = \frac{\lambda}{2H} \frac{\rho \mathcal{U}^2}{2}.$$

In dimensionless variables, we have

$$\mathcal{U}' = \int\_0^1 v'(y) dy = \frac{\mathcal{Q}}{2}, \quad \Lambda \equiv \frac{\text{Re} \cdot \lambda}{8} = \frac{\Pi}{\mathcal{U}'^2}.$$

where Λ is the reduced friction factor.

To analyse flows on the plane (*U*0 , Λ), we omit the prime indexes. Calculations reveal that, starting from the pressure gradient law

$$
\Pi(t) = \Pi\_0 (1 + \sin \pi t),
$$

the curve

$$\mathcal{U} = \mathcal{U}(t), \quad \Lambda(t) = \frac{\Pi(t)}{\mathcal{U}^2(t)}, \quad 0 < t < 1,$$

becomes as is shown in Figure 7. The top part of this curve corresponds to the time interval 0 < *t* < 1/2. Along this part, both *U* and Π grow, Π<sup>0</sup> < Π < 2Π0, whereas Λ decreases. The lower part of the curve corresponds to the time interval 1/2 < *t* < 1. Along this part, both *U* and Π decrease, whereas Λ grows. How can one calculate the friction factor Λ corresponding to a designed mean velocity *U*∗ ? The answer depends on the history; one can attain *U*<sup>∗</sup> by increasing *U* or by decreasing *U*. Given *U*<sup>∗</sup> lying between *Umin* = *U*|*t*=<sup>0</sup> and *Umax* = *U*|*t*=1/2, we choose *t*<sup>1</sup> and *t*<sup>2</sup> in such a away that

$$0 < t\_1 < 1/2 < t\_2 < 1, \quad \mathcal{U}(t\_1) = \mathcal{U}(t\_2) = \mathcal{U}\_\*.$$

With Λ*<sup>i</sup>* standing for Λ(*ti*), one can conclude from Figure 7 that Λ<sup>1</sup> > Λ<sup>2</sup> despite the fact that both Λ<sup>1</sup> and Λ<sup>2</sup> correspond to the same *U*∗. Thus,

$$\Lambda\_1 = \lim\_{\substack{\mathcal{U} \nearrow \mathcal{U}\_\bullet}} \Lambda(\mathcal{U}) = \Lambda|\_{\mathcal{U}\_\bullet -} > \Lambda|\_{\mathcal{U}\_\bullet +} = \lim\_{\substack{\mathcal{U} \searrow \mathcal{U}\_\bullet}} \Lambda(\mathcal{U}) = \Lambda\_2.$$

Returning to the issue of oil transportation, one can attain the productive regime in two ways by switching from a faster or from a slower flux. After switching to a productive regime, the developed friction loss is less in the first case. Such a conclusion agrees with known in situ data [9].

Consider the stress response to a change in velocity gradient. It follows from the dimensional steady-state Equation (22) that the shear stress *S*<sup>12</sup> is given by the formula

$$S\_{12} = p\_{\underline{x}} y\_{\prime} \quad \tilde{\tau} \equiv -\left. \mathcal{S}\_{12} \right|\_{y=H} = -p\_{\underline{x}} H\_{\prime 1}$$

where *τ*˜ is the stress at the upper plane *y* = *H*. Let us calculate the curve *τ*˜ = *τ*˜(*γ*˙ <sup>1</sup>) where *γ*˙ <sup>1</sup> = −*v<sup>y</sup> y*=*H* . Observe that *γ*˙ <sup>1</sup> does not stay for the the shear rate in the micropolar fluid theory. We pass to the dimensionless variables

$$\tau = \frac{\mathfrak{k}}{2\mu\_s \Omega} = \Pi, \quad \dot{\gamma} = -v\_{y'}'\Big|\_{y'=1} \quad \dot{\gamma}\_1 = \Omega \dot{\gamma}.$$

Performing calculation of the parametric curve

$$
\Pi = \Pi\_0 (1 + \sin \pi t), \quad \dot{\gamma} = \dot{\gamma}(t), \quad 0 < t < 1,\tag{31}
$$

we arrive at the hysteresis loop *τ* = *τ*(*γ*˙), which is shown in Figure 8. Thus, there is no one-to-one correspondence between velocity gradient and shear stress in shear flows. Such an effect has been seen in worm-like micelle solutions [39].

**Figure 7.** Hysteresis loop on the (*U*, Λ) plane, where *U* is the mean velocity and Λ is the friction factor. The data are the same as in Figure 2.

**Figure 8.** Hysteresis of the rheological stress–strain curve Π = Π(*γ*˙).

Let us introduce the apparent viscosity *η<sup>a</sup>* = *τ*˜/*γ*˙ <sup>1</sup>. Figure 9 depicts how its dimensionless replica *η* = *τ*/*γ*˙ varies with time for the case Π = const. Sustained temporal oscillations of macroscopic viscosity are observed in [40] for the rod-like suspension.

**Figure 9.** Apparent viscosity versus time. The case Π = Π<sup>0</sup> = const.

As far as the function *η* = *η*(*γ*˙) is concerned, the shear thinning nature of suspensions of rod-like particles is clearly depicted on Figure 10 in agreement with observations in [1].

**Figure 10.** Apparent viscosity *η* versus velocity gradient *γ*˙ for the case Π = Π<sup>0</sup> = const.

#### **5. Discussion**

We address rheology of suspensions of rodlike particles. To take into account both particle–fluid and particle–particle interactions, we treat the suspension as a Cosserat continuum and apply the micropolar fluid theory approach. Assuming that local stress depends on the rods directions, we include the micro-inertia tensor into constitutive laws as an independent variable jointly with the rate of strain tensors. The micropolar fluid theory allows for particle's rotation obeying the angular momentum conservation law. The Cauchy stress tensor loses symmetry and the couple stress tensor is of importance. This is why one should formulate two stress-rate of strain rheological equations for the stress tensor and the couple stress tensor. Unlike a Newtonian fluid, a micropolar fluid is characterized by two rates of strain tensors, through which the linear velocity gradient and the angular velocity gradient are expressed. The impact of variation of rate of strains and the microinertia onto the local stress state in the rheological equations is manifested through the viscosities. This is why, in addition to the usual shear viscosity, we also introduce skewsymmetric and anisotropic viscosities. The derived governing equations are proved to be consistent with basic principles of thermodynamics. By performing calculations of simple one-dimensional pressure driven flows between two parallel planes, we establish that the skew-symmetric and the anisotropic viscosities underlie a number of important properties, which include gradient banding instability, coexistence of isotropic and nematic phases, sustained temporal oscillations of macroscopic viscosity, shear thinning and hysteresis. Keeping in mind data for oil transport in pipelines, we also establish that the total flow rate depends not only on the pressure gradient, but on the history of its variation as well.

**Funding:** The research in Section 2 on rheology of anisotropic micropolar fluids is funded by the Government of the Russian Federation (Grant No. 14.W03.31.0002). The theoretical research and computations in Sections 3–5 concerning flows between two planes are funded by Russian Science Foundation (Grant No. 20-19-00058; Funder ID: 10.13039/501100006769).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

