*Article* **Rotational Particle Separation in Solutions: Micropolar Fluid Theory Approach**

**Vladimir Shelukhin 1,2**

<sup>1</sup> Lavrentyev Institute of Hydrodynamics, 630090 Novosibirsk, Russia; shelukhin@list.ru

<sup>2</sup> Novosibirsk State University, 630090 Novosibirsk, Russia

**Abstract:** We develop a new mathematical model for rotational sedimentation of particles for steady flows of a viscoplastic granular fluid in a concentric-cylinder Couette geometry when rotation of the Couette cell inner cylinder is prescribed. We treat the suspension as a micro-polar fluid. The model is validated by comparison with known data of measurement. Within the proposed theory, we prove that sedimentation occurs due to particles' rotation and rotational diffusion.

**Keywords:** suspensions; micro-polar fluids; yield stress

### **1. Introduction**

The classical water-based drilling muds contain only water and clay and their performances are generally poor. Polymers used in drilling fluids improve stability and cutting removal. Currently, the polymer-based drilling fluids represent 15 to 18% of the total cost (about 1 million) of petroleum well drilling [1]. The reason is that such fluids appear to have load carrying capabilities, or, in other words, a yield stress, associated with the solid-like state and which primarily arises from the colloidal forces between the smallest suspended particles. Furthermore, in such systems, when the agitation is increased, a fluid-like state is recovered. Even in the fluid state, these materials usually show highly nonlinear behavior. The complex rheology related to the change of the internal state in the suspensions is still rather poorly understood.

Due to scale separation between colloidal and non-colloidal particles, polymer-based drilling muds with cuttings and other materials (concrete casting, foodstuff transport, etc.) can be considered as suspensions of noncolloidal particles embedded in a yield stress fluid. Substantial progress in the understanding of the behavior of such materials can thus be made by studying the impact of adding noncolloidal particles to a yield stress fluid of known properties [2]. Here, we develop a new mathematical model for suspensions embedded in a polymer-based fluid. To validate the approach, we consider rotation flows between two concentric cylinders when the external cylinder is fixed and the inner one rotates with a prescribed speed. The principal feature of such flows is the particles' migration toward the external boundary.

It is proved by Svedberg [3] that the particles' migration effect occurs also in pure colloidal suspensions. Moreover, sedimentation happens at high rotations. Currently, Svedberg's ultracentrifugation is known as an effective tool for studies of interaction between macromolecules of colloidal systems [4]. In this method, a highly disperse colloidal solution is enclosed in a wedge-shape cell rotating about an axis coinciding with the apex of the wedge [5]. Samples are centrifuged at speeds to produce sedimentation and shallow concentration gradients. A great progress in the study of such a flow was achieved by applying the diffusive Lamm equation based on the special empirical sedimentation coefficient [6]. Here, we restrict ourself to rotation flows of suspensions between two concentric cylinders paying attention to comparison with known laboratory experiments. We use methods of mechanics of a continuum by applying conservation laws only and

**Citation:** Shelukhin, V. Rotational Particle Separation in Solutions: Micropolar Fluid Theory Approach. *Polymers* **2021**, *13*, 1072. https:// doi.org/10.3390/polym13071072

Academic Editor: Célio Pinto Fernandes

Received: 17 February 2021 Accepted: 24 March 2021 Published: 29 March 2021

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

**Copyright:** © 2021 by the authors. 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/).

not involving the concept of sedimentation coefficient inherent in the Lamm equation. However, we also obtain the sedimentation effect.

We take into account particles rotation and rotational diffusion. To this end, we apply the theory of micropolar fluids which allows for particles rotation and microinertia. 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 good for colloidal solutions. However, the micropolar fluid theory turned out to be useful in the study of suspension of noncolloidal particles and also including granular fluids. As is proved in [7], it is due to particle rotation that the Segré–Silberberg effect occurs [8]. 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.

Important sectors handling granular fluids include civil engineering (bitumen, concrete, embankments, ballast trains, soil stability), mining (extraction, transport), the chemical industry (fuel and catalysts are often deployed in the form of grains in order to maximize the surface of exchange), the pharmaceutical industry (from the handling of powders for the manufacture of medicine to the handling of drugs themselves) and the food industry (animal food, cereals), to name but a few.

Results from laboratory experiments, numerical simulations, and theoretical approaches from other fields have enriched and renewed our understanding of granular fluids. In many rheological papers on shear flows, particle rotation is ignored. Instead, the authors apply the theory which states that microstructural non-homogeneous particles distribution is due to anisotropy. To this end, the rheological equations involve not only shear stress but normal stress differences as well [9]. The micro-polar fluid theory stands out among other approaches since it handles particle rotations and micro-inertia effects within the mechanics of continua. Such a theory finds many applications in granular flows [10], electrorheology [11,12], ferrofluids [13], visco-elastic micro-polar fluids [14], and liquid crystals [15]. Some results concern blood rheology [16–18]. Particularly, explanations were provided for phenomena including the Fåhræus–Lindqvist effect, the Fåhræus effect, and the plasma skimming effect [19–22]. Any new theory involves additional unknown rheological constants. This restricts applications and further progress. In the micro-polar theory, one such constant is a relative viscosity. What is it? First, we comment on the ordinary fluid viscosity.

In simple shear flows of typical Newtonian or non-Newtonian fluids, viscosity *ηs* [*Pa* · *s*] is given by the formula *η<sup>s</sup>* = *τ*/*γ*˙ or

$$
\pi = \eta\_{\text{s}} \dot{\gamma}\_{\text{\textquotedblleft}} \tag{1}
$$

where *τ* is the shear stress and *γ*˙ is the shear rate. The rates' definition will be given below.

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. Instead of Equation (1), one writes for simple shear flows the following equation:

$$
\pi = \sqrt{(\eta\_{\rm s}\dot{\gamma})^2 + (\eta\_{\rm sk}\dot{\gamma}\_a)^2} \tag{2}
$$

where *γ*˙ *<sup>a</sup>* is a shear rate related to particle rotation, and *ηsk* is an additional viscosity. In what follows, we call *ηsk* skew-symmetrical viscosity motivated by rheological relationships which will be discussed later.

When the volume concentration of particles *φ* is equal to zero, the viscosity *η<sup>s</sup>* of the interstitial fluid is assumed to be known. Clearly, *η<sup>s</sup>* is the usual shear viscosity and it is well known how it can be measured. The skew-symmetric viscosity manifests itself when *φ* > 0, and poor knowledge of it hampers applications of the micro-polar fluid models. The above description of possible applications supports the view that it is of importance to know how skew-symmetric viscosity depends on the particle concentration. The goal of the first sections of the present paper is to get a better insight of how *ηsk* depends on *φ*. In Sections 1–5, we validate some formula for the skew-symmetric viscosity and apply it in Section 6 to the problem of particles migration for the steady flows of yield stress granular fluids in a concentric-cylinder Couette geometry when rotation of the Couette cell inner cylinder is prescribed.

#### **2. Micro-Polar Fluids**

Here, we remind notions concerning a micropolar fluid within the theory of the Cosserat continuum. Such a fluid exhibits micro-rotational effects and micro-rotational inertia; the fluid can support couple stresses and body couples and possesses a non-symmetric stress tensor. The theory of micro-polar fluids goes back to [23–25], where gyration tensor, inertial spin, conservation of micro-inertia, and objectivity of micro-deformation rate tensors are derived and discussed. For an overview of developed theories, we refer the reader to [26,27].

In the Cosserat continuum [28], each material point is treated as a rigid body in the following sense. To such a point with the Lagrange coordinates (*t*, *ξ*), one can assign a position vector **x**(*t*, *ξ*) in the three-dimensional Euclidean space and three orthogonal directors **d***i*(*t*, *ξ*), *i* = 1, 2, 3. Rotation of the vectors **d***<sup>i</sup>* is governed by a rotation orthogonal tensor Q(*t*, *ξ*). The rotation velocity tensor

$$
\Omega(t, \mathbf{x}) = \mathbf{Q}t\mathbf{Q}^\*
$$

is skew-symmetric, and it enjoys the representation formula

$$
\Omega \cdot \mathbf{h} = \omega \times \mathbf{h} \quad \forall \quad \mathbf{h} \in \mathbb{R}^3 \quad (\Omega \cdot \mathbf{h})\_{\mathbf{i}} \equiv \Omega\_{\mathbf{i}\mathbf{j}} h\_{\mathbf{j}\mathbf{k}}
$$

where *ω*(*t*, **x**) is the micro-rotational velocity vector,

$$2\boldsymbol{\omega} = \mathbf{e}\_{i} \times (\boldsymbol{\Omega} \cdot \mathbf{e}\_{i}) = \boldsymbol{\epsilon} : \boldsymbol{\Omega}.$$

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

$$\mathfrak{c}\langle \mathbf{a}, \mathbf{b}, \mathbf{c} \rangle = \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}), \quad \mathbf{e}\_{i} \times \mathbf{e}\_{j} = \mathfrak{e}\_{\mathrm{si}\dagger} \mathbf{e}\_{s\prime} \quad \mathfrak{e}\_{\mathrm{si}\dagger} \equiv \mathfrak{c}\langle \mathbf{e}\_{s\prime} \mathbf{e}\_{i\prime}, \mathbf{e}\_{j} \rangle, \quad (\mathfrak{e} : \Omega)\_{i} \equiv \mathfrak{e}\_{i\dagger} \Omega\_{\mathrm{jk}}.$$

Given a 3 × 3-matrixes *A* and *B*, we use the notation *A* ∗ for the adjoint matrix such that

$$\mathbf{a} \cdot (A \cdot \mathbf{b}) = \mathbf{b} \cdot (A^\* \cdot \mathbf{a}) \quad \forall \mathbf{a}, \mathbf{b} \in \mathbb{R}^3, \quad (A^\*)\_{ij} = A\_{ji}\nu$$

and the scalar product *A* : *B* is defined by *A* : *B* = *AijBij*.

With **v**(*t*, **x**) standing for the velocity of the mass center of the Cosserat material point (*t*, *ξ*), the micropolar fluid is characterized by two rates of strain tensors *B* and *A*:

$$B = \nabla \mathbf{v} - \Omega, \quad A = \nabla \omega. \tag{3}$$

Here, we use the notations (∇**v**)*ij* = *∂vi*/*∂x<sup>j</sup>* , (∇**v**) ∗ *ij* = *∂vj*/*∂x<sup>i</sup>* . An instant stress state of such a fluid is characterized by the couple stress tensor N(*t*, **x**) in addition to the Cauchy stress tensor T(*t*, **x**). Let *S* stand for the viscous part of the stress tensor, T = −*p* I + S. In what follows, we use the symmetric and skew-symmetric parts of a matrix *D*:

$$D\_s = \frac{D + D^\*}{2}, \quad D\_d = \frac{D - D^\*}{2}.$$

In the typical stress–strain relation

$$T = -pI + 2\eta\_s \cdot (\nabla \mathbf{v})\_{s\prime} \tag{4}$$

of Newtonian or non-Newtonian fluids, the scalar factor *η<sup>s</sup>* is the viscosity, *p* is the pressure, with *I* and *T* being the identity and stress tensors. The tensor *T* is symmetric as is well

known in the classical fluid mechanics theory, i.e., *T* ∗ = *T*. Such a symmetry results from momentum and angular momentum laws. On the other hand, the angular momentum law is valid automatically if we postulate symmetry of *T*. This is why nobody invokes such a law in applications.

We remind readers that the constitutive laws of a simple micropolar fluid are [26]

$$T = -pI + S, \quad S = 2\eta\_s B\_s + 2\eta\_{sk} B\_{a\nu} \quad N = 2\gamma A,\tag{5}$$

where *η<sup>s</sup>* , *ηsk* are the viscosities and *γ* is angular viscosity. The first rheological equation in (5) suggests that the contributions of the symmetric part *B<sup>s</sup>* = (∇**v**)*<sup>s</sup>* and skew-symmetric part *B<sup>a</sup>* of the rate of strain tensor *B* into local stress state are different. It is proved in [25] that both the rate of strain tensors *B* and *A* are objective.

Let us introduce the relative angular velocity

$$
\omega\_r = \omega - \text{rot}\,\mathbf{v}/2.
$$

Observe that in fluid mechanics shear stress and shear rates in (1) are defined as follows:

$$
\tau = \sqrt{\mathcal{S} : \mathcal{S}/2}, \quad \dot{\gamma} = \sqrt{2B\_{\text{s}} : B\_{\text{s}}}, \quad \dot{\gamma}\_{\text{d}} = \sqrt{2B\_{\text{d}} : B\_{\text{d}}} = 2|\omega\_{r}|.
$$

Due to the identity 2(∇**v**)*<sup>a</sup>* · **h** = rot **v** × **h** ∀**h**, we have

$$B\_{\mathfrak{a}} \cdot \mathbf{h} = -\omega\_r \times \mathbf{h} \quad \forall \mathbf{h}.$$

Thus, the skew-symmetric viscosity characterizes how relative micro-rotations contribute into the local stress state.

Observe that in the Cosserat continua the Cauchy stress tensor is not symmetric and the vector

$$\mathbf{t} = \mathbf{e}\_i \times (T \cdot \mathbf{e}\_i) = \boldsymbol{\epsilon} : T$$

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

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

The momentum and the angular momentum laws are

$$\frac{\partial(\rho \mathbf{v})}{\partial t} + \text{div}(\rho \mathbf{v} \otimes \mathbf{v}) = -\nabla p + \text{div } \mathbf{S} + \rho \mathbf{f},\tag{7}$$

$$J\left(\frac{\partial(\rho\omega)}{\partial t} + \text{div}(\omega \otimes \rho\mathbf{v})\right) = \text{div}\,N - \varepsilon : \mathbb{S},\tag{8}$$

where *ρ* is the density, *J* is the micro-inertia scalar, **f** is the mass force vector and

$$(\operatorname{div} \mathbf{N})\_{\mathbf{i}} \equiv \partial \mathbf{N}\_{\mathbf{i}\mathbf{j}} / \partial \mathbf{x}\_{\mathbf{j}} \quad (\boldsymbol{\omega} \otimes \mathbf{v})\_{\mathbf{i}\mathbf{j}} = \omega\_{\mathbf{i}} v\_{\mathbf{j}} . \mathbf{x}$$

The density *ρ* satisfies the mass conservation law

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

#### **3. Skew-Symmetrical Viscosity of Dilute Suspensions of Rigid Particles**

Let us consider Couette-like steady flows of suspensions between two parallel planes in the *x*-direction when the upper plane *y* = *h* is fixed and the lower plane *y* = 0 moves in the *x*-direction with the velocity *V*. The volume particle concentration *φ* is assumed to be fixed.

We outline the method for determination of the skew-symmetrical viscosity *ηsk*. Assume that the stress applied to the moving plate can be measured. By continuity, one can

tell the fluid stress in the nearby fluid region. On the other hand, one can calculate such a fluid stress by one or another mathematical model. First, we determine the stress *S m* 1 at the moving plane by applying the micro-polar fluid theory. Then, we calculate the stress *S v* 1 at the moving plane by applying the Navier–Stokes theory. We equate the two stresses and derive the skew-symmetrical viscosity *ηsk* from the equality *S m* <sup>1</sup> = *S v* 1 .

Let us first treat the suspension as a micropolar fluid with the prescribed volume particle concentration *φ*. The above assumptions upon the flows suggest that the velocity vector **v**, the micro-rotational velocity vector *ω* and the pressure *p* depend on the vertical variable *y* only:

> **v** = *v*(*y*)(1, 0, 0) *T* , *ω* = *ω*(*y*)(0, 0, 1) *T* , *p* = *p*(*y*), 0 < *y* < *h*.

For such flows, the matrices ∇**v** and Ω become

$$
\nabla \mathbf{v} = \begin{pmatrix} 0 & v\_y & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \Omega = \begin{pmatrix} 0 & -\omega & 0 \\ \omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \frac{\partial v}{\partial y} = v\_y.
$$

Hence,

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

Projections of the momentum and the angular momentum equations on the *x* and *z*-directions become

$$0 = \frac{\partial S\_{12}}{\partial y},$$

$$0 = \frac{\partial N\_{32}}{\partial y} + \mathbf{S}\_{21} - \mathbf{S}\_{12\prime} \tag{11}$$

respectively, where the tensor components are given by the formulas

$$S\_{21} = \mathbf{e}\_y \cdot S \langle \mathbf{e}\_x \rangle, \quad S\_{12} = \mathbf{e}\_x \cdot S \langle \mathbf{e}\_y \rangle, \quad N\_{32} = \mathbf{e}\_z \cdot N \langle \mathbf{e}\_y \rangle.$$

Constitutive laws (5) take the form

$$\mathcal{S}\_{12} = (\eta\_{\rm s} + \eta\_{\rm sk}) \frac{\partial v}{\partial y} + 2\eta\_{\rm sk}\omega, \quad \mathcal{S}\_{21} = (\eta\_{\rm s} - \eta\_{\rm sk}) \frac{\partial v}{\partial y} - 2\eta\_{\rm sk}\omega. \tag{12}$$

We note that system (10)–(11) does not contain pressure. As is well known, it can be restored from the momentum equation projected on the *y*-direction.

The viscosity *ηsk*(*φ*) vanishes when *φ* → 0. The same is true for the relative viscosity

$$
\varepsilon(\phi) = \frac{\eta\_{sk}(\phi)}{\eta\_s},
$$

which we represent via the expansion series

$$
\varepsilon(\phi) = \Lambda \phi + \Lambda\_2 \phi^2 + \dotsb \cdot \mathbb{1}
$$

Whereas the velocity **v** satisfies the no-slip boundary conditions, we require that

$$
\omega = \mathfrak{a}(\phi) \mathbf{rot} \,\mathbf{v}/2 \quad \text{at} \quad y = 0 \quad \text{and} \quad y = h,\tag{13}
$$

where *α*(*φ*) = *α*0*φ*. The latter condition implies that the micro-rotations agree with macrorotations at the boundary [29]. For the Couette-like flows, we arrive at the following boundary-value problem in the domain 0 < *y* < *h*:

$$\frac{\partial}{\partial y}\left[ (1+\varepsilon(\phi)) \frac{\partial v}{\partial y} + 2\varepsilon(\phi)\omega \right] = 0,\tag{14}$$

$$2\gamma \frac{\partial^2 \omega}{\partial y^2} - 2\eta\_s \varepsilon(\phi) \left(\frac{\partial v}{\partial y} + 2\omega\right) = 0,\tag{15}$$

$$
v|\_{y=0} = V, \quad v|\_{y=h} = 0, \quad \omega|\_{y=0,h} = -\frac{a\_0 \phi}{2} \frac{\partial v}{\partial y}\Big|\_{y=0,h}.\tag{16}$$

We solve the boundary value problem (14)–(16) looking for (*v*, *w*) as the asymptotic expansion series

$$v(y,\phi) = v^0(y) + v^1(y)\phi + \cdots, \quad \omega(y,\phi) = \omega^0(y) + \omega^1(y)\phi + \cdots \tag{17}$$

Setting these series in (14)–(16), one can write each equality in (14)–(16) in the form

$$
\phi^0(\cdot \cdot \cdot)\_0 + \phi^1(\cdot \cdot \cdot)\_1 + \dots = 0.
$$

The coefficients *v <sup>i</sup>* and *<sup>ω</sup><sup>j</sup>* are determined from the conditions (· · ·)*<sup>k</sup>* <sup>=</sup> 0 for any *<sup>k</sup>*.

Particularly, if *k* = 0, we derive the following boundary value problems for the functions *v* 0 (*y*), *ω*<sup>0</sup> (*y*):

$$\frac{\partial^2 v^0}{\partial y^2} = 0, \quad v^0|\_{y=0} = V, \quad v^0|\_{y=h} = 0,\tag{18}$$

$$0 = \frac{\partial^2 \omega^0}{\partial y^2}, \quad \omega^0|\_{y=0} = \omega^0|\_{y=h} = 0. \tag{19}$$

Similarly, we find that the function *v* 1 satisfies the boundary value problem

$$
\frac{
\partial
}{
\partial y
}
\left[
\frac{
\partial v^1
}{
\partial y
} + \Lambda \left(
\frac{
\partial v^0
}{
\partial y
} + 2\omega^0
\right)
\right] = 0, \quad v^1|\_{y=0} = v^1|\_{y=h} = 0. \tag{20}
$$

Solving these problems, we find that

$$
\omega^0 = 0, \quad v^1 = 0, \quad v^0 = V(1 - y/h).
$$

Starting from the definitions (12) related to the micro-polar fluid theory, we can write the expansion series for the stress *S*12(*φ*) as follows:

$$\mathcal{S}\_{12}(\phi) = \mathcal{S}\_{12}^0 + \mathcal{S}\_{12}^1 \phi + o(\phi).$$

Clearly,

$$S\_{12}^{0} = \eta\_s \frac{\partial v^0}{\partial y}, \quad S\_{12}^{1} = \eta\_s \frac{\partial v^1}{\partial y} + \eta\_s \Lambda \left(\frac{\partial v^0}{\partial y} + 2\omega^0\right).$$

Now, we can calculate the relative stress at the moving plane:

$$\frac{S\_{12}^m(\phi)}{S\_{12}^m(0)} \equiv \frac{S\_{12}(\phi)}{S\_{12}(0)}\Big|\_{y=0} = 1 + \phi \Lambda + o(\phi). \tag{21}$$

It is assumed that both the stresses *S m* <sup>12</sup>(*φ*) and *S <sup>m</sup>*(0)<sup>12</sup> are measured at the same velocity *V* of the moving plate.

Let us consider flows within the same Couette geometry, starting from the Navier– Stokes theory. In such a theory, the stress tensor *S* is symmetric. Denoting *S* = *S*12, one can find the velocity *v*(*y*) by solving the boundary-value problem

$$0 = \frac{\partial S}{\partial y'}, \quad S = \eta(\phi) \frac{\partial v}{\partial y'}, \quad v|\_{y=0} = V, \quad v|\_{y=h} = 0,\tag{22}$$

where *η*(*φ*) is the effective viscosity.

Let *S v* 1 stand for the stress *S* at the moving plate, *S v* <sup>1</sup> = *S*|*y*=0. Given *S v* 1 , one can derive from (22) the following formula for the apparent viscosity *η*(*φ*):

$$\eta(\phi) = -\frac{hS\_1^v}{V} \tag{23}$$

Observe that *S v* 1 has negative values since *v*(*y*) is a decreasing function of *y*. It follows from (23) that

$$\frac{\eta(\phi)}{\eta(0)} = \frac{S\_1^v(\phi)}{S\_1^v(0)}, \quad \eta(0) \equiv \eta\_s. \tag{24}$$

It is assumed that both the stresses *S v* 1 (*φ*) and *S*1(0) *<sup>v</sup>* are measured at the same velocity *V* of the moving plate.

For dilute suspensions, the left-hand side of (24) is given by the Einstein formula [30]

$$\frac{\eta(\phi)}{\eta\_s} = 1 + E\phi + o(\phi), \quad E \simeq 2.5. \tag{25}$$

We equate the relative stresses: *S m* <sup>12</sup>(*φ*)/*S m* <sup>12</sup>(0) = *S v* 1 (*φ*)/*S v* 1 (0). Now, it follows from (21), (24) and (25) that

$$1 + E\phi = 1 + \phi\Lambda + o(\phi).$$

Hence, Λ = *E*. By the above arguments, we conclude that the skew-symmetric viscosity for dilute suspensions satisfies the representation formula

$$
\eta\_{\rm sk}(\phi) / \eta\_{\rm s} = E\phi + o(\phi),
\tag{26}
$$

where *E* ' 2.5 is the Einstein factor.

**Remark 1.** *By the same asymptotic arguments, we can conclude that for the Couette flows between two concentric cylinders formula* (26) *becomes*

$$
\eta\_{\rm sk}(\phi) / \eta\_{\rm s} = \mathbf{G} \cdot E\phi + o(\phi) \,, \tag{27}
$$

*where G is the geometrical factor equal to* (*R*2/*R*1) 2 *, with R*<sup>1</sup> *being the smaller radius.*

One more conclusion from the above arguments is that there is a correlation between the Navier–Stokes theory and the micro-polar fluid theory:

$$
\eta(\phi)/\eta(0) = 1 + \eta\_{sk}(\phi)/\eta\_s + o(\phi), \quad \eta\_s = \eta(0), \tag{28}
$$

where *η<sup>s</sup>* , *ηsk* are the micro-polar fluid viscosities of the suspension and *η*(*φ*) is the apparent viscosity of the same suspension described by the Navier–Stokes rheology. The law (28) is verified by the asymptotic series argument for dilute suspensions, with the left-hand side given by the Einstein law *η*(*φ*) = *η*(0)(1 + *Eφ*)+ *o*(*φ*).

On the other hand, there is an extended Krieger–Douhgerty empirical closure [31]

$$
\eta(\phi)/\eta(0) = (1 - \phi/\phi^\*)^{-E\phi^\*}, \quad E = 2.5,\tag{29}
$$

for dense suspensions, where *φ* ∗ is a maximal volume concentration. Such a closure suggests that, by setting (29) in (28), we can define the skew-symmetric viscosity *ηsk*(*φ*) as follows:

$$(1 - \phi/\phi^\*)^{-E\phi^\*} = 1 + \eta\_{sk}(\phi)/\eta\_s. \tag{30}$$

In the next sections, we verify this empirical formula by studying flows of dense suspensions paying attention to particle rotation.

#### **4. Flows of Suspensions of Rigid Particles in the Herschel–Bulkley Fluid**

By the theoretical approach of Chateau et al. [32], Ovarlez et al. [9] showed that the flows of yield stress suspensions in a concentric-cylinder Couette geometry can be modeled by a Herschel–Bulkley behavior of same index as their interstitial fluid. The theory was proved to be in an agreement with the laboratory experiments based on the magnetic resonance imaging techniques [9].

We are going to address the same experiments as in [9] to verify formula (30). First, we extend the constitutive laws (5) to allow for the yield stress rheology. According to [33], the Cosserat–Bingham fluid rheology is defined as follows:

$$\mathbf{S}^{\prime} = \begin{cases} \ 2\eta\_{s}B\_{\mathrm{s}} + 2\eta\_{s\mathrm{k}}B\_{\mathrm{a}} + \tau\_{\ast}\frac{B\_{\mathrm{0}}}{|B\_{\mathrm{0}}|^{\prime}} & \text{if} \quad \mathsf{B}(\mathbf{x},t) \neq \mathbf{0},\\\ \mathsf{S}\_{\mathrm{p}}(\mathbf{x},t), & \text{if} \quad \mathsf{B}(\mathbf{x},t) = \mathbf{0}, \end{cases} \tag{31}$$

$$\mathbf{N}\_{\perp} = \begin{cases} \ 2\gamma A + \pi\_{\mathbb{R}} \frac{A}{|A|}, & \text{if} \quad A(\mathbf{x}, t) \neq \mathbf{0}, \\\ \mathbf{N}\_{p}(\mathbf{x}, t), & \text{if} \quad A(\mathbf{x}, t) = \mathbf{0}, \end{cases} \tag{32}$$

where

$$B\_0 = B\_s + \varepsilon B\_{a\nu} \quad \varepsilon(\phi) = \frac{\eta\_{sk}(\phi)}{\eta\_s}$$

and *τ*<sup>∗</sup> and *τ<sup>n</sup>* are yield stresses; the unknown plug tensors S*<sup>p</sup>* and N*<sup>p</sup>* obey the restrictions

$$|\mathbf{S}\_{p}| \le \mathfrak{r}\_{\star \prime} \quad |\mathbf{N}\_{p}| \le \mathfrak{r}\_{\mathbf{n}}.$$

It is proved in [34] that the formulation (31) is equivalent to the inclusion *S* ∈ *∂V*∗(*B*0), where the scalar potential *V*∗(*D*) is defined for any matrix *D* ∈ *R* <sup>3</sup>×<sup>3</sup> by the formula *V*∗(*D*) = *η<sup>s</sup>* |*D*| <sup>2</sup> <sup>+</sup> *<sup>τ</sup>*∗|*D*|. We remind that the subdifferential formulation *<sup>S</sup>* <sup>∈</sup> *<sup>∂</sup>V*∗(*B*0) implies that

$$S: (D - B\_0) \le V\_\*(D) - V\_\*(B\_0) \quad \text{for all} \quad D \in \mathbb{R}^{3 \times 3}.$$

Similarly, the constitutive law (32) is equivalent to the inclusion *N* ∈ *∂Vn*(*A*) with *Vn*(*D*) = *γ*|*D*| <sup>2</sup> <sup>+</sup> *<sup>τ</sup>n*|*D*|, <sup>∀</sup>*<sup>D</sup>* <sup>∈</sup> *<sup>R</sup>* 3×3 . The meaning of the plug zone |*N*(*x*, *t*)| ≤ *τ<sup>n</sup>* is discussed in [35].

Let *T*<sup>0</sup> be a characteristic time. We denote the dimensionless second invariant of the rate of strain tensor *B*<sup>0</sup> by *I*: *I* = *T*0|*B<sup>s</sup>* |. To transform the Cosserat–Bingham fluid constitutive laws (31) and (32) into the Cosserat–Herschel–Bulkley fluid rheological equations, we assume that

$$\eta\_{\rm s} = \eta\_{\rm s0} \mathbf{I}^{n-1} \quad \text{and} \quad \frac{\eta\_{\rm sk}(\boldsymbol{\phi})}{\eta\_{\rm s}} = \boldsymbol{\varepsilon}(\boldsymbol{\phi}), \quad \text{where} \quad \boldsymbol{\varepsilon}(\boldsymbol{\phi}) = (1 - \boldsymbol{\phi}/\boldsymbol{\phi}^{\*})^{-\mathbb{E}\boldsymbol{\phi}^{\*}} - 1. \tag{33}$$

Observe that the classical Hershel–Bulkley model results from the constitutive laws (31)–(33) if the particle volume concentration *φ* vanishes.

We consider steady axially symmetric flows of a suspension between two coaxial cylinders centered on the *z*-axis. The inner cylinder of the radius *R*<sup>1</sup> rotates with the angular velocity Ω0[s −1 ] and the external cylinder of the radius *R*<sup>2</sup> is fixed. The volume particle concentration *φ* is assumed invariable along the radial coordinate. The case of variable *φ* will be addressed in Section 6.

In what follows, we use the unit vectors **e***<sup>r</sup>* , **e***ϕ*, **e***<sup>z</sup>* of the cylindrical coordinate system. The assumption of axially symmetry of flows suggests that the velocity vector **v**, the microrotational velocity vector *ω* and the pressure *p* depend on the radial variable *r* only:

$$\mathbf{v} = v(r)\mathbf{e}\_{\varphi} \quad \omega = \omega(r)\mathbf{e}\_{z\prime} \quad p = p(r), \quad \text{rot}\,\mathbf{v} = \left(\frac{\partial v}{\partial r} + \frac{v}{r}\right)\mathbf{e}\_{z\prime} \quad \text{rot}\,\omega = -\frac{\partial \omega}{\partial r}\mathbf{e}\_{\Phi}. \tag{34}$$

For such flows, the matrices ∇**v** and Ω in the cylindrical coordinate system become

$$
\nabla \mathbf{v} = \begin{pmatrix} 0 & -v/r & 0 \\ \partial v & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \Omega = \begin{pmatrix} 0 & -\omega & 0 \\ \omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad \nabla \omega = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ \partial \omega & 0 & 0 \end{pmatrix},
$$

where we denoted *∂v*/*∂r* by *∂v* for simplicity. Here,

$$(\nabla \mathbf{v})\_{\vec{i}\vec{j}} = \mathbf{e}\_i \cdot \mathbf{v}'(\mathbf{x}) \langle \mathbf{e}\_{\vec{j}} \rangle, \quad \mathbf{i} = (r, \varphi, z), \quad \mathbf{v}'(\mathbf{x}) \langle \mathbf{a} \rangle = \frac{d}{d\lambda} \mathbf{v}(\mathbf{x} + \lambda \mathbf{a}) \big|\_{\lambda=0}.$$

Hence,

$$B = \begin{pmatrix} 0 & -v/r + \omega & 0 \\ \partial v - \omega & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix},$$

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

As for the rate of strain tensor *A*, we have that *Azr* = *∂ω*/*∂r* and *Aij* = 0, otherwise. Projections of the momentum Equation (7) and the angular momentum law (8) on the vectors **e***<sup>ϕ</sup>* and **e***<sup>z</sup>* become

$$0 = \frac{\partial S\_{\varphi r}}{\partial r} + \frac{S\_{\varphi r} + S\_{r\varphi}}{r},\tag{35}$$

$$0 = \frac{\partial N\_{zr}}{\partial r} + \frac{N\_{zr}}{r} + S\_{qr} - S\_{r\varphi\nu} \tag{36}$$

respectively, where the tensor components are given by the formulas

$$S\_{\varphi r} = \mathbf{e}\_{\varphi} \cdot \mathcal{S} \langle \mathbf{e}\_{r} \rangle, \quad S\_{r\varphi} = \mathbf{e}\_{r} \cdot \mathcal{S} \langle \mathbf{e}\_{\varphi} \rangle, \quad N\_{zr} = \mathbf{e}\_{z} \cdot N \langle \mathbf{e}\_{r} \rangle.$$

Given a vector **e**, we apply the notation (*S*h**e**i)*<sup>i</sup>* = *Sije<sup>j</sup>* . Observe that the other components of *S* and *N* are equal to zero. In what follows, we use the equation

$$\frac{\rho v^2}{r} = \frac{\partial p}{\partial r} \, \tag{37}$$

resulting from projection of the momentum equation (7) onto the vector **e***<sup>r</sup>* .

We calculate that

$$I = 2^{-1/2} T\_0 \sqrt{(\partial v/\partial r - v/r)^2 + \varepsilon^2 (\partial v/\partial r + v/r - 2\omega)^2}.$$

The constitutive laws (31)–(32) become

$$S\_{r\varphi} = \left(2\eta\_s + \frac{\tau\_\* T\_0}{I}\right) \left[\frac{1}{2}\left(\frac{\partial v}{\partial r} - \frac{v}{r}\right) - \frac{\varepsilon}{2}\left(\frac{\partial v}{\partial r} + \frac{v}{r} - 2\omega\right)\right] \quad \text{if} \quad I \neq 0,\tag{38}$$

$$S\_{\varphi r} = \left(2\eta\_s + \frac{\tau\_\* T\_0}{I}\right) \left[\frac{1}{2}\left(\frac{\partial v}{\partial r} - \frac{v}{r}\right) + \frac{\varepsilon}{2}\left(\frac{\partial v}{\partial r} + \frac{v}{r} - 2\omega\right)\right] \quad \text{if} \quad I \neq 0,\tag{39}$$

$$|S\_{r\varphi}|^2 + |S\_{\varphi r}|^2 \le \tau\_\*^2 \quad \text{if} \quad I = 0,\tag{40}$$

$$N\_{\rm zr} = 2\gamma \frac{\partial \omega}{\partial r} + \tau\_{\rm n} \text{sign}\left(\frac{\partial \omega}{\partial r}\right) \quad \text{if} \quad \frac{\partial \omega}{\partial r} \neq 0, \quad |N\_{\rm zr}| \le \tau\_{\rm n} \quad \text{if} \quad \frac{\partial \omega}{\partial r} = 0,\tag{41}$$

The boundary condition (13) for the angular velocity *ω* and the no-slip condition for *v* become

$$
\omega|\_{r=R\_i} = \frac{a\_0 \phi}{2} \left(\frac{\partial v}{\partial r} + \frac{v}{r}\right)\Big|\_{R\_i}, \quad v|\_{r=R\_1} = R\_1 \Omega, \quad v|\_{r=R\_2} = 0. \tag{42}
$$

To study numerically the boundary-value problem (33)–(42) in the annulus *R*<sup>1</sup> < *r* < *R*2, we pass to dimensionless variables:

$$\sigma' = \frac{r}{R\_1}, \quad v' = \frac{v}{V}, \quad \omega' = \frac{\omega}{\omega\_0}, \quad S\_{r\rho}' = \frac{S\_{r\rho}}{S^0}, \quad S\_{\varphi r}' = \frac{S\_{\varphi r}}{S^0}, \quad N\_{zr}' = \frac{N\_{zr}}{N\_0}, \quad \gamma\_1 = \frac{\gamma}{R\_1^2 \eta\_{s0}}, \quad \omega' = \frac{\gamma}{S^0}, \quad \omega' = \frac{\gamma}{S^0}, \quad \omega' = \frac{\gamma}{S^0}, \quad \omega' = \frac{\gamma}{S^0}, \quad \omega' = \frac{\gamma}{S^0}, \quad \omega' = \frac{\gamma}{S^0}, \quad \omega' = \frac{\gamma}{S^0}$$

with

$$V = R\_1 \Omega\_0, \omega\_0 = \Omega\_0, \ T\_0 = \frac{1}{\Omega\_0} [\mathbf{s}], \quad \mathbb{S}\_0 = \eta\_{s0} \Omega\_0, \ \tau\_{s1} = \frac{\tau\_\*}{\eta\_{s0} \Omega\_0}, \ \tau\_{\eta 1} = \frac{\tau\_{\eta} R\_1}{\gamma \Omega\_0}, \ N\_0 = R\_1 \eta\_{s0} \Omega\_0.$$

Observe that the dimensionless yield stress *τ*∗<sup>1</sup> is the inverse of the Bingham number for the Couette flows:

$$
\pi\_{\*1} = \frac{1}{Bn'} \quad Bn = \frac{\eta\_{s0} \Omega\_0}{\tau\_\*}.
$$

In new variables,

$$I = 2^{-1/2} \sqrt{(\partial' v'/\partial r' - v'/r')^2 + \epsilon^2 (\partial' v'/\partial r' + v'/r' - 2\omega')^2}, \quad \frac{\eta\_s}{\eta\_{s0}} = I^{n-1}.\tag{43}$$

$$S'\_{r\rho} = \left[\frac{1}{2}\left(\frac{\partial' v'}{\partial r'} - \frac{v'}{r'}\right) - \frac{\varepsilon}{2}\left(\frac{\partial' v}{\partial r'} + \frac{v'}{r'} - 2\omega'\right)\right] \left(2I^{n-1} + \frac{\tau\_{\ast 1}}{I}\right) \quad \text{if} \quad I \neq 0,\tag{44}$$

$$S\_{qr}' = \left[\frac{1}{2}\left(\frac{\partial'v'}{\partial r'} - \frac{v'}{r'}\right) + \frac{\varepsilon}{2}\left(\frac{\partial'v'}{\partial r'} + \frac{v'}{r'} - 2\omega'\right)\right]\left(2I^{n-1} + \frac{\tau\_{\ast 1}}{I}\right) \quad \text{if} \quad I \neq 0,\tag{45}$$

$$|S\_{r\rho}'|^2 + |S\_{\varphi r}'|^2 \le \pi\_{\*1}^2 \quad \text{if} \quad I = 0,\tag{46}$$

$$N\_{zr}' = \gamma\_1 \left( 2 \frac{\partial' \omega'}{\partial r'} + \tau\_{\text{n1}} \text{sign} \frac{\partial' \omega'}{\partial r'} \right) \quad \text{if} \quad \frac{\partial' \omega'}{\partial r'} \neq 0, \quad |N\_{zr}'| \le \tau\_{\text{n1}} \quad \text{if} \quad \frac{\partial' \omega'}{\partial r'} = 0,\tag{47}$$

$$0 = \frac{\partial' S'\_{qr'}}{\partial r'} + \frac{S'\_{qr} + S'\_{rq}}{r'},\tag{48}$$

$$0 = \frac{\partial' \mathcal{N}\_{\mathbb{Z}r}'}{\partial r'} + \frac{\mathcal{N}\_{\mathbb{Z}r}'}{r'} + \mathcal{S}\_{qr}' - \mathcal{S}\_{r\varphi}' \tag{49}$$

$$
\omega'|\_{r'=1,a} = \frac{a\_0 \phi}{2} \left( \frac{\partial' v'}{\partial r'} + \frac{v'}{r'} \right) \Big|\_{r'=1,a'} \quad v'|\_{r'=1} = 1, \quad v'|\_{r'=a} = 0. \tag{50}
$$

#### **5. Skew-Symmetric Viscosity versus Particles Concentration**

Here, we apply the mathematical model developed in the previous section to justify formula (30) for the skew-symmetric viscosity. To perform calculations, we fix parameters of the interstitial Hershel–Bilkley fluid. Rheological constitutive law of such a fluid results from (31) by setting *φ* = 0:

$$\mathbf{S}\_{\perp} = \begin{cases} 2\eta\_{\text{s}} \mathbf{B}\_{\text{s}} + \tau\_{\ast} \frac{\mathbf{B}\_{\text{s}}}{|\mathbf{B}\_{\text{s}}|}, & \text{if} \quad \mathbf{B}\_{\text{s}}(\mathbf{x}, t) \neq 0, \\\ \mathbf{S}\_{\text{p}}(\mathbf{x}, t), & \text{if} \quad \mathbf{B}\_{\text{s}}(\mathbf{x}, t) = \mathbf{0}, \end{cases} \quad \text{where} \quad \eta\_{\text{s}} = \eta\_{\text{s}0} (T\_0 |\mathbf{B}\_{\text{s}}|)^{\eta - 1}. \tag{51}$$

We denote

$$
\pi\_y = \sqrt{2}\pi\_{\ast \prime} \quad \eta\_{HB} = 2^{5/4} \eta\_{s0} T\_0^{-1/2} \tag{52}
$$

and set *n* = 1/2, *R*<sup>1</sup> = 4[cm], *R*<sup>2</sup> = 6[cm]. Then, it follows from (51) that

$$
\sqrt{2S:S} = \tau\_{\mathcal{Y}} + \eta\_{HB} (\sqrt{2B\_{\mathcal{S}}:B\_{\mathcal{S}}})^{1/2}.\tag{53}
$$

The concentrated emulsion obeying Equation (53) was considered in [9] with *τ<sup>y</sup>* = 22 [Pa] and *ηHB* =5.3 [Pa·s 1/2]. Given the angular velocity Ω [rpm] of the rotating inner cylinder, we calculate that *T*<sup>0</sup> = 1/Ω<sup>0</sup> = (60/Ω)[s]. We consider the same fluid as in [9], hence one can define the consistency *ηs*<sup>0</sup> and the yield stress *τ*<sup>∗</sup> as follows:

$$\eta\_{\rm s0} = 2^{-5/4} \eta\_{HB} [\text{Pa} \cdot \text{s}^{1/2}] \left(\frac{60}{\Omega}\right)^{1/2} [\text{s}^{1/2}], \quad \tau\_\* = \tau\_y / \sqrt{2} [\text{Pa}]. \tag{54}$$

Observe that consistency depends on Ω because of the special choice the characteristic time *T*<sup>0</sup> and the definition of the dimensionless invariant *I* of the rate of strain tensor *B<sup>s</sup>* .

It is proved in [35] that the rotation yield stress *τ<sup>n</sup>* causes the appearance of clusters of particles, with each cluster being a plug zone which rotates as a rigid body. Conglomerates of particles were not observed in [9] for the Couette flows of suspensions between two rotating cylinders; this is why *τ<sup>n</sup>* and *τn*<sup>1</sup> can be neglected. As for the dimensionless angular viscosity *γ*<sup>1</sup> and the boundary-value dimensionless parameter *α*0, we variate them to fit experimental data.

Approximate solutions of the system (43)–(50) can be obtained by regularization [36]. Given a small positive *δ*, we substitute the dimensionless invariant *I* in (44) and (45) by *I<sup>δ</sup>* , where

$$I\_\delta = 2^{-1/2}\sqrt{(\partial'v'/\partial r' - v'/r')^2 + \varepsilon^2(\partial'v'/\partial r' + v'/r' - 2\omega')^2 + \delta^2}\,\_1$$

with *δ* & 0.

First, we tune the model (43)–(50) by setting *φ* = 0 and addressing the pure interstitial Herschel–Bulkley fluid with *τy*(0) = 22[Pa] and *ηHB*(0) = 5.3[Pa · s 1/2] as in [9]. Equations become

$$I = 2^{-1/2} \sqrt{(\partial^\prime v^\prime / \partial r^\prime - v^\prime / r^\prime)^2}, \quad \frac{\eta\_s}{\eta\_{s0}} = I^{n-1}, \quad n = 1/2. \tag{55}$$

$$S'\_{r\!\!\!\!} = S'\_{\!\!\!\!\!\!\!\!\!} = \frac{1}{2} \left( \frac{\partial' v'}{\partial r'} - \frac{v'}{r'} \right) \left( 2I^{n-1} + \frac{\pi\_{\ast 1}}{I} \right) \quad \text{if} \quad I \neq 0,\tag{56}$$

$$2|\mathcal{S}'\_{r\rho}|^2 \le \tau\_{\*1}^2 \quad \text{if} \quad I = 0,\tag{57}$$

$$0 = \frac{\partial' S'\_{\varphi r'}}{\partial r'} + \frac{2S'\_{\varphi r}}{r'},\tag{58}$$

$$|v'|\_{r'=1} = 1, \quad v'|\_{r'=\mathcal{R}\_2/\mathcal{R}\_1} = 0. \tag{59}$$

Observe that in such a case the model (55)–(59) depends on one parameter *τ*∗1(0) only. Given the angular velocity Ω[s −1 ], we find from (52) the value of the Bingham number *Bn*(0) by the formula

$$Bn(0) = \frac{\eta\_{HB}(0)\Omega\_0^{1/2}}{2^{3/4}\tau\_{\mathcal{Y}}(0)}.\tag{60}$$

Figures 1 and 2 depict very good agreement of calculation results for *φ* = 0 with experiment data [9] for different angular velocities Ω0[s −1 ] but in the case of passage to the effective Bingham number *Bne*(0):

$$B n\_{\varepsilon}(0) = 1.5 \cdot B n(0),\tag{61}$$

We think that such a discrepancy between the measured and effective Bingham numbers is due to the following reasons. Real 3*D*-flows are described by 1*D*-equations. The gravitation, the height of the annulus region, and the lateral boundaries effect are not taken into account. It may be that viscoelastic fluid property is also of importance, which calls for more adequate modeling.

Figure 1 corresponds to Ω = 2, Ω = 5 and Ω = 100[rpm]. The same agreement between calculation results and experiment data is observed for Ω = 10, 20 and 50[rpm], but we omit pictures to save the space. Figure 2a combines all velocity profiles for Ω = 2, 5, 10, 20, 50 and 100[rpm]; it fits the laboratory experiments exposed in Figure 2b. Why does the velocity profile become less steep as Ω increases? Equations (60) and (61) answer the question. In fact, there is a motionless plug zone of the Herschel–Bulkley fluid near the exterior cylinder. Our approximate solutions based on the regularization approach do not catch the plug zone well. The bigger the plug zone, the steeper the velocity curve. However, the dimensionless yield stress *τ*∗1(0) stipulates the size of the plug zone and, at the same time, it varies inversely with the angular velocity Ω.

Now, we consider suspensions assuming that they can be modeled by a Herschel– Bulkley behavior of the same index as their interstitial fluid, with their consistency and their yield stress depending on the particle volume fraction. To determine the function *τ*∗1(*φ*) in Equations (55)–(59), we apply correlations proposed in [9] for *n* = 1/2:

$$\frac{\eta\_{HB}(\phi)}{\eta\_{HB}(0)} = \left(\frac{\tau\_y(\phi)}{\tau\_y(0)}\right)^{3/2} \cdot (1-\phi)^{-1/2} \text{.}\tag{62}$$

$$\frac{\tau\_y(\phi)}{\tau\_y(0)} = \sqrt{(1-\phi)(1-\phi/\phi^\*)^{-2.5\phi^\*}}\tag{63}$$

In our notations (see (33)), equality (63) becomes

$$\frac{\tau\_y(\phi)}{\tau\_y(0)} = \sqrt{(1-\phi)(\varepsilon(\phi)-1)}.$$

It follows from the definitions (60) and (61) that

$$\frac{\pi\_{\*1}(\phi)}{\pi\_{\*1}(0)} = \left(\frac{1-\phi}{1+\varepsilon(\phi)}\right)^{1/4}.\tag{64}$$

With the function *τ*∗1(*φ*) given by (64), we solve Equations (55)–(59) and find an agreement with experiment data [9]. Calculations and laboratory data depicted in Figure 3 are related to *φ* = 0.3 for Ω = 2, 5, 10, 20, 50 and 100[rpm].

Let us return to the general system (43)–(50) which describes micropolar fluid. We fix *τ*∗<sup>1</sup> by Equations (64), (60) and (61). To choose the dimensionless parameters *γ*<sup>1</sup> and *α*0, we apply the method of least squares based on minimizing the function

$$F(\gamma\_1, \mathfrak{a}\_0) = \sum\_{i,j,k} |v'\_{\gamma\_1, \mathfrak{a}\_0}(r\_{i\nu} \phi\_{j\nu} \Omega\_k) - v'\_{data}(r\_{i\nu} \phi\_{j\nu} \Omega\_k)|^2.$$

Here, *v* 0 *data* is the measured velocity at different locations *r<sup>i</sup>* for different volume concentrations *φ<sup>j</sup>* and different angular velocities Ω*<sup>k</sup>* , *v* 0 *γ*1 ,*α*<sup>0</sup> (*r*, *φ*, Ω) is the calculated velocity, with *γ*<sup>1</sup> and *α*<sup>0</sup> being prescribed. Calculations reveal that the optimal *γ*<sup>1</sup> and *α*<sup>0</sup> take values *γ* ∗ <sup>1</sup> = 10.93 and *α* ∗ <sup>0</sup> = 0.79 provided *τy*(0) = 22[*Pa*] and *ηHB*(0) = 5.3[*Pa* · s 1/2]. Below, we provide results of calculations with the chosen data *τ*∗1, *γ* ∗ 1 and *α* ∗ 0 .

Figures 4 and 5a concern calculations for Ω = 100[rpm] when *φ* takes values 0, 0.1 and 0.3. Measured data in Figure 5b borrowed from [9] confirm agreement with calculations. The same is true for Ω = 5[rpm] when *φ* takes values 0, 0.1 and 0.3 as shown in Figures 6 and 7. In Figures 8 and 9, *φ* is fixed equal to 0.3 with Ω taking on the values 2, 5, 10, 20, 50, and 100[rpm]. Calculations agree with the measured data from [9].

The micro-polar fluid rheology equations (5) predict particle rotation. Figure 10 depicts profiles of the dimension angular velocity *w*(*r*) when *φ* is fixed equal to 0.3 with Ω taking on the values 2, 5, 10, 20, 50, and 100[rpm].

Although the dimensionless angular viscosity *γ*<sup>1</sup> = 2 5/4*γR* −2 1 *η* −1 *HB*(*φ*)Ω is determined, we can not identify the dimensional angular viscosity *γ*. Indeed, the tuning step (61) implies that we substituted *τy*(0)/*ηHB*(0) by 0.2 · *τy*(0)/*ηHB*(0). However, in doing so, it is impossible to know individual reduced values both of *τy*(0) and *ηHB*(0). Hence, we don't know the reduced value of *ηHB*(*φ*).

Let us comment on some discrepancy between calculations and data of measurement. One can see in Figures 4b and 7 that, with increasing the rotational velocity at a constant volume concentration *φ* = 0.1, the theoretical results are shown to better agree with experiment. It is a problem of calculations. The reason is that we consider the viscoplastic fluid and the plug zone (with zero velocity and low shear stress) near the external cylinder being

larger at low rotations. The governing system of equations becomes degenerate in such a zone. Mathematical theory of degenerate systems of differential equations and corresponding numerical methods built into Wolfram Mathematica are not well developed yet. There is one more difficulty related to the Hershel–Bulkley fluid viscosity (with the power *n* = 1/2) becoming infinite when the velocity gradient vanishes somewhere. To get over these difficulties, we apply the regularization technique and substitute the invariant *I* of the rate of strain tensor *B*<sup>0</sup> in Equations (38)–(40) by its non-vanishing approximation *I<sup>δ</sup>* .

A comparison of the results in Figure 1 (upper curve), and Figure 4 at the rotational velocity Ω = 100 rpm shows that, with increasing the volume concentration, the calculated results better agree with experiments at intermediate concentrations *φ*=0.1; at lower (*φ* = 0) and at higher concentrations (*φ* = 0.3), the deviations increase. The reason is that there are no data of measurement for the dimensionless parameters *γ*<sup>1</sup> and *α*0. To choose them, we apply the method of least squares based on minimizing the functional *F*(*γ*1, *α*0). As it happened, a discrepancy between the measured data and calculations for different concentrations and angular velocities is due to the optimal choice of these unknown parameters.

The above calculations confirm that Equation (30) for the skew-symmetric viscosity *ηsk*(*φ*) can be of use.

Observe that comparison with experiments for colloidal suspensions is contained in Figures 1, 2 and 5 since, in the case *φ* = 0, the suspension becomes a pure colloidal fluid.

**Figure 1.** Calculated dimensionless velocity profiles (solid lines) versus the radial variable for pure interstitial Herschel–Bulkley fluid without particles, *φ* = 0. The lines from the bottom upward correspond to Ω = 2, Ω = 5 and Ω = 100, respectively. Stars, balls, and triangles stand for measurement data [9] in the cases Ω = 2, Ω = 5 and Ω = 100, respectively.

**Figure 2.** Dimensionless velocity profiles versus the radial variable in pure interstitial Herschel– Bulkley fluid without particles, *φ* = 0, for Ω = 2, 5, 10, 20, 50, 100[rpm] from the bottom upwards. (**a**) calculations, (**b**) measured data [9] .

**Figure 3.** Dimensionless velocity profiles of the Herschel–Bulkley fluid with the apparent *ηHB*(*φ*) and *τy*(*φ*) for *φ* = 0.3 and Ω = 2, 5, 10, 20, 50, 100[rpm] from the bottom upwards. (**a**) calculations, (**b**) measured data [9].

**Figure 4.** The solid line corresponds to a dimensionless velocity profile versus the radial variable for Ω = 100[rpm]. Dots stand for experimental data [9]. (**a**) *φ* = 0.3, (**b**) *φ* = 0.1.

**Figure 5.** Dimensionless velocity profiles for Ω = 100[rpm]. The curves from the bottom upwards correspond to *φ* = 0, *φ* = 0.1 and *φ* = 0.3. (**a**) Calculations, (**b**) measured data [9].

**Figure 6.** Dimensionless velocity profiles for Ω = 5[rpm]. The curves from the bottom upwards correspond to *φ* =0,0.1 and 0.3. (**a**) calculations, (**b**) measured data [9].

**Figure 7.** The solid line corresponds to dimensionless velocity profile versus the radial variable for Ω = 5[rpm] and *φ* = 0.1. Dots stand for experimental data [9].

**Figure 8.** Calculated dimensionless velocity profiles versus the radial variable for *φ* = 0.3. The lines from bottom upward correspond to Ω = 2[rpm], Ω = 5[rpm] and Ω = 50[rpm], respectively. Triangles, balls and stars are the measured data from [9] corresponding to Ω = 2[rpm], Ω = 5[rpm] and Ω = 50[rpm], respectively.

**Figure 9.** (**a**) Calculated dimensionless velocity profiles for *φ* = 0.3. The curves correspond to Ω = 2, 5, 10, 20, 50, 100[rpm] from the bottom upwards. (**b**) Measured [9] dimensionless velocity profiles of the Herschel–Bulkley fluid with *φ* = 0.3 and Ω taking values 2, 5, 10, 20, 50, 100[rpm] from the bottom upwards.

**Figure 10.** Calculated angular velocity *w*[rpm] profiles for *φ* = 0.3. The curves correspond to Ω = 2, 5, 10, 20, 50, 100[rpm] from the top down.

#### **6. Rotational Sedimentation**

Here, we consider more general mathematical model allowing for non-uniform particle distribution. We introduce the mass concentration of particles as follows:

$$
\mathcal{L} = \frac{\mathfrak{p}\_s \mathfrak{\phi}}{\rho}, \quad \rho = \mathfrak{p}\_s \mathfrak{\phi} + \mathfrak{p}\_f (1 - \mathfrak{\phi}), \tag{65}
$$

where *ρ* is the total density, *ρ*¯*<sup>s</sup>* is the particle density and *ρ*¯*<sup>f</sup>* is the density of the interstitial fluid. Given *c*, one can restore from (65) the volume concentration and the total density by the formulas

$$\phi = \frac{\bar{\rho}\_f c}{\bar{\rho}\_f c + \bar{\rho}\_s (1 - c)}, \quad \rho(c) = \frac{\bar{\rho}\_f \bar{\rho}\_s}{\bar{\rho}\_f c + \bar{\rho}\_s (1 - c)}. \tag{66}$$

Due to these formulas, any given function of the volume concentration like the relative viscosity *ε*(*φ*) can be defined in terms of the mass concentration *c*. It is explained in [7] that *c* satisfies the conservation law

$$\frac{\partial(\rho c)}{\partial t} + \text{div}(\rho c \mathbf{v} + \mathbf{l}) = 0,\tag{67}$$

where **l** is the concentration flux obeying the generalized Fick equation

$$
\rho \mathbf{1} = -D \nabla c - D\_p \nabla p + D\_\omega \mathbf{rot} \,\omega \times \boldsymbol{\omega}\_r. \tag{68}
$$

The scalar parameters *D*[cm2/s], *Dp*[cm<sup>3</sup> · s/g], and *<sup>D</sup><sup>ω</sup>* [cm<sup>2</sup> · s], stand for the diffusion, barodiffusion, and spin diffusion coefficients.

Observe that, instead of (5), the couple stress tensor *N* is prescribed by the rheological equation

$$N = 2\gamma A + \frac{D\_{\omega}}{2}\epsilon : (1 \times \omega\_r) \tag{69}$$

to meet the entropy production law [7] where the skew-symmetric matrix *e* : **a** is defined by the formula

$$(\epsilon : \mathbf{a})\_{\vec{l}\vec{j}} = a\_k \epsilon\_{ik\vec{j}\prime} \quad \mathbf{a} = a\_{\vec{l}} \mathbf{e}\_{\vec{l}\prime}$$

in any orthogonal basis {**e***i*} 3 1 . It is due to spin diffusion that the Ségre-Silberberg effect is explained within the micropolar theory [7].

Due to the identity rot *ω* × **b** = 2(∇*ω*)*<sup>a</sup>* · **b**, ∀**b** ∈ *R* 3 , it follows from (68) and (69) that

$$\left(\rho + \frac{D\_{\omega}^{2}|\omega\_{r}|^{2}}{4\gamma}\right)\mathbf{l} = -D\nabla c - D\_{p}\nabla p - \frac{D\_{\omega}}{2\gamma}\mathbf{N}\_{a} \cdot \omega\_{r}\mathbf{l}$$

in agreement with the definition of rotary diffusion [37]: "Just as the translational diffusion coefficient is calculated in terms of the drag force, so the rotary diffusion coefficient is expressed in terms of the moment of the forces on a particle executing a rotary movement."

It follows from (9) that for steady flows the mass conservation law becomes div**v** = 0. With the above definitions, we arrive at the following conservation laws for steady flows:

$$\operatorname{div}(\rho \mathbf{v} \otimes \mathbf{v}) = -\nabla p + \operatorname{div} \mathcal{S} + \rho \mathbf{f},\tag{70}$$

$$J\text{div}[\omega \otimes (\rho \varepsilon \mathbf{v} + \mathbf{1})] = \text{div}\, N - \varepsilon : \mathbf{S} \,, \tag{71}$$

$$\operatorname{div}(\rho c \mathbf{v} + \mathbf{1}) = 0,\tag{72}$$

with tensors *S*, *N* and the flux **l** given by Equations (51),(69) and (68) respectively.

For the flows in a concentric-cylinder Couette geometry, we calculate that

$$1 = l\mathbf{e}\_{r\prime} \quad \rho l = -D \frac{\partial \mathcal{c}}{\partial r^{\prime}} - D\_p \frac{\partial p}{\partial r^{\prime}} + \frac{D\_{\omega} \omega\_r}{2} \frac{\partial \omega}{\partial r^{\prime}}.$$

Under the assumption that *c* = *c*(*r*), we arrive at the formula ∇*c* = **e***r∂c*/*∂r*. Due to equation div**v** = 0, we obtain that, for the rotation flows, the equation div(*ρc***v**) = 0 holds. Now, it follows from (72) that

$$0 = \operatorname{div} \mathbf{l} = \frac{1}{r} \frac{\partial (rl)}{\partial r} \quad \text{and} \quad rl = \text{const.}$$

At the same time, the now-flow boundary conditions **l** · **n**|*R<sup>i</sup>* = 0 imply that *l* = 0 and **l** = 0. Thus, the particles mass concentration obeys the equation *l* = 0 or

$$\frac{\partial \mathcal{C}}{\partial r} = -\frac{D\_p}{D} \frac{\partial p}{\partial r} + \frac{D\_{\omega} \omega\_r}{2D} \frac{\partial \omega}{\partial r}.$$

Due to (37), the latter equation can be written as

$$\frac{\partial c}{\partial r} = -\frac{D\_p}{D} \frac{\rho v^2}{r} + \frac{D\_\omega}{2D} \left(\omega - \frac{1}{2} \left(\frac{\partial v}{\partial r} + \frac{v}{r}\right)\right) \frac{\partial \omega}{\partial r}.\tag{73}$$

In what follows, we use the representations *D<sup>p</sup>* = *D*∗ *p c*(1 − *c*), *D<sup>ω</sup>* = *D*<sup>∗</sup> *<sup>ω</sup>c*(1 − *c*) since both *D<sup>p</sup>* and *D<sup>ω</sup>* vanish at *c* = 0 and *c* = 1. Given a mean value *c*<sup>0</sup> of *c*, we set the following condition:

$$\int\_{R\_1}^{R\_2} r c(r) \, dr = \frac{c\_0 (R\_2^2 - R\_1^2)}{2}. \tag{74}$$

As for the particle migration, there is one more approach based on the Fick law. This is known as the Lamm equation for a highly disperse colloidal solution enclosed in a wedge-

shape cell rotating at an angular velocity Ω about an axis coinciding with the apex of the wedge [5].

Let us return to Equation (73). One more consequence of the equality **l** = 0 is that the rheological Equation (69) reduces to *N* = 2*γA* in the Couette geometry.

Let us summarize the mathematical model. We look for functions *Srϕ*, *Sϕ<sup>r</sup>* , *Nzr*, *v*, *ω*, *c*, obeying the Equations (35), (36), (38)–(41) and (73), with the given function *ε*(*c*).

We introduce dimensionless parameters

$$
\rho\_0 = \frac{\overline{\rho}\_s}{\overline{\rho}\_f}, \quad D\_1 = \frac{D\_p^\* R\_1^2 \Omega\_0^2 \overline{\rho}\_s}{D}, \quad D\_2 = \frac{D\_\omega^\* \Omega\_0^2}{2D}.
$$

In new notations,

$$
\varepsilon(c) = \left(1 - \frac{c/\phi\_\*}{c + \rho\_0(1 - c)}\right)^{-E\phi\_\*} - 1.
$$

Let us pass to dimensionless variables. Then, Equations (73) and (74) become

$$\frac{1}{c(1-c)}\frac{\partial'c}{\partial r'} = -\frac{D\_1v'^2}{r'(c+\rho\_0(1-c))} + D\_2\frac{\partial'\omega'}{\partial r'}\left(\omega' - \frac{\partial'v'/\partial r' + v'/r'}{2}\right),\tag{75}$$

$$\int\_{1}^{a} r' \mathbf{c}(r') \, dr' = \frac{c\_0(a^2 - 1)}{2}, \quad a = \frac{R\_2}{R\_1}. \tag{76}$$

We summarize the governing equations as follows. We look for the dimensionless functions *v* 0 (*r* 0 ), *ω*0 (*r* 0 ) and *c*(*r* 0 ) which satisfy Equations (43)–(50), (75) and (76). Observe that these equations are not decoupled since the relative viscosity *ε*(*c*) in Equations (44) and (45) depends on the particle concentration *c*.

Figure 11 depicts results of calculations of concentration along the radial variable. We apply the Wolfram Mathematica solver for ordinary differential equations. Agreement between calculations and experiment [9] is achieved for Ω = 10<sup>2</sup> rpm with the choice *D*∗ *<sup>ω</sup>*/(2*D*) = <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> [s 2 ]. The sedimentation effect happens when we increase Ω to <sup>14</sup> <sup>×</sup> <sup>10</sup><sup>3</sup> [rpm]. We prove that such an effect is due to the rotational diffusion *D<sup>ω</sup>* since the particle separation does not happen when *D<sup>ω</sup>* = 0.

**Remark 2.** *It is known that, for many practical purposes, the polymers or colloidal particles can be regarded as rigid particles. Examples are triblock Janus particles which can be modeled as crosslinked polystyrene spheres whose poles are patched with sticky alkyl groups, and their middle band is covered with negative charges [38]. This is why the above results on rotational sedimentation can also be applied to polymer flows. As was proved by Svedberg, such flows are of great importance in the studies of the polymer structure. Our contribution is that we propose an alternative approach to the Lamm equation based on the empirical notion of the sedimentation coefficient [6]. The advantages are that we apply the conservation laws of continuum mechanics and take into account the shape of the particles. Moreover, we prove that the polymer sedimentation is due to its rotation. This result suggests a new direction of laboratory studies on polymer flows.*

**Figure 11.** Profiles of mass concentration *c*. Both the solid line that is based on calculations and dots standing for experimental data [9] correspond to Ω = 10<sup>2</sup> rpm. Agreement between calculations and experiment is achieved by the choice *D*∗ *<sup>ω</sup>*/(2*D*) = <sup>5</sup> <sup>×</sup> <sup>10</sup>−<sup>5</sup> [s 2 ]. The dashed line corresponding to calculations reveals the sedimentation effect when we increase <sup>Ω</sup> to 14 <sup>×</sup> <sup>10</sup><sup>3</sup> [rpm].

#### **7. Discussion**

We address the rotational sedimentation of particles for steady flows of yield stress granular fluids in a concentric-cylinder Couette geometry. Apart from the Lamm equation approach, we do not use the empirical sedimentation coefficient. Instead, we apply conservation laws of the micro-polar equations which allow for particle rotation. We prove that it is due to the rotational diffusion that the particle sedimentation occurs at high angular velocity of the Couette cell inner cylinder. To validate the mathematical model, we perform a comparison with published data of measurements by choosing the relative viscosity related with the particle rotation. First, we justify analytically this choice for dilute suspensions starting from the Einstein correlation for the apparent viscosity. As for dense suspensions, we apply the Krieger–Douhgerty empirical closure for the apparent viscosity. Though we performed calculations for steady flows, the developed approach allows for unsteady flows and non-spherical particles due to the micro-inertia tensor involved into the angular momentum conservation law.

**Funding:** This research was funded by RUSSIAN SCIENCE FOUNDATION Grant No. 20-19-00058.

**Data Availability Statement:** Not applicable.

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

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Polymers* Editorial Office E-mail: polymers@mdpi.com www.mdpi.com/journal/polymers

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-6667-2