Flows of Linear Polymer Solutions and Other Suspensions of Rod-like Particles: Anisotropic Micropolar-Fluid Theory Approach

Abstract

We formulate equations governing flows of suspensions of rod-like particles. Such suspensions include linear polymer solutions, FD-virus, and worm-like micelles. To take into account the particles that form and their rotation, we treat the suspension as a Cosserat continuum and apply the theory of micropolar fluids. Anisotropy of suspensions is determined through the inclusion of the microinertia tensor in the rheological constitutive equations. We check that the model is consistent with the basic principles of thermodynamics. In addition to anisotropy, the theory also captures gradient banding instability, coexistence of 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.

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,3,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 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 $\mathit{\xi }$, the position vector $\mathbf{x}(t,\mathit{\xi })$ at the time instant t in the three dimensional Euclidean space jointly with orthogonal director-vectors ${\mathbf{d}}_i(t,\mathit{\xi })$, $i=1,2,3$, are assigned to such a point to treat it as a rigid body. Orientation of ${\mathbf{d}}_i$ is controlled by an orthogonal tensor $\mathrm{Q}(t,\mathit{\xi })$:

$${\mathbf{d}}_i(t,\mathit{\xi })=\mathrm{Q}(t,\mathit{\xi }){\mathbf{d}}_i(0,\mathit{\xi }),{\mathrm{Q}}^{\ast }\mathrm{Q}=\mathrm{Q}{\mathrm{Q}}^{\ast }=I.$$

Here, I is the identity matrix with the elements ${\delta }_j^i$, ${\mathrm{Q}}^{\ast }$ is the adjoint matrix, ${\left({\mathrm{Q}}^{\ast }\right)}_{ij}=$ ${\mathrm{Q}}_{ji}$. The skew-symmetric tensor $\Omega (t,\mathbf{x})={\mathrm{Q}}_t{\mathrm{Q}}^{\ast }$ defines the particle’s rotation with the angular velocity

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

where ${\left\{{\mathbf{e}}_i\right\}}_1^3$ is an orthonormal basis in ${\mathbb{R}}^3$ and $ϵ$ is the Levi–Civita third order tensor,

$$\epsilon (\mathbf{a},\mathbf{b},\mathbf{c})=\mathbf{a}·(\mathbf{b}\times \mathbf{c}),{\mathbf{e}}_i\times {\mathbf{e}}_j={ϵ}_{sij}{\mathbf{e}}_s,{ϵ}_{sij}=\epsilon ({\mathbf{e}}_s,{\mathbf{e}}_i,{\mathbf{e}}_j),{(ϵ:\Omega )}_i={ϵ}_{ijk}{\Omega }_{jk}.$$

Given the velocity field $\mathbf{v}(t,\mathbf{x})=$ ${\mathbf{x}}_t(t,\mathit{\xi })$, we introduce the rate of strain tensors [22].

$$B=\nabla \mathbf{v}-\Omega ,A=\nabla \mathit{\omega },$$

where ${(\nabla \mathbf{v})}_{ij}=\partial v_i/\partial x_j$. Observe that both B and A are objective relative to the change of frame of references.

Let $\rho$, $\mathbf{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 $\dot{\rho }$ (or $d\rho /dt$) related to the velocity field $\mathbf{v}$ as follows

$$\dot{\rho }=\frac{\partial \rho }{\partial t}+v_i\frac{\partial \rho }{\partial x_i}\mathrm{or}\dot{\rho }=\frac{\partial \rho }{\partial t}+(\mathbf{v}·\nabla )\rho .$$

(1)

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

$$\dot{\rho }+\rho \mathrm{d}\mathrm{i}\mathrm{v}\mathbf{v}=0,$$

(2)

$$\rho \dot{\mathbf{v}}=\mathrm{d}\mathrm{i}\mathrm{v}T+\rho \mathbf{f},$$

(3)

$$\rho \dot{\mathbf{s}}=\mathrm{d}\mathrm{i}\mathrm{v}N-\mathit{ϵ}:T+\rho \mathit{l},$$

(4)

where $\mathbf{f}$ is the mass force density, $\mathit{l}$ is the mechanical couple density, and

$${\left(\mathrm{d}\mathrm{i}\mathrm{v}T\right)}_i=\partial T_{ij}/\partial x_j.$$

Observe that the stress tensor T is not symmetric. Given an orthonormal basis ${\left\{{\mathbf{e}}_i\right\}}_1^3$, the vector

$$\mathbf{t}={\mathbf{e}}_i\times (\mathrm{T}·{\mathbf{e}}_i)=\mathit{\epsilon }:\mathrm{T}$$

does not depend on the choice of the basis and it is a stress symmetry defect measure in the sense that the equality $\mathbf{t}=0$ implies ${\mathrm{T}}^{\ast }=\mathrm{T}$ and vice versa. By the definition of $\mathbf{t}$, we have the formula

$$\mathbf{t}·\mathit{\omega }=\mathrm{T}:\Omega .$$

(5)

The internal specific spin is defined by the formula $\mathbf{s}=$ $J\mathit{\omega }$, where the symmetric inertia tensor $J\left[{\mathrm{cm}}^2\right]$ obeys the identity [22]

$$\dot{J}-\Omega J+J\Omega =0.$$

(6)

Before proceeding to constitutive laws, we address the thermodynamics issue. Given a specific internal energy e, the total energy $E=\mathrm{e}+\mathbf{v}·\mathbf{v}/2+\mathbf{s}·\mathit{\omega }/2$ satisfies the equation [32]

$$\rho \dot{E}=\mathrm{d}\mathrm{i}\mathrm{v}(T^{\ast }\mathbf{v}+N^{\ast }\mathit{\omega }-\mathbf{q})+\rho \mathbf{f}·\mathbf{v}+\rho \mathit{l}·\mathit{\omega },$$

(7)

where $\mathbf{q}$ is the heat flux obeying the Fourier law $\mathbf{q}=-$ $\varkappa \nabla \theta$, with $\varkappa$ standing for the heat conductivity. Generally, internal energy e depends on $\rho$, $\eta$ and J, $e=e(\rho ,\eta ,J)$, where $\eta$ is the specific entropy. It is common knowledge that absolute temperature and pressure are defined by the derivatives $\theta =\partial e/\partial \eta$, $p={\rho }^2\partial e/\partial \rho$ respectively [33]. We calculate that

$$\dot{\mathrm{e}}={\mathrm{e}}_{\rho }\dot{\rho }+{\mathrm{e}}_{\eta }\dot{\eta }+{\nabla }_{\mathrm{J}}\mathrm{e}:\dot{J},\mathrm{where}{\nabla }_{\mathrm{J}}\mathrm{e}:\dot{J}=\left({\dot{J}}_{ij}\frac{\partial }{\partial J_{ij}}\right)\mathrm{e}.$$

From the rheological point of view, the internal energy $e(\rho ,\eta ,J)$ should be an isotropic function of J. Hence, ${\nabla }_{\mathrm{J}}\mathrm{e}$ is also an isotropic function of J; it implies that [34]

$${\nabla }_{\mathrm{J}}\mathrm{e}={\alpha }_{0}I+{\alpha }_{1}J+{\alpha }_2{J}^2,$$

(8)

where the scalar functions ${\alpha }_i$ depend on invariants of J. Now, it follows from (6) and (8) that ${\nabla }_{\mathrm{J}}\mathrm{e}:\dot{J}=0.$

We use Equation (6) to calculate that

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

Hence,

$$\rho \dot{E}=\frac{p}{\rho }\dot{\rho }+\rho \theta \dot{\eta }+\rho \dot{\mathbf{v}}·\mathbf{v}+\rho \dot{\mathbf{s}}·\mathit{\omega }=-p\mathrm{div}\mathbf{v}+\rho \dot{\mathbf{v}}·\mathbf{v}+\rho \dot{\mathbf{s}}·\mathit{\omega }.$$

(9)

Multiplying Equations (3) and (4) by $\mathbf{v}$ and $\mathit{\omega }$, respectively, we arrive at the energetic equality

$$\rho \mathbf{v}·\dot{\mathbf{v}}+\rho \mathit{\omega }·\dot{\mathbf{s}}=\mathrm{div}(T^{\ast }\mathbf{v}+N^{\ast }\mathit{\omega })-T:B-N:A+\rho \mathbf{f}·\mathbf{v}+\rho \mathit{l}·\mathit{\omega }$$

(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\mathrm{div}\mathbf{v}+$ $S:B$, where $T:B=T_{ij}{B}_{ij}$. Now, it follows from (7), (9) and (10) that the entropy equation

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

(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\ge 0$ can be formulated as

$$S:B+N:A\ge 0.$$

(12)

In what follows, we use the notations

$$B_s=\frac{B+B^{\ast }}{2},B_a=\frac{B-B^{\ast }}{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}JB,N=\frac{2\nu }{{\sigma }^2}A+2{\nu }_{an}AJ,$$

(13)

where ${\mu }_s,{\mu }_a,{\mu }_{an},\nu ,{\nu }_{an}\left[\mathrm{g}/(\mathrm{cm}·\mathrm{s})\right]$ are the viscosities and $\sigma \left[{\mathrm{cm}}^{-1}\right]$ is the specific particles surface area. The first rheological equation in (13) suggests that the contributions of the symmetric part $B_s$ and skew-symmetric part $B_a$ 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{\left|B^{\ast }{\mathbf{e}}_j\right|}^2,$$

where ${\mathbf{e}}_j$ and ${\lambda }_j$ are the eigenvectors and the eigenvalues of J. Observe that ${\lambda }_j\ge 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{\left|B^{\ast }{\mathbf{e}}_j\right|}^2\ge 0.$$

Similarly, one can verify that

$$AJ:J=\sum {\lambda }_j{\left|A{\mathbf{e}}_j\right|}^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 $\left|y\right|<H$ between two parallel planes under the prescribed pressure gradient $\nabla p=$ $(p_x,0,0)$, $p_x\left(t\right)<0$, Figure 1a. In this case, $v_2=v_3=0$, $v_1=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, $\mathit{\omega }=(0,0,\omega )$.

Let us describe the micro-inertia tensor J. First, we consider a cylinder $V_0$ stretched along the y-axis with the height h and the radius r, Figure 1b. By definition, the inertia tensor $J\left(V_0\right)$ of $V_0$ is given by the formula

$$J\left(V_0\right)=\underset{V_0}{\int }{\left|\mathit{\zeta }\right|}^2·I-\mathit{\xi }\otimes \mathit{\xi }d\xi \mathrm{or}{J}_{ij}\left(V_0\right)=\underset{V_0}{\int }{\left|{\xi }_k\right|}^2{\delta }_{ij}-{\xi }_i{\xi }_{j}d\xi ,$$

where I is the identity matrix and $\mathbf{a}\otimes \mathbf{b}$ stands for the tensor product of two vectors $\mathbf{a}$ and $\mathbf{b}$, ${(\mathbf{a}\otimes \mathbf{b})}_{ij}=$ $a_i{b}_j$. Calculations reveal that

$$J\left(V_0\right)=\left(\begin{array}{ccc}{r}^2/4+h^2/3& 0& 0\\ 0& r^2/2& 0\\ 0& 0& r^2/4+h^2/3\end{array}\right).$$

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

$$J^0=\underset{r\to 0}{lim}J\left(V_0\right)=j_0\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right).j_0=h^2/3.$$

Let $J\left(V\right)$ be the inertia tensor of the cylinder V, which is produced by rotation of $V_0$ around the z-axis by the angle $\phi$ counted from the axis y counter-clockwise, see Figure 1a. By definition of the spin $\mathbf{s}$, we find that

$$\mathbf{s}=J\left(V\right)·\mathit{\omega }=\underset{V}{\int }\mathbf{x}\times (\mathit{\omega }\times \mathbf{x})dx=Q_{\phi }{J}^0{Q}_{\phi }^{\ast }·\mathit{\omega },$$

where $Q_{\phi }$ is the orthogonal matrix such that

$$\Omega ={\dot{Q}}_{\phi }{Q}_{\phi }^{\ast },\Omega ·\mathbf{h}=\mathit{\omega }\times \mathbf{h}\forall \mathbf{h},{\left(Q^{\ast }\right)}_{ij}=Q_{ji},$$

(14)

$$Q_{\phi }=\left(\begin{array}{ccc}cos\phi & -sin\phi & 0\\ sin\phi & cos\phi & 0\\ 0& 0& 1\end{array}\right),\Omega =\left(\begin{array}{ccc}0& -\omega & 0\\ \omega & 0& 0\\ 0& 0& 0\end{array}\right),$$

(15)

with “dot” standing for the material derivative (1) related to the velocity vector $\mathbf{v}$. Thus, $J\left(V\right)=Q_{\phi }J\left(V_0\right)Q_{\phi }^{\ast }$. In the limit as $r\to 0$, we obtain the inertia tensor $J\left(\phi \right)$ of the rod particle with the position angle $\phi$:

$$J\left(\phi \right)=Q_{\phi }{J}^0{Q}_{\phi }^{\ast }=j_0\left(\begin{array}{ccc}{cos}^2\phi & sin\phi cos\phi & 0\\ sin\phi cos\phi & {sin}^2\phi & 0\\ 0& 0& 1\end{array}\right),\frac{\partial \phi }{\partial t}={\phi }_t=\omega .$$

(16)

Given an initial distribution of particle’s angles ${\phi }_0\left(y\right)$, we denote the initial micro-inertia tensor by $J_0\left(y\right)=$ $J\left({\phi }_0\left(y\right)\right)$:

$${J|}_{t=0}=J_0.$$

(17)

For the described one-dimensional flows, the material derivative $\dot{J}$ reduces to the time derivative $J_t$. With the use of Equation (15), (6) can be written as follows:

$$\frac{\partial }{\partial t}{J}_{11}=-2\omega J_{12},\frac{\partial }{\partial t}{J}_{12}=\omega (J_{11}-J_{22}),\frac{\partial }{\partial t}{J}_{22}=2\omega J_{12},$$

(18)

and $J_{ij}=0$ otherwise. Observe that

$$\frac{\partial J_{ij}}{\partial t}=\frac{\partial J_{ij}}{\partial \phi }\frac{\partial \phi }{\partial t}=J_{ij}^{\prime }\omega ,J_{ij}^{\prime }=\frac{\partial J_{ij}}{\partial \phi }.$$

Hence, the system (18) is equivalent to

$$J_{11}^{\prime }=-2J_{12},J_{12}^{\prime }=J_{11}-J_{22},J_{22}^{\prime }=2J_{12},J_{ij}{|}_{\phi ={\phi }_0}=J_0.$$

(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=\left(\begin{array}{ccc}0& v_y+\omega & 0\\ -\omega & 0& 0\\ 0& 0& 0\end{array}\right),B^{\ast }=\left(\begin{array}{ccc}0& -\omega & 0\\ v_y+\omega & 0& 0\\ 0& 0& 0\end{array}\right),$$

(20)

$$B_s=\left(\begin{array}{ccc}0& v_y/2& 0\\ v_y/2& 0& 0\\ 0& 0& 0\end{array}\right),B_a=\left(\begin{array}{ccc}0& v_y/2+\omega & 0\\ -v_y/2-\omega & 0& 0\\ 0& 0& 0\end{array}\right),$$

$$j_0^{-1}JB=\left(\begin{array}{ccc}-\omega cos\phi sin\phi & {cos}^2\phi (v_y+\omega )& 0\\ -\omega {sin}^2\phi & cos\phi sin\phi (v_y+\omega )& 0\\ 0& 0& 0\end{array}\right),$$

$$A=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& {\omega }_y& 0\end{array}\right),j_0^{-1}AJ=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ {\omega }_{y}cos\phi sin\phi & {\omega }_y{sin}^2\phi & 0\end{array}\right),$$

(21)

Let us denote

$${\epsilon }_1=\frac{{\mu }_a}{{\mu }_s},{\epsilon }_{20}=\frac{{\mu }_{an}{j}_0{\sigma }^2}{{\mu }_s},{\epsilon }_{30}=\frac{{\nu }_{an}{j}_0{\sigma }^2}{\nu },B^0=B_s+{\epsilon }_1{B}_a+{\epsilon }_{20}JB.$$

Calculations reveal that matrix $B^0$ is equal to

$$\left(\begin{array}{ccc}-{\epsilon }_{20}\omega cos\phi sin\phi & \frac{v_y(1+{\epsilon }_1+2{\epsilon }_{20}{cos}^2\phi )}{2}+\omega ({\epsilon }_1+{\epsilon }_{20}{cos}^2\phi )& 0\\ \frac{v_y(1-{\epsilon }_1)}{2}-\omega ({\epsilon }_1+{\epsilon }_{20}{sin}^2\phi )& {\epsilon }_{20}cos\phi sin\phi (v_y+\omega )& 0\\ 0& 0& 0\end{array}\right).$$

We consider incompressible fluids with the assumption $\rho =$ const. For one-dimensional flows, the incompressibility condition $\mathrm{div}\mathbf{v}=0$ is satisfied automatically. Other conservation laws (3) and (4) become

$${\phi }_t=\omega ,\rho v_t=-p_x+\frac{\partial S_{12}}{\partial y},\rho j_0{\omega }_t=\frac{\partial N_{32}}{\partial y}+S_{21}-S_{12}.$$

(22)

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

$$S_{ij}=2{\mu }_s{B}_{ij}^0,N_{32}=2\frac{\nu }{{\sigma }^2}{A}_{32}+2{\nu }_{an}{\left(AJ\right)}_{32}.$$

(23)

Observe that

$$S_{21}-S_{12}=2{\mu }_s(B_{21}^0-B_{12}^0),B_{21}^0-B_{12}^0=-v_y({\epsilon }_1+{\epsilon }_{20}{cos}^2\phi )-\omega (2{\epsilon }_1+{\epsilon }_{20}).$$

We formulate boundary conditions at $\left|y\right|=H$ as follows:

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

(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 $\mathit{\omega }$ depends linearly on the macro-rotation $\nabla \times \mathbf{v}/2$ [25].

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

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

$$y^{\prime }=\frac{y}{H},v^{\prime }=\frac{v}{V},t^{\prime }=\frac{t}{T},{\omega }^{\prime }=\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

$${\phi }_{t^{\prime }}={\omega }^{\prime },\frac{Re}{2}{v}_{t^{\prime }}^{\prime }=\Pi +\frac{\partial S_{12}^{\prime }}{\partial y^{\prime }},\frac{Re{j}_0}{2H^2}{\omega }_{t^{\prime }}^{\prime }=\gamma \frac{\partial N_{32}^{\prime }}{\partial y}+S_{21}^{\prime }-S_{12}^{\prime }.$$

(25)

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

$${\phi }_{t^{\prime }}={\omega }^{\prime },0=\Pi +\frac{\partial S_{12}^{\prime }}{\partial y^{\prime }},$$

(26)

$$0=\gamma \frac{\partial }{\partial y^{\prime }}\left[{\omega }_{y^{\prime }}^{\prime }(1+{\epsilon }_{30}{sin}^2\phi )\right]-[v_{y^{\prime }}^{\prime }({\epsilon }_1+{\epsilon }_{20}{cos}^2\phi )+{\omega }^{\prime }(2{\epsilon }_1+{\epsilon }_{20})],$$

(27)

where

$$S_{12}^{\prime }=\frac{v_{y^{\prime }}^{\prime }(1+{\epsilon }_1+2{\epsilon }_{20}{cos}^2\phi )}{2}+{\omega }^{\prime }({\epsilon }_1+{\epsilon }_{20}{cos}^2\phi ).$$

Because of the symmetry conditions

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

we consider flows only in the upper half-layer $0<y^{\prime }<1$. In such a case the initial and boundary conditions take the form

$${\phi |}_{t^{\prime }=0}={\phi }_0\left(y^{\prime }\right),v^{\prime }\left(1\right)=0,v_{y^{\prime }}^{\prime }\left(0\right)=0,{\omega }^{\prime }\left(1\right)=-0.5\alpha v_{y^{\prime }}^{\prime }\left(1\right),{\omega }^{\prime }\left(0\right)=0.$$

(28)

To perform numerical solution, one should fix the dimensionless parameters $\Pi$, ${\epsilon }_1$, ${\epsilon }_{20}$, ${\epsilon }_{30}$, $\gamma$, $\alpha$.

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,3,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 $\mathbf{v}=$ $(v,0,0)$, y is the gradient direction along which the flow has non-zero derivative $\partial v/\partial y$. The z-axis is the vorticity direction along the non-zero macro-vorticity vector $\nabla \times \mathbf{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

$${\epsilon }_1=1,{\epsilon }_{20}=2,{\epsilon }_{30}=2,\gamma =1.3,\alpha =0.3,{\phi }_0=0.$$

(29)

Intervals where $\phi \left(y\right)=$ const or $\omega \left(y\right)=$ const correspond to the nematic phase. The profiles of the intrinsic angular velocity $\omega$ at Figure 2b and Figure 3 imply appearance and instability of the nematic phase. Figure 4b depicts the phase transition from the isotropic phase to the nematic phase.

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 ${\phi }_0\left(y\right)=0$ to a spatially heterogeneous orientation (like ${\phi }_0\left(y\right)=$ $4y+$ $9y^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).

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 $\Pi$ and the total fluid flux $Q=2{\int }_0^{1}vdy$. Given a time-dependent pressure gradient $\Pi \left(t\right)$, one can calculate the corresponding total flux $Q\left(t\right)$. Let us consider the parametric line

$$\Pi ={\Pi }_0(1+sin\pi t),Q=Q\left(t\right),0<t<1,$$

(30)

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

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

$$A=({\Pi }_{\ast },Q_A),B=({\Pi }_{\ast },Q_B).$$

Clearly, there are $t_A$ and $t_B$ such that

$$0<t_A<1/2<t_B<1,\Pi \left(t_A\right)=\Pi \left(t_B\right)={\Pi }_{\ast },Q_A<Q_B,Q_i=Q\left(t_i\right),$$

with $i=A,B.$

Let us choose the points $C=({\Pi }_C,Q_C)$ and $D=({\Pi }_D,Q_D)$ in such a way that

$${\Pi }_C<{\Pi }_{\ast }<{\Pi }_D.$$

When the value of $\Pi$ goes from the low value ${\Pi }_C$ to ${\Pi }_{\ast }$, the value of Q changes from $Q_C$ to

$$Q_A=\underset{\Pi \nearrow {\Pi }_{\ast }}{lim}Q(\Pi ).$$

When the value of $\Pi$ goes from the upper value ${\Pi }_D$ to ${\Pi }_{\ast }$, the value of Q changes from $Q_D$ to

$$Q_B=\underset{\Pi \searrow {\Pi }_{\ast }}{lim}Q(\Pi ).$$

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

Similarly, we consider determination of $\Pi$ starting from values of Q. Again, one should know a prehistory of Q. Indeed, let us consider a total flux $Q_{\ast }$, which is between ${Q|}_{t=0}$ and ${Q|}_{t=1/2}$. Let us consider the intersection of the horizontal line $Q=Q_{\ast }$ with the hysteresis loop $\widehat{O,P,L}$. In this way, we arrive at the points N and M:

$$N=({\Pi }_N,Q_{\ast }),M=({\Pi }_M,Q_{\ast }).$$

Clearly, there are $t_N$ and $t_M$ such that

$$0<t_N<1/2<t_M<1,Q\left(t_N\right)=Q\left(t_M\right)=Q_{\ast },{\Pi }_M<{\Pi }_N,{\Pi }_i=\Pi \left(t_i\right),$$

with $i=N,M.$ Let us choose points $R=({\Pi }_R,Q_R)$ and $S=({\Pi }_S,Q_S)$ in such a way that $Q_R<Q_{\ast }<Q_S.$

If Q goes from the lower value $Q_R$ to $Q_{\ast }$ then $\Pi$ changes from ${\Pi }_R$ to

$${\Pi }_N=\underset{Q\nearrow Q_{\ast }}{lim}\Pi \left(Q\right).$$

If Q goes from the upper value $Q_S$ to $Q_{\ast }$ then $\Pi$ changes from ${\Pi }_S$ to

$${\Pi }_M=\underset{Q\searrow Q_{\ast }}{lim}\Pi \left(Q\right).$$

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 $\lambda$ are defined as follows:

$$U=\frac{1}{2H}\underset{-H}{\overset{H}{\int }}v\left(y\right)dy,\left|p_x\right|=\frac{\lambda }{2H}\frac{\rho U^2}{2}.$$

In dimensionless variables, we have

$$U^{\prime }={\int }_0^1{v}^{\prime }\left(y\right)dy=\frac{Q}{2},\Lambda \equiv \frac{\mathrm{R}\mathrm{e}·\lambda }{8}=\frac{\Pi }{U^{\prime }{}^2},$$

where $\Lambda$ is the reduced friction factor.

To analyse flows on the plane $(U^{\prime },\Lambda )$, we omit the prime indexes. Calculations reveal that, starting from the pressure gradient law

$$\Pi \left(t\right)={\Pi }_0(1+sin\pi t),$$

the curve

$$U=U\left(t\right),\Lambda \left(t\right)=\frac{\Pi \left(t\right)}{U^2\left(t\right)},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 $\Pi$ grow, ${\Pi }_0<\Pi <2{\Pi }_0$, whereas $\Lambda$ decreases. The lower part of the curve corresponds to the time interval $1/2<t<1$. Along this part, both U and $\Pi$ decrease, whereas $\Lambda$ grows. How can one calculate the friction factor $\Lambda$ corresponding to a designed mean velocity $U^{\ast }$? The answer depends on the history; one can attain $U_{\ast }$ by increasing U or by decreasing U. Given $U_{\ast }$ lying between $U_{min}{=U|}_{t=0}$ and $U_{max}{=U|}_{t=1/2}$, we choose $t_1$ and $t_2$ in such a away that

$$0<t_1<1/2<t_2<1,U\left(t_1\right)=U\left(t_2\right)=U_{\ast }.$$

With ${\Lambda }_i$ standing for $\Lambda \left(t_i\right)$, one can conclude from Figure 7 that ${\Lambda }_1>{\Lambda }_2$ despite the fact that both ${\Lambda }_1$ and ${\Lambda }_2$ correspond to the same $U_{\ast }$. Thus,

$${\Lambda }_1=\underset{U\nearrow U_{\ast }}{lim}\Lambda \left(U\right)=\Lambda {{|}_{U_{\ast }-}>\Lambda |}_{U_{\ast }+}=\underset{U\searrow U_{\ast }}{lim}\Lambda \left(U\right)={\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_{12}$ is given by the formula

$$S_{12}=p_{x}y,\tilde{\tau }\equiv -S_{12}{|}_{y=H}=-p_{x}H,$$

where $\tilde{\tau }$ is the stress at the upper plane $y=H$. Let us calculate the curve $\tilde{\tau }=\tilde{\tau }\left({\dot{\gamma }}_1\right)$ where ${\dot{\gamma }}_1=-v_y{|}_{y=H}$. Observe that ${\dot{\gamma }}_1$ does not stay for the the shear rate in the micropolar fluid theory. We pass to the dimensionless variables

$$\tau =\frac{\tilde{\tau }}{2{\mu }_s\Omega }=\Pi ,\dot{\gamma }=-v_{y^{\prime }}^{\prime }{|}_{y^{\prime }=1},{\dot{\gamma }}_1=\Omega \dot{\gamma }.$$

Performing calculation of the parametric curve

$$\Pi ={\Pi }_0(1+sin\pi t),\dot{\gamma }=\dot{\gamma }\left(t\right),0<t<1,$$

(31)

we arrive at the hysteresis loop $\tau =\tau \left(\dot{\gamma }\right)$, 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].

Let us introduce the apparent viscosity ${\eta }_a=\tilde{\tau }/{\dot{\gamma }}_1$. Figure 9 depicts how its dimensionless replica $\eta =\tau /\dot{\gamma }$ varies with time for the case $\Pi =\mathrm{const}$. Sustained temporal oscillations of macroscopic viscosity are observed in [40] for the rod-like suspension.

As far as the function $\eta =\eta \left(\dot{\gamma }\right)$ is concerned, the shear thinning nature of suspensions of rod-like particles is clearly depicted on Figure 10 in agreement with observations in [1].

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 micro-inertia 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 skew-symmetric 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.