Recursive Settling of Particles in Shear Thinning Polymer Solutions: Two Velocity Mathematical Model

Abstract

Processing of the available experimental data on particles settling in shear-thinning polymer solutions is performed. Conclusions imply that sedimentation should be recursive, since settling also occurs within the sediment. To capture such an effect, a mathematical model of two continua has been developed, which corresponds to experimental data. The model is consistent with basic thermodynamics laws. The rheological component of this model is a correlation formula for gravitational mobility. This closure is justified by comparison with known experimental data available for particles settling in vertical vessels. In addition, the closure is validated by comparison with analytical solutions to the Kynch one-dimensional equation, which governs dynamics of particle concentration. An explanation is given for the Boycott effect and it is proven that sedimentation is enhanced in a 2D inclined vessel. In tilted vessels, the flow is essentially two-dimensional and the one-dimensional Kynch theory is not applicable; vortices play an important role in sedimentation.

1. Introduction

Particulate fluids are common in both natural and industrial processes. Fiber-reinforced polymers, detergents, blood, and drilling muds are some examples where particles are present in a suspending fluid [1,2]. Sedimentation of suspensions in complex fluids is of interest in technological operations, such as oil and gas exploration. Many drilling muds are polymer-based fluids with shear-thinning rheology manifesting itself through a loss of apparent viscosity with increasing strain-rates. Operational stops occur frequently during oil-drilling processes. Interruption of drilling mud pumping causes the particles to settle through the annular space. This may give rise some undesirable phenomena in drilling operations, such as a stuck pipe. The understanding of suspension dynamics in complex fluids is also important for applications in particle manipulation in microfluidic devices [3]. The methods and numerical algorithms developed in the present paper contribute to enhancing knowledge of the processes and their optimization in the fields of oil exploration and microfluidic devices.

Laboratory study based on the gamma-ray attenuation technique reveals significant differences in the behavior of particles settling in Newtonian and non-Newtonian fluids [4]. Particularly, control of the concentration of solids proves that sediment is formed faster in shear-thinning fluid when compared to a Newtonian fluid of similar apparent viscosity.

In order to develop a theory of settling, it is first necessary to understand the interaction of a single particle with a carrier fluid and the mutual interaction between two particles.

Clearly, the dynamics of a particles can be strongly affected by the rheology of the interstitial fluid. In Ref. [5], some problems involving rigid spherical particles in shear-thinning fluids are considered in the absence of inertia; for settling particles, analytical formulas are derived and differences in comparison to a Newtonian fluid are demonstrated. Experiments on sedimentation in non-Newtonian shear thinning fluids were considered in Ref. [6] in order to validate a certain empirical dependence of the particle settling velocity on the volume concentration of particles. It should be noted that, based on one-velocity continuum models, the existing tools of computational fluid dynamics make it possible to solve sedimentation problems, taking into account many complex processes. It is worth mentioning the works [7,8,9], where the sedimentation of particles in stirred vessels is analyzed in the case of two-phase turbulent mixing flows, taking into account chemical reactions, crystallization, and dissolution processes.

In this work, a different approach is used. In order to take into account particle–fluid and particle–particle interactions, particles are considered as a phase, which can be described within the continuum mechanics approach. As a result, the entire suspension becomes a two-velocity continuum, with the fluid and the solid phases enjoying some rheological constitutive laws. Such a method was validated in the recent paper concerning the Boycott sedimentation effect, which stated that enhanced sedimentation occurs in a tilted vessel [10].

In a great number of papers [11], particles and the carrier fluid are also assumed to be two different phases. Let us highlight the differences and similarities between our model and other two-phase approaches. First of all, the equations used are consistent with the laws of thermodynamics. The article by A.S. Baumgarten and K. Kamrin [12] also refers to the two-phase model, which is consistent with thermodynamics, but the method and equations are different. To achieve compliance with the laws of thermodynamics in the case of concentrated suspensions, the ideas that were proposed for the mathematical description of the superfluidity of liquid helium II in the papers of L. D. Landau and I. M. Khalatnikov [13,14] are applied.

An important point of the theory is that interaction forces between phases for reversible processes without dissipation effects can be uniquely identified while reconciling the energy conservation equation, both with other conservation laws and the basic thermodynamic principles. As far as dissipative processes are concerned, the interaction forces are introduced, fitting the general deGroot–Mazur principles of irreversible thermodynamics [15]. The Landau–Khalatnikov thermodynamic method finds applications in the description of multi-phase flows [16], particularly in studying fluid-saturated poroelastic media, both with the use of a two-velocity [17] and a three-velocity continuum [18]. Recently, it was established in Ref. [19] that the same approach can be applied in building up conservation laws for suspensions with rotating particles. In the present paper, we extend the Landau–Khalatnikov approach by formulating the diffusion equation for the mass concentration of particles involving a generalized Fick’s law for the concentration flux vector in such a way that it depends not only on the concentration gradient and gravitation vector scaled by a mobility factor, but on the pressure gradient, temperature gradient, and gradient of modulus of the slip velocity as well.

The principle goal of the present paper is to adjust the method of [10] for the description of particles settling in a polymer solution. To this end, particles are assumed to be suspended in a non-Newtonian power-law fluid, in contrast to [10], where such a fluid is considered to be Newtonian. One more rheological feature of this work is that the gravitational mobility in Fick’s law is chosen in such a way that the numerical results on sedimentation are consistent with experiments on particles settling in polymer solutions [4]. In fact, a suitable correlation is proposed for the dependence of mobility on particle concentration. The proposed arguments are based on the Kynch theory, in which the sedimentation of a dense suspension is considered as a concentration wave [20]. A zone with a high concentration of particles is first formed at the bottom of a vessel; then, this zone propagates upward like a shock wave. Some successful applications of this theory can be found in Refs. [21,22]. The calculations performed in this work show that the sedimentation of particles in polymer solutions is recursive, since settling is also observed inside the sediment. Sedimentation in a tilted vessel is considered and it is proven that the Boycott effect also holds in the polymer solution. A comparison with Newtonian carrier fluid is performed.

2. Mathematical Model

Let us consider a joint flow of two continua when an arbitrary volume V contains a fluid (index f) and a granular phase (index s). Volume, mass, and pressure of the fluid and the granular phases are denoted by $V_f,m_f,p_f$, and $V_s,m_s,p_s$, respectively. It is assumed that the granular phase is a mixture of dry particles and a carrier fluid such as proppant and gel, for example. In this case, $V_s=V_M+V_p$ and $m_s=m_M+m_p$. The particles are “frozen” in the carrier liquid; i.e., the granular phase is characterized by just one speed ${\mathbf{v}}_s$, one viscosity, and one stress tensor. In what follows, the indexes f and s stand for fluid and solid phases, respectively.

The quantities

$$\rho =\frac{m}{V},{\rho }_s=\frac{m_s}{V},{\rho }_f=\frac{m_f}{V},{\rho }_p=\frac{m_p}{V},{\varphi }_j=\frac{V_j}{V},{\rho }_M=\frac{m_M}{V},c=\frac{m_p}{m}.$$

(1)

are assigned to the unit volume. Here, $c={\rho }_p/\rho$ is the particle mass concentration and ${\varphi }_j$ is the volume fraction of the j-phase with $j=f,p,M$. It follows from the above definitions that the partial densities ${\rho }_j$ are related to the material densities ${\overline{\rho }}_j$ by the following formulas

$${\rho }_j={\varphi }_j{\overline{\rho }}_j,{\overline{\rho }}_j\equiv \frac{m_j}{V_j},{\varphi }_f+{\varphi }_s=1,{\varphi }_s={\varphi }_p+{\varphi }_M,\rho ={\rho }_f+{\rho }_p+{\rho }_M.$$

Generally, the phase pressures $p_s$ and $p_f$ are different. However, as in Ref. [19], it is assumed that $p_s=p_f=p$. Such a hypothesis works well when the surface tension at the boundaries separating the phases is negligible.

Let ${\mathbf{v}}_i$, $T_i$, $\mathbf{l}$, k stand for the velocity, the viscous part of the stress tensor, the particle concentration flux vector, and the interface friction coefficient, respectively.

In what follows, we use the tensor notations. Given two vectors $\mathbf{a}$ and $\mathbf{b}$, we define the scalar product $\mathbf{a}·\mathbf{b}=$ $a_i{b}_i$. The tensor product $\mathbf{a}\otimes \mathbf{b}$ is a matrix such that ${(\mathbf{a}\otimes \mathbf{b})}_{ij}=$ $a_i{b}_j$. The matrix $A^{\ast }$ stands for the adjoint matrix of A, i.e., ${\left(A^{\ast }\right)}_{ij}={\left(A\right)}_{ji}$. The i-th component of the vector $\mathrm{div}A$ is defined by the formula ${\left(\mathrm{div}A\right)}_i=\partial A_{ik}/\partial x_k$.

Neglecting the rotation of particles and thermal effects, it follows from [19] that the mathematical model

$$\frac{\partial \left({\rho }_s{\mathbf{v}}_s\right)}{\partial t}+\mathrm{div}({\rho }_s{\mathbf{v}}_s\otimes {\mathbf{v}}_s)=-\frac{{\rho }_s}{\rho }\nabla p-\frac{{\rho }_s{\rho }_f}{2\rho }\nabla {\mathbf{u}}^2-k\mathbf{u}+\mathrm{div}{T}_s+{\rho }_s\mathbf{g},$$

(2)

$$\frac{\partial \left({\rho }_f{\mathbf{v}}_f\right)}{\partial t}+\mathrm{div}({\rho }_f{\mathbf{v}}_f\otimes {\mathbf{v}}_f)=-\frac{{\rho }_f}{\rho }\nabla p+\frac{{\rho }_s{\rho }_f}{2\rho }\nabla {\mathbf{u}}^2+k\mathbf{u}+\mathrm{div}{T}_f+{\rho }_f\mathbf{g},$$

(3)

$$\frac{\partial \left(\rho c\right)}{\partial t}+\mathrm{div}(c\mathbf{j}+\mathbf{l})=0,$$

(4)

$${\rho }_{st}+\mathrm{div}\left({\rho }_s{\mathbf{v}}_s\right)=0,{\rho }_{ft}+\mathrm{div}\left({\rho }_f{\mathbf{v}}_f\right)=0,$$

(5)

for six unknown functions ${\rho }_s$, ${\rho }_f$, p, c, ${\mathbf{v}}_s$, ${\mathbf{v}}_f$ can be derived. Here, $p=p\left(\rho \right)$ is the prescribed state equation, $\mathbf{g}$ is the gravitation vector, and

$${(\nabla p)}_i=\frac{\partial p}{\partial x_i},\mathbf{j}={\rho }_s{\mathbf{v}}_s+{\rho }_f{\mathbf{v}}_f,\rho ={\rho }_s+{\rho }_f,\mathbf{u}={\mathbf{v}}_s-{\mathbf{v}}_f,\mathrm{div}\mathbf{l}\equiv \frac{\partial l_i}{\partial x_i},{\mathbf{u}}^2=\mathbf{u}·\mathbf{u}.$$

Rheological assumptions are formulated as follows. Given a velocity field $\mathbf{v}\left(x\right)$, the corresponding rate of strain tensor D is defined by the formula

$$D=\frac{\nabla \mathbf{v}+{(\nabla \mathbf{v})}^{\ast }}{2},{(\nabla \mathbf{v})}_{ij}=\frac{\partial v_i}{\partial x_j}.$$

The fluid phase is considered to be a non-Newtonian fluid. This implies that

$$T_f=2{\eta }_f{D}_f,$$

(6)

with the power-law viscosity

$${\eta }_f={\eta }_f^0{\dot{\gamma }}^{n-1},$$

(7)

where the power n satisfies the shear thinning condition $0<n<1$. Here, ${\eta }_f^0$ is the consistency and $\dot{\gamma }$ is the dimensionless shear strain:

$$\dot{\gamma }=\frac{\sqrt{2D_f:D_f}}{\omega },D:D=D_{ij}{D}_{ij},$$

(8)

with $\omega$ being the reference frequency. The shear stress $\tau =\sqrt{T_f:T_f/2}$ satisfies the equality

$$\tau =\omega {\eta }_f^0{\dot{\gamma }}^n.$$

(9)

Given the volume fraction of the solid phase ${\varphi }_s$, the rheology of the solid phase is defined by the Newtonian law

$$T_s\left({\varphi }_s\right)=2{\eta }_s\left({\varphi }_s\right)D_s.$$

(10)

Here,

$${\eta }_s\left({\varphi }_s\right)={\eta }_s^0{\left(1-\frac{{\varphi }_s}{{\varphi }_s^{\ast }}\right)}^{-2.5}$$

(11)

is the viscosity given by the Krieger–Douhgerty empirical closure [23], with ${\varphi }_s^{\ast }$ and ${\eta }_s^0$ being the maximal reference value of ${\varphi }_s$ and the consistency, respectively.

One more rheological equation is the Fick law [10]:

$$\mathbf{l}=-\left({\gamma }_3\nabla c+{\gamma }_1\nabla p+{\gamma }_4\nabla {\mathbf{u}}^2\right)+\rho cB\mathbf{g}.$$

(12)

Due to the mass conservation laws (5), Equations (2) and (3) reduce to

$${\rho }_s\left(\frac{\partial {\mathbf{v}}_s}{\partial t}+{\mathbf{v}}_s·\nabla {\mathbf{v}}_s\right)=-\frac{{\rho }_s}{\rho }\nabla p-\frac{{\rho }_s{\rho }_f}{2\rho }\nabla {\mathbf{u}}^2-k\mathbf{u}+\mathrm{div}{T}_s+{\rho }_s\mathbf{g},$$

(13)

$${\rho }_f\left(\frac{\partial {\mathbf{v}}_f}{\partial t}+{\mathbf{v}}_f·\nabla {\mathbf{v}}_f\right)=-\frac{{\rho }_f}{\rho }\nabla p+\frac{{\rho }_s{\rho }_f}{2\rho }\nabla {\mathbf{u}}^2+k\mathbf{u}+\mathrm{div}{T}_f+{\rho }_f\mathbf{g},$$

(14)

where ${\left(\mathbf{v}·\nabla \mathbf{v}\right)}_i=v_j\partial v_i/\partial x_j.$ Such equations are of use in the numerical calculations performed below.

3. Incompressible Isothermal Flows

Let us formulate a hypothesis of incompressibility. It is assumed that the mud volume fraction ${\varphi }_M$ is negligible and the densities of materials ${\overline{\rho }}_f$, ${\overline{\rho }}_p$ and ${\overline{\rho }}_M$ are constants. Then, it follows from (1) that

$${\rho }_s\approx c\rho ,{\rho }_f\approx (1-c)\rho ,{\varphi }_s=\frac{c}{R_0+c(1-R_0)},{\varphi }_f=\frac{R_0(1-c)}{R_0+c(1-R_0)},$$

(15)

where $R_0={\overline{\rho }}_s/{\overline{\rho }}_f$. Observe that the total density $\rho$ and the partial densities ${\rho }_j$ are not constant in contrast to the densities of the materials. By the incompressibility assumption, one can easily derive the following formulas:

$$\frac{{\rho }_s}{{\overline{\rho }}_f}=c[1+(R_0-1){\varphi }_s\left(c\right)]\equiv r_s\left(c\right),\frac{{\rho }_f}{{\overline{\rho }}_f}=1-{\varphi }_s\left(c\right)\equiv r_f\left(c\right),\rho =\frac{{\overline{\rho }}_s}{R_0+c(1-R_0)}.$$

The functions $r_s\left(c\right)$ and $r_f\left(c\right)$ are dimensionless partial densities.

One more consequence of the incompressibility assumption is that the volumetric mean velocity is divergence-free:

$$\mathrm{div}\mathbf{v}=0,\mathbf{v}\equiv {\varphi }_s\left(c\right){\mathbf{v}}_s+{\varphi }_f\left(c\right){\mathbf{v}}_f.$$

(16)

Equation (4) is equivalent to

$$\frac{\partial c}{\partial t}+\tilde{\mathbf{v}}·\nabla c+{\rho }^{-1}\left(c\right)\mathrm{div}\mathbf{l}=0,\tilde{\mathbf{v}}\equiv c{\mathbf{v}}_s+(1-c){\mathbf{v}}_f,$$

(17)

where $\tilde{\mathbf{v}}$ is the mean mass velocity. Thus, a mathematical model for four unknown functions p, c, ${\mathbf{v}}_s$ ${\mathbf{v}}_f$ is derived that obeys the Equations (13)–(17). The parameters ${\eta }_s$, ${\eta }_f$, k, ${\gamma }_j$ are assumed to be known functions of the mass concentration c.

Under the incompressibility hypothesis, pressure is no longer a thermodynamic parameter and does not satisfy the equation of state. It is now included in the list of unknown functions, as in the case of Navier–Stokes models of a viscous incompressible fluid. Densities can be restored from equalities (15).

The diffusion coefficients ${\gamma }_j$ vanish when any phase disappears. As for the friction, we use the correlation formula

$$k\left(c\right)=\frac{3}{4}{C}_D\frac{c{\overline{\rho }}_f\left|\mathbf{u}\right|}{d_p},$$

(18)

proposed in Ref. [24], where $d_p$ is the particle diameter and $C_D$ is the particle/fluid friction:

$$C_D=\left\{\begin{array}{cc}\frac{24}{R{e}_p}\left(1+0.15R{e}_p^{0.678}\right)\hfill & \mathrm{if}R{e}_p<1000,\hfill \\ 0.44\hfill & \mathrm{if}R{e}_p>1000,\hfill \end{array}\right.R{e}_p=\frac{d_p{\overline{\rho }}_f\left|\mathbf{u}\right|}{{\eta }_f^0}.$$

For the case of sedimentation, the particle Reynolds number $R{e}_p$ is very small, $R{e}_p<<1$, and one can use the following approximation $C_D=24/R{e}_p,k=18{\eta }_f^{0}c/\left(d_p^2\right).$

With g being the gravitation acceleration, the formula $\mathbf{g}=$ $-g{\mathbf{e}}_y$ is valid where ${\mathbf{e}}_y={(0,1)}^T$. The domain ${\Omega }_0$ denotes the vertical cell $\{0<x<a_1,0<y<a_2\},$ with the y-axis directed upwards. In what follows, flows are considered in the tilted cell ${\Omega }_{\alpha }$, with the inclination angle $\alpha$ measured from ${\mathbf{e}}_y$.

Given the reference values $a_1$, V, $l_0$, $t_0$, and $\overline{p}$, the dimensionless variables are defined as follows

$$x^{\prime }=\frac{x}{a_1},y^{\prime }=\frac{y}{a_1},{\mathbf{v}}^{\prime }=\frac{\mathbf{v}}{V},p^{\prime }=\frac{p}{\overline{p}},{\mathbf{l}}^{\prime }=\frac{\mathbf{l}}{l_0},t^{\prime }=\frac{t}{t_0},$$

(19)

with the assumptions

$$t_0=\frac{a_1}{V},\omega =\frac{V}{a_1},l_0={\overline{\rho }}_{f}V,\overline{p}={\overline{\rho }}_{f}g{a}_1.$$

As a result of passage to dimensionless variables, the following dimensionless parameters and functions appear:

$$\mathrm{Re}=\frac{a_{1}V{\overline{\rho }}_f}{{\eta }_f^0},k_1=\frac{18{\eta }_f^0{a}_1}{{\overline{\rho }}_{f}V{d}_p^2},\mathrm{Fr}=\frac{g{a}_1}{V^2},$$

(20)

$$\lambda =\mathrm{Fr}·\mathrm{Re}=\frac{{\overline{\rho }}_f{a}_1^{2}g}{{\eta }_f^{0}V},\beta =k_1·\mathrm{Re}=18\frac{a_1^2}{d_p^2},$$

(21)

$${\Gamma }_1\left(c\right)=\frac{{\gamma }_1\left(c\right)\overline{p}}{l_0{a}_1},{\Gamma }_3=\frac{{\gamma }_3}{a_1{l}_0},{\Gamma }_4\left(c\right)=\frac{{\gamma }_4\left(c\right)V^2}{a_1{l}_0},{\Gamma }_5\left(c\right)=\frac{cr\left(c\right)B\left(c\right)g}{V}.$$

(22)

In the new variables, the rheological equations become $\dot{\gamma }=\sqrt{2D_f^{\prime }:D_f^{\prime }}$,

$$T_f=\frac{{\eta }_f^{0}V}{a_1}{T}_f^{\prime },T_f^{\prime }=2{\dot{\gamma }}^{n-1}{D}_f^{\prime },T_s=\frac{{\eta }_f^{0}V}{a_1}{T}_s^{\prime },T_s^{\prime }=2\frac{{\eta }_s}{{\eta }_s^0}\frac{{\eta }_s^0}{{\eta }_f^0}{D}_s^{\prime }.$$

(23)

Omitting the primes, we find that the functions ${\mathbf{v}}_f^{\prime }(x^{\prime },y^{\prime },t^{\prime })$, ${\mathbf{v}}_s^{\prime }(x^{\prime },y^{\prime },t^{\prime })$, $c(x^{\prime },y^{\prime },t^{\prime })$, and $p(x^{\prime },y^{\prime },t^{\prime })$ satisfy the equations

$$r_s\mathrm{Re}\frac{d_s{\mathbf{v}}_s}{dt}=-\frac{r_s\lambda }{r}\nabla p+\frac{{\eta }_{s0}}{{\eta }_f^0}\mathrm{div}\left(\frac{2{\eta }_s}{{\eta }_{s0}}{D}_s\right)-\beta c\mathbf{u}-\frac{r_s{r}_f\mathrm{Re}}{2r}\nabla {\mathbf{u}}^2-\lambda r_s·{\mathbf{e}}_y,$$

(24)

$$r_f\mathrm{Re}\frac{d_f{\mathbf{v}}_f}{dt}=-\frac{r_f\lambda }{r}\nabla p+\mathrm{div}\left(2{\dot{\gamma }}^{n-1}{D}_f\right)+\beta c\mathbf{u}+\frac{r_f{r}_s\mathrm{Re}}{2r}\nabla {\mathbf{u}}^2-\lambda r_f{\mathbf{e}}_y,$$

(25)

$$\mathrm{div}\mathbf{v}=0,\mathbf{v}\equiv {\varphi }_s\left(c\right){\mathbf{v}}_s+{\varphi }_f\left(c\right){\mathbf{v}}_f.$$

(26)

$$\frac{dc}{dt}+\frac{\mathrm{div}\mathbf{l}}{r}=0,\mathbf{l}=-\left({\Gamma }_3\nabla c+{\Gamma }_1\nabla p+{\Gamma }_4\nabla {\mathbf{u}}^2\right)-{\Gamma }_5{\mathbf{e}}_y,$$

(27)

in the domain

$${\Omega }_{\alpha }=\{x=Xcos\alpha +Ysin\alpha ,y=-Xsin\alpha +Ycos\alpha \},$$

(28)

where $0<X<1$, $0<Y<h=a_2/a_1$. Let us formulate boundary and initial conditions:

$$\partial {\Omega }_{\alpha }:{\mathbf{v}}_s=0,{\mathbf{v}}_f=0,\mathbf{l}·\mathbf{n}=0,\nabla p·\mathbf{n}=0,$$

(29)

$$t=0:{\mathbf{v}}_s={\mathbf{v}}_{s0},{\mathbf{v}}_f={\mathbf{v}}_{f0},c=c_0.$$

(30)

The diffusion coefficients ${\gamma }_1\left(c\right)$ and ${\gamma }_4\left(c\right)$ vanish at $c=0$ and $c=1$. This is why it is reasonable to set ${\Gamma }_i={\Gamma }_i^{0}c(1-c),$ where ${\Gamma }_i^0$ are constants, $i=1,4$.

4. Gravitation Mobility of Particles in a Solution of Polymers

First, we consider settling in a Newtonian fluid where $n=1$ and ${\eta }_f={\eta }_f^0$. For simplicity, it is assumed that $\rho =\mathrm{const}$ and that the gravitation diffusion is dominant; i.e., ${\gamma }_i=0$$\forall i$. In dimension variables, it follows from (4) and (12) that concentration obeys the equation

$$\frac{\partial c}{\partial t}+\mathrm{div}\left[c(\tilde{\mathbf{v}}-gB\left(c\right){\mathbf{e}}_y)\right]=0.$$

(31)

On the other hand, some authors [25] apply the equation

$$\frac{\partial c}{\partial t}+\mathrm{div}\left(c{\mathbf{v}}_p\right)=0$$

(32)

while addressing sedimentation. Here, ${\mathbf{v}}_p=\tilde{\mathbf{v}}+{\mathbf{v}}_{slip}$ is the particle velocity and

$${\mathbf{v}}_{slip}=-V_{St}H\left(c\right){\mathbf{e}}_y,V_{St}=\frac{2({\overline{\rho }}_p-{\overline{\rho }}_f)g{(d_p/2)}^2}{9{\eta }_f^0},$$

(33)

where $V_{St}$ is the Stokes settling velocity in Newtonian fluid and $H\left(c\right)$ is the hindered settling function [25]. Comparing (31) and (32), we find that

$$B\left(c\right)=\frac{V_{St}H\left(c\right)}{g}=2({\overline{\rho }}_p-{\overline{\rho }}_f){\left(\frac{d_p}{2}\right)}^2\frac{H\left(c\right)}{9{\eta }_f^0}.$$

(34)

As shown in Section 6, the Richardson–Zaki closure

$$H_{rz}\left(c\right)={(1-c)}^m$$

(35)

for the function $H\left(c\right)$ meets the settling in a Newtonian fluid with m = 4.65 [26].

As for shear thinning non-Newtonian fluids, it is proven in Section 6 that the shear thinning closure

$$H_{sh}\left(c\right)=a{(1-c)}^k+b{c}^{m-1}(1-c)-d{(1-c)}^l.$$

(36)

for the function $H\left(c\right)$ with the special choice (46) of the constants $a,b,d,k,m$ and l fits the experiment data in Ref. [4] better than any hindered settling function of the form (35). In what follows, the functions $B_{rz}\left(c\right)$ and $B_{sh}\left(c\right)$ are defined by Equation (34) with the hindered settling function $H\left(c\right)$ given by Equations (35) or (36), respectively.

5. Settling in 2D Tube

Here, we apply the mobility closure Equation (34) with the hindered settling function (36) to particles settling in polymer solutions governed by the rheology of the shear thinning fluids. Particularly, we consider sedimentation in inclined vessels on the bases of the two-velocity model (24)–(27) and establish the Boycott effect.

First, sedimentation in a vertical vessel is considered. Calculations reveal that the concentration structures are different for the mobilities $B_{rz}$ or $B_{sh}$. Figure 1 is obtained with the use of the mobility $B_{rz}$ and shows that there are two concentration waves smoothed by the diffusion effect, which propagate in opposite directions. The downward wave starts from the top of the vessel and leaves no particles behind, whereas the wave going up describes an increasing zone of sediment. As for the non-Newtonian fluid with the shear thinning mobility $B_{sh}$, there are three waves, as shown in Figure 2. At the initial stage, only two waves moving towards each other are observed. Then, a third wave appears, going from bottom to top. This wave indicates that settling is also observed in the sediment. Thus, sedimentation of particles in polymer solutions is recursive, since settling is also observed inside the sediment. If all dissipative effects can be neglected, the third concentration wave appears from the very beginning; see Figure 3.

A comparison with experiment data [4] is performed in Section 6 in the case of small diffusion coefficients. Figure 4 and Figure 5 confirm qualitative agreement when we simulate settling in the Newtonian fluid starting from Equations (24)–(27) with $n=1$, using the mobility $B_{rz}\left(c\right)$. The main feature of this case is monotonicity in the following sense. Both the experiment data and our calculations reveal that there is a critical height such that at each level above this height, the concentration decreases with time, and at each level below this height, the concentration increases with time.

Simulation of settling in non-Newtonian shear thinning fluids is based on Equations (24)–(27), with $n=0.34$ and $B\left(c\right)=B_{sh}\left(c\right)$. Qualitative agreement is depicted in Figure 6 and Figure 7. In such a case, the experimental data and our calculations show that the monotonicity property is preserved only for sufficiently high or sufficiently low vertical locations. There are intermediate locations in which concentration first increases with time and then decreases.

For both the fluxes $B_{rz}$ and $B_{sh}$, calculations reveal that in spite of no-slip boundary conditions for velocities, horizontal-strata formation in the vertical vessel occurs in the sediment; see Figure 8 and Figure 9.

All the calculations in the present paper are carried out with the use of the open-source PDE Solver FreeFEM++ based on the finite element method. To obey the restriction that the mean volume velocity is divergence-free, the method of artificial compressibility is applied. At each time step, the concentration field is calculated by the Galerkin-characteristics method, which is implemented in FreeFEM++ through the “convect” function. Next, the Navier–Stokes equations are solved to define velocity and pressure. To tackle nonlinearity, iterations are carried out until convergence. Then, a transition is made to the next time step. A weak formulation of the problem and a detailed description of the algorithm are given in Ref. [10]. The ParaView open-source package visualization application is used to visualize the results.

In experiments with blood, A. Boycott (1920) noticed that erythrocyte particles settled faster in an inclined test tube than in a vertical one. Since that time, many attempts have been made to explain this effect in terms of the theory of particle motion in a viscous fluid. For this purpose, additional hydrodynamic forces were introduced, such as the Archimedes force, the Magnus force, etc. Forces were even used that obviously depended on time and on prehistory, such as the Basset–Boussinesq force [27]. In our recent work [10], the Boycott effect was explained without the use of additional hydrodynamic forces, but only due to the gravitational component of the concentration flux in Fick’s law. As far as we know, the issue of particle settling in non-Newtonian fluids has not been considered theoretically and experimentally in the case of inclined vessels.

Let us consider the issue of particle deposition in a 2D tilted tube within the model (24)–(27). Note that this model was used in Ref. [10] in a particular case under the rheological assumptions that $B=B_{rz}$ and $n=1$ in (7). Now, we set $B=B_{sh}$. As in Ref. [4], calculations are performed for the case $n=0.34$. The data $B_{rz}$ and $n=1$ correspond to a Newtonian fluid, whereas the data $B_{sh}$ and $n=0.34$ correspond to a non-Newtonian shear thinning fluid.

Figure 8 and Figure 9 depict calculated snapshots of concentration in the tilted 2D tube with the inclination angle $\alpha ={30}^{°}$. The above clear area extends differently in the Newtonian and non-Newtonian cases. First, the settling is faster in the Newtonian fluid than in the non-Newtonian one.

Let us introduce the reduced volume of the clear fluid region

$$V\left(t\right)=\frac{1}{a_1{a}_2}\underset{{\Omega }_{\alpha }}{\int }{1}_{c(x,y,t)=0}dxdy,$$

(37)

where $1_A(x,y)$ is the characteristic function of the set A. Figure 10 corresponds to the Newtonian fluid and shows how $V\left(t\right)$ depends on time for the vertical and inclined cells. Enhanced sedimentation is observed due to inclination in agreement with the Boycott effect [28]. The result of inclination in the case of the non-Newtonian fluid is shown in Figure 11. This implies that the Boycott effect also takes place in the non-Newtonian shear thinning fluid.

Streamlines of the mean volume velocity $\mathbf{v}$ of the non-Newtonian fluid, which define transport of the concentration c in the vertical cell, are shown in Figure 12 for some time instant. There are two identical macro-vortices rotating in opposite directions. Due to the lateral particle migration, the lower boundary of the clear fluid zone c = 0 is horizontal at any time. In the case of the tilted cell, one of these vortices becomes dominated and it is in excess of each vertical vortices. Such vortex pattern was observed in experiments by Kinosita [29] and Hill et al. [30].

The recent paper [10] is dedicated to the settling of particles in a Newtonian fluid inside a 2D tube. In this paper, the comparison with experimental data and calculations is given in great detail. One of the results of this work is an explanation of the Boycott effect occurring in the tilted vessel.

6. Kinematic Sedimentation Equation

One can process the Moreira [4] experimental data in Figure 6, which show the dependence of particle concentration on time at different heights of a vertical vessel. As a result of the graphical work of converting one data set to another, concentration profiles can be obtained at various time instants. It turns out that the profiles have a smoothed three-wave structure; one wave goes from top to bottom, and the other two from bottom to top. It should be recalled here that there is a classical Kynch sedimentation theory based on a one-dimensional equation for concentration waves [20]. In this approach, only vertical velocities are taken into account, while transverse velocities are considered negligible. The theory has many applications, including sedimentation with absorption [22]. The Kynch equation can be explained in different ways. In fact, it can be derived from the 3D-system (24)–(27). To do this, it suffices to neglect dissipative effects.

Indeed, let us address vertical flows along the variable y, $0<y<h$, in the case of neutrally buoyant particles. It results from the assumption ${\overline{\rho }}_s={\overline{\rho }}_f$ that the mean volume velocity and the mean mass velocities coincide. Because of the equation $\mathrm{div}\tilde{\mathbf{v}}=0$, the vertical component $\tilde{v}$ of $\tilde{\mathbf{v}}$ satisfies the equation $\partial \tilde{v}/\partial y=0$. At the same time, $\tilde{v}=0$ at the bottom $y=0$. Hence, $\tilde{v}=0$ and Equation (31) reduces to the Kynch equation

$$\frac{\partial c}{\partial t}+\frac{\partial F\left(c\right)}{\partial y}=0,\mathrm{with}F=-gcB\left(c\right),$$

(38)

where $gB=-v_{slip}$ and $v_{slip}=-V_{St}H\left(c\right)$. The solutions of the Kynch equation (38) are determined mainly by the empirical gravitational mobility $B\left(c\right)$. A review of the literature shows that the theoretical three-wave structure can be obtained with a special choice of $B\left(c\right)$, [31,32,33].

To describe such a choice, we first consider settling in a Newtonian fluid. Normally, the slip velocity should be determined experimentally. The Richardson–Zaki correlation

$$v_{slip}=-V_{St}{H}_{rz}\left(c\right),H_{rz}\left(c\right)={(1-c)}^m,$$

(39)

meets experiments with $m=4.65$ in the case of Newtonian fluid flows with small terminal Reynolds numbers $R{e}_t<<1$ [26].

Let us set $V_{St}=1$ for simplicity. This condition is equivalent to switching to another time scale. Given a vanishing diffusion $\epsilon$, the settling problem reduces to the following initial boundary value problem:

$$\frac{\partial c}{\partial t}+\frac{\partial F\left(c\right)}{\partial y}=\epsilon \frac{{\partial }^{2}c}{\partial y^2},F=-cH\left(c\right),0<y<1,$$

(40)

with initial and boundary conditions

$${c|}_{t=0}=c_0{,l|}_{y=0,y=1}=0,$$

(41)

where $l=\epsilon c_y-c{v}_{slip}$. For a more general class of functions $F\left(c\right)$ than $F_{rz}\left(c\right)$, we refer the reader to [21] for the study of Equation (40).

First, we consider the flux $F_{rz}\left(c\right)$ and perform calculations of the concentration profile versus the vertical variable y, $0<y<1$, at different time instances. The boundary conditions

$${c|}_{y=0}=c_{bot}{,c|}_{y=1}=c_{top}$$

(42)

are alternative to (41). The results in Figure 13 and Figure 14 are obtained for the conditions (42). One can observe two discontinuity waves propagating towards each other. There are no particles left behind a wave moving from above. The upward wave describes an increase in the height of the sediment. The data for Figure 13 and Figure 14 differ only in the bottom value of concentration. However, the corresponding concentration waves have different structures. When $c_{bot}=0.45$, the rising shock-wave starting from the bottom is followed by a center rarefaction wave; see Figure 13. This is not the case when $c_{bot}=1$; see Figure 14. Concluding the discussion of sedimentation with the function $F_{rz}\left(c\right)$, we note that experiments on particles settling in Newtonian fluids sometimes lead to a three-wave structure [31]. This implies that the flux $F_{rz}\left(c\right)$ is not a unique choice for settling in Newtonian fluids.

Figure 4 and Figure 5 show good qualitative agreement of our calculations with available experiment data for Newtonian fluids [4], provided we use the flux $F_{rz}\left(c\right)$ in the boundary value problem (40), (41). Each curve in these figures corresponds to a time variation of concentration at a fixed value of the vertical variable. Although a quantitative comparison with these experiments is not possible due to the 1D assumption of the present simulations, the results of calculations generally reproduce the experimental trends.

According to the experiment data in the case of a Newtonian fluid (Figure 4), sedimentation in the 30.1 cm-long tube proceeds monotonously in the following sense. The mean concentration at the level y decreases in time at each point y above the height $y\simeq$ 6 cm and it increases in time at each point y below the height $y\simeq$ 6 cm. The settling loses this monotonicity property as far as the polymer solutions are concerned; see Figure 6. Indeed, one can see that there is an intermediate vertical layer $h_1<y<h_2$ such that the concentration $c(y,t)$ first increases and then decreases in time for the height y from this layer.

To address such issues of monotonicity in the case of shear thinning fluids, we argue as in Ref. [21] and look for the settling velocity correlation in the form

$$v_{slip}\left(c\right)=-V_{St}{H}_{sh}\left(c\right),$$

(43)

$$F_{sh}\left(c\right)=-V_{St}c{H}_{sh}\left(c\right).$$

(44)

with

$$H_{sh}\left(c\right)=a{(1-c)}^k+b{c}^{m-1}(1-c)-d{(1-c)}^l.$$

(45)

It is assumed for simplicity that $V_{St}=1$. Applying the least squares method for minimization of a discrepancy functional, we find that the values

$$a=b=1,d=0.95,n=4,m=10,l=8$$

(46)

meet the data in Figure 6 related to the shear thinning fluids [4]. The discrepancy functional determines how close the data are in Figure 6 and Figure 7. Figure 7 is based on solving the initial boundary value problem (40), (41) with $F=F_{sh}$. Remember that the restriction $V_{St}=1$ is equivalent to switching to another time scale. Clearly, the data (46) depend on the time scale.

We performed calculations of concentration profile versus the vertical variable y, $0<y<1$, at different time instances; see Figure 3. Contrary to the case with the flux $F_{rz}\left(c\right)$, one can observe three discontinuity waves. One wave goes top-down, and the other two rise up, following one after the other. This implies that one more sedimentation occurs within the sediment.

It is worth to remark that loss of monotonicity is also observed for the flux $F_{rz}$. Figure 15 and Figure 16 show that, for some initial data, there is monotonicity, but for other initial data, this property disappears. However, the flux $F_{sh}$ is preferred for simulation of settling in non-Newtonian fluids, as it is consistent with the effect of two-fold particle sedimentation.

The physical meaning of violation of the monotonicity property is as follows. There is a middle layer of the vessel where, at first, the concentration increases due to the sediment going from bottom to top. Then, the concentration in this layer decreases due to a sufficiently pure liquid going from top to bottom. Thus, this effect is the result of the interaction of waves traveling in opposite directions.

Let us summarize the arguments formulated in favor of the new choice (45) for the gravitational mobility $B\left(c\right)=V_{St}{H}_{sh}\left(c\right)/g$. First of all, our method has a physical background. For a wide variety of coefficients, the equation (45) guarantees the experimentally observed three-wave structure of the concentration profile. From top to bottom, there is a discontinuity wave, followed by a rarefaction wave. From the bottom-up, there are two waves of discontinuity. In addition, the choice of (45) is consistent with the violation of the monotonicity property, which is also observed in experiments. An optimal quantitative choice of coefficients in Equation (45) is achieved by minimizing the discrepancy functional. By comparing the data in Figure 6 and Figure 7, this functional determines how experimental and calculated data are close.

7. Conclusions

A two-continua mathematical model for description of particles settling in a shear thinning polymer solution is developed. A correlation formula for the gravitational mobility in Fick’s law is derived and numerical results on sedimentation are consistent with experiments on settling of particles in polymer solutions. In fact, we justify a suitable correlation for the dependence of mobility on particle concentration. The main arguments are based on the Kynch (1952) theory, in which sedimentation of dense suspension is considered as a concentration wave. A zone with a high concentration of particles is first formed at the bottom of a vessel; then, this zone propagates upward like a shock wave. From our calculations, sedimentation of particles in polymer solutions is recursive, since settling is also observed inside the sediment. Sedimentation in a tilted vessel is addressed and it is proven that the Boycott effect also holds in a polymer solution. Comparison with Newtonian carrier fluid is performed.