2. Problem formulation

Following the same assumptions as Bearon et al. (Reference Bearon, Bees and Croze2012), Croze et al. (Reference Croze, Sardina, Ahmed, Bees and Brandt2013, Reference Croze, Bearon and Bees2017) and Fung et al. (Reference Fung, Bearon and Hwang2020), in the present study, we are only concerned with the axisymmetric and axially uniform solution of the suspension in an infinitely long vertical pipe. The flow variables, therefore, consist of the streamwise (downward) velocity $U(r;t)$ and cell concentration $N(r;t)$ as functions of the radial position $r$ and the time $t$. Following Fung et al. (Reference Fung, Bearon and Hwang2020), the equations of motion are given by

(2.1a)\begin{gather} \frac{\partial U}{\partial t}=- G+\frac{1}{Re}\frac{1}{r}\frac{\partial }{\partial r} {\left(r \frac{\partial U}{\partial r} \right)}+Ri (N-1), \end{gather}

(2.1b)\begin{gather} \frac{\partial N}{\partial t}=-\frac{1}{r}\frac{\partial }{\partial r}(r N \langle e_{r} \rangle)+\frac{1}{D_R}\frac{1}{r}\frac{\partial }{\partial r}\left(r {D}_{rr} \frac{\partial N}{\partial r}\right), \end{gather}

with the boundary conditions and compatibility conditions:

(2.1c)\begin{equation} U(1) = 0, \quad \left. \left[\langle e_{r} \rangle N-\frac{{D}_{rr}}{D_R}\frac{\partial N}{\partial r} \right] \right|_{r=1} = 0, \quad \left. \frac{\partial U}{\partial r} \right|_{r=0} = 0, \left. \quad \frac{\partial N}{\partial r} \right|_{r=0} = 0. \end{equation}

Here, the equations have been non-dimensionalised by the radius of the pipe $h^*$, the cell swimming speed $V_s^*$ and the averaged cell concentration $N^*$. The dimensionless parameters of interest are

(2.2a–d)\begin{equation} Ri= \frac{\displaystyle N^* v^* g'^*h^*}{\displaystyle V_c^{*2}}, \quad Re=\frac{\displaystyle V_c^*h^*}{\displaystyle \nu^*}, \quad \lambda=\frac{\displaystyle 1}{\displaystyle 2 B^* D_R^*}, \quad D_R=\frac{D_R^* h^*}{V_c^*}, \end{equation}

where $Re$ is the Reynolds number, $Ri$ the Richardson number, $\lambda$ the inverse of the gyrotactic time scale normalised by rotational diffusivity and $D_R$ the dimensionless rotational diffusivity, sometimes also known as the swimming Péclet number. Here, $\nu ^*$ is the kinematic viscosity of the fluid, $v^*$ the cell volume, $g'^*$ the reduced gravity, $B^*$ the gyrotactic time scale and $D_R^*$ the isotropic rotational diffusivity, which models the randomness in the rotational motion of swimmers. Also, $G={\partial p/\partial z -Ri}$ is the dimensionless driving pressure gradient that excludes the hydrostatic pressure from the negative buoyant cells. As mentioned in § 1, the net cell swimming $\langle e_{r} \rangle$ and effective cell diffusivity $D_{rr}$ in the radial direction are modelled using the GTD theory. Hence, $\langle e_{r} \rangle$ and $D_{rr}$ are functions of the local shear rate $S=-\partial _r U /D_R$, the strength of gyrotaxis ($\lambda$) and rotational diffusion ($D_R$). The shear rate $S$ is sometimes known as the rotary Péclet number (Hinch & Leal Reference Hinch and Leal1972). Since the total number of cells is preserved over the given control volume, we impose the normalisation condition for the cell concentration,

(2.3)\begin{equation} \int_0^1 N(r;t) r\,\textrm{d}r=\tfrac{1}{2}. \end{equation}

Unlike a stationary suspension, the current set-up allows for a net axial flow rate. Therefore, there is an extra degree of freedom, which can be specified either by the pressure gradient $G$ or by the flow rate $Q$ given by

(2.4)\begin{equation} \int_0^1 U(r;t) r\,\textrm{d}r=\frac{Q}{2 {\rm \pi}}. \end{equation}

If $Q$ is prescribed, then the pressure gradient is obtained as $G=(2/Re) \partial _r U|_{r=1}$.

Lastly, the steady version of (2.1) can be further simplified into

(2.5a)\begin{equation} N(r)= N(0) \exp \left(D_R \int_0^r \frac{\langle e_{r} \rangle}{D_{rr}}\,\textrm{d}r\right) \end{equation}

and

(2.5b)\begin{equation} -\frac{1}{Re}\mathcal{D}^2{U}= Ri \left(N(0) \exp \left(D_R \int_0^r \frac{\langle e_{r} \rangle}{D_{rr}} \,\textrm{d}r\right)-1\right)- G, \end{equation}

where $N(0)$ is determined by (2.3) and $G$ is determined by (2.4). Here, we denote the radial Laplace operator by $\mathcal {D}^2=(1/r) \partial _r (r \partial _r)$.

3. Numerical solutions

Following the same approach as Fung et al. (Reference Fung, Bearon and Hwang2020, § 2.6), we solve for the steady solution to (2.1) numerically, where the $r$ dependence is discretised using a Chebyshev collocation method. After that, the resulting algebraic equations are solved using a Newton–Raphson method. The computation is performed with the number of nodes $N_r=250$. Numerical convergence is achieved, as the numerical results only differ from those of $N_r=175$ by $0.2\,\%$ on average. We then perform a pseudo-arclength continuation of the solutions by varying $Ri$ but at a fixed flow rate $Q$. The values of all parameters are taken from table 2 of Fung et al. (Reference Fung, Bearon and Hwang2020), which will not be repeated here for brevity. In Fung et al. (Reference Fung, Bearon and Hwang2020, § 3.1 and figure 2), we have briefly shown that there exist multiple branches of solutions at each $Q$ as long as $Q$ is near zero. Here, the same bifurcation diagram, but with a larger range of $Ri$ and $Q$, is presented. In figures 1(a) and 2, the axial velocity $U(0)$ is used to represent the state of the steady solutions. Figure 1(a) shows how the solution from (2.5) changes with $Ri$ at $Q=0$, while figure 2 shows the same at other $Q$ in the range $[-2.5,2.5]$. We have also plotted some of the steady solutions $U(r)$ when $|U(0)|=1$ and $Q=0$ in figure 1($b$–$i$, blue lines), for reasons that will become apparent later in § 4.1. The respective value of $Ri$ for each $U(r)$ can be found in each panel of figure 1($b$–$i$), as well as figure 1(a).

Figure 1.

Comparison between numerical solutions of (2.5) and the linear and weakly nonlinear analysis in § 4. $(a)$ Bifurcation diagram ($U(0)$ against $Ri$) of the $Q=0$ solutions from (2.5). Here, the line type represents the stability of the solution: , stable; , unstable. The bifurcation point ($Ri_{c,n}$) from (4.5) is denoted by magenta crosses ($\times$), while the slope of the bifurcation curve calculated from (4.8) is indicated by the short magenta dot-dashed segment (). The dotted blue lines (), representing $U(0)=\pm 1$, are added to locate the $Ri$ values for ($b$–$i$). ($b$–$i$) Comparisons between the nonlinear steady solution (blue lines) and the corresponding (appropriately normalised) linear instability mode $\hat {u}(r)$ at $Q=0$ and $Ri$ as indicated in $(a)$.

Figure 2.

Contours of $U(0)$ against $Ri$ at each given $Q \in [-2.5,2.5]$. Here, each contour line in the $U(0)$–$Ri$ plane indicates the value of $Q$: black, $Q=0$; blue to green, $Q \in [0.1,2.5]$ with $0.2$ increment; red to yellow, $Q \in [-2.5,-0.1]$ with $-0.2$ increment. The line type represents the stability of the solution: , stable; , unstable. The thick dotted black lines () show the bifurcation by varying $Q$ for given $G=0$ (see § 5).

Focusing on $Q=0$, from figure 1(a), it is apparent that there exist multiple branches of steady solutions other than the uniform-suspension solution represented by $U(0)=0$. In fact, as we increase $Ri$, there seems to be a countable but infinite number of branches emerging via a sequence of transcritical bifurcations. Each branch coincides with the uniform suspension at a certain $Ri$, which forms the bifurcation point. We shall define $Ri_{c,n}$ as the Richardson number of each bifurcation point from the left in figure 1(a), where $n$ is the index. For convenience, we shall also index the $Q=0$ branches accordingly too. For example, the first bifurcation at $Ri\ (=Ri_{c,1} \approx 190)$ is connected to the first $Q=0$ branch (see the black line crossing the leftmost $Ri_c$ in figure 1a). At each bifurcation point, there is also an exchange of stability, as we examine the stability of the solution near each bifurcation (see § 4.1). However, except for the first bifurcation, the stability exchange takes place not with the most unstable mode but with a less unstable one. When a small downflow or upflow is applied (i.e. $Q \neq 0$), each transcritical bifurcation point turns into a saddle-node point via an imperfect bifurcation. This is similar to the finding of Fung et al. (Reference Fung, Bearon and Hwang2020) for the first branch solution, but the same happens to all the other branches.

4. The origin of an infinite number of solution branches

In this section, we will examine the linear stability of a uniform suspension in the vertical pipe. We will further restrict the perturbation to be axisymmetric, parallel and axially uniform, given the nature of the solutions of interest. To be consistent with our numerical results in § 3 and figure 1(a), in this section we will fix the flow rate when adding the perturbation, i.e. $Q=0$. We will demonstrate that the multiple transcritical bifurcations of the computed solutions are the result of the pipe geometry, which also confines the suspension in a domain with finite horizontal extent. We will then extend it to the weakly nonlinear regime, similarly to previous work by Bees & Hill (Reference Bees and Hill1999). The resulting amplitude equation demonstrates that all the bifurcations in figure 1(a) are indeed transcritical.

4.1. Linear stability analysis

We first consider a perturbation to the stationary uniform suspension in a cylindrical pipe with infinite depth, i.e. $U=\epsilon u_1(r,t)+\epsilon ^2 u_2(r,t)+O(\epsilon ^3)$ and $N=1+\epsilon n_1(r,t)+\epsilon ^2 n_2(r,t)+O(\epsilon ^3)$. Given that $\langle e_{r} \rangle (S)$ is odd and $D_{rr}(S)$ is even with respect to the radial shear rate $S$, this yields

(4.1a,b)\begin{align} \langle e_{r} \rangle=\epsilon\beta \frac{\partial u_1}{\partial r}+\epsilon^2\beta \frac{\partial u_2}{\partial r}+{\textit{O}}(\epsilon^3) \quad \!\text{and}\! \quad \frac{D_{rr}}{D_R}=D+\left.\frac{\epsilon^2}{2 D_R^3}\frac{\partial^2 D_{rr}}{\partial S^2}\right|_{S=0} \left(\frac{\partial u_1}{\partial r}\right)^2+{\textit{O}}(\epsilon^3), \end{align}

where $\beta =-(\partial _S \langle e_{r} \rangle |_{S=0})/D_R$ and $D=D_{rr}|_{S=0}/D_R$ are constants that depend only on $D_R$ and $\lambda$. We also note that $\beta /D=\lambda /2$ (see Bearon et al. Reference Bearon, Bees and Croze2012). At $\textit{O}(\epsilon )$, the perturbed equations for linear stability are then obtained as

(4.2a)\begin{equation} \frac{\partial u_1}{\partial t}=\frac{1}{Re}\mathcal{D}^2 u_1+Ri\ n_1- {G_1} \end{equation}

and

(4.2b)\begin{equation} \frac{\partial n_1}{\partial t}=-\beta \mathcal{D}^2 u_1+D \mathcal{D}^2 n_1, \end{equation}

with boundary conditions

(4.2c)\begin{equation} u_1(1) = 0, \quad \left. \left[\beta \frac{\partial u_1}{\partial r}-D\frac{\partial n_1}{\partial r} \right] \right|_{r=1} = 0. \end{equation}

Because we have fixed $Q=0$, the perturbed pressure gradient $G_1$ is found such that the flow rate is not altered by the perturbation, i.e.

(4.2d)\begin{equation} \int_0^1 u_1 r\,\textrm{d}r=0. \end{equation}

While (4.2) can be solved numerically, we shall proceed to focus on the special case when one of the stability modes is neutrally stable. Introducing a neutrally stable and stationary normal mode (i.e. $u_1(r,t)=\hat {u}(r)$ and $n_1(r,t)=\hat {n}(r)$), (4.2) is then simplified into a single equation:

(4.3)\begin{equation} \mathcal{D}^2 \hat{u} +\kappa^2 \hat{u}=G_1, \end{equation}

where $\kappa ^2= Ri\,Re\,\beta /D= Ri\,Re\,\lambda /2$. The left-hand side of (4.3) is the Bessel differential equation, which admits the Fourier–Bessel series as solutions. Substituting the boundary conditions, the mode shapes of $\hat {u}(r)$ should take the form

(4.4)\begin{equation} \hat{u}(r)=\frac{G_1}{\kappa^2}\left( 1- \frac{{\rm J}_0(\kappa r)}{{\rm J}_0(\kappa)} \right), \end{equation}

where ${\rm J}_m(r)$ is the $m$th Bessel function of the first kind. Enforcement of (4.2d) into (4.4) subsequently leads to an infinite number of discrete values of $\kappa \ (=\kappa _{c{,n}})$ satisfying ${\rm J}_2(\kappa _{c,n})=0$, at which (4.4) becomes a neutrally stable solution to (4.2). Here, $n$ indicates the $n$th zeros of ${\rm J}_2(r)$. In this case, $G_1$ in (4.4) becomes an arbitrary real constant, as (4.2d) is satisfied for any $G_1$. The values of $\kappa _{c{,n}}$ also yield the critical values of

(4.5)\begin{equation} Ri_{c{,n}}=2\kappa_{c{,n}}^2/(Re\,\lambda) \end{equation}

for the neutral stability of each mode. These values of $Ri_{c{,n}}$ calculated from $\kappa _{c,n}$ match perfectly with the bifurcation points computed numerically in the previous section, as shown in figure 1(a). Finally, it should be mentioned that (4.5) is equivalent to (3.14) in Pedley et al. (Reference Pedley, Hill and Kessler1988) and (31) in Bees & Hill (Reference Bees and Hill1999), where a continuous set of the critical values of the parameters equivalent to $Ri$ and $\kappa$ are obtained from linear stability analysis. However, in the present study, the introduction of a finite domain in the radial direction results in discrete values of $\kappa _{c{,n}}$ and $Ri_{c{,n}}$ with the corresponding eigenmode in the form of a cylindrical harmonic (i.e. Bessel functions) that satisfies the given boundary conditions.

4.2. Weakly nonlinear analysis

We further proceed to perform a weakly nonlinear analysis close to $Ri_{c{,n}}$. In the previous study by Bees & Hill (Reference Bees and Hill1999), where a uniform suspension is considered in an unbounded domain, it was assumed that the leading nonlinear term would appear at the third order due to translational invariance of the suspension in the horizontal direction. However, in the present study, such invariance is broken due to the pipe's cylindrical geometry. Hence, there is no reason that the leading nonlinear term would emerge at the third order. In this study, we therefore start by assuming $Ri-Ri_{c{,n}}=\epsilon {\rm \Delta} Ri$ (${\rm \Delta} Ri$ is the normalised distance from the bifurcation point) with a slow time scale $T=\epsilon t$. Given the perturbation form introduced in § 4.1, at $\textit{O}(\epsilon ^2)$, we get

(4.6a)\begin{equation} \partial_T u_1-{\rm \Delta} Ri \, n_1 = -\partial_t u_2-G_2 + \frac{1}{Re} \mathcal{D}^2 u_2 +Ri_{c{,n}} \, n_2, \end{equation}

where $G_2=(2/Re) u_2'(1)$ and

(4.6b)\begin{equation} \partial_T n_1+\beta \frac{1}{r}\partial_r ( r(\partial_r u_1) n_1 ) = -\partial_t n_2-\beta \mathcal{D}^2 u_2 + D \mathcal{D}^2n_2. \end{equation}

We also introduce amplitude $A(T)$ for the linear perturbation:

(4.7a,b)\begin{equation} u_1(r,t,T)=A(T)\hat{u}(r),\quad n_1(r,t,T)=A(T)\hat{n}(r), \end{equation}

where the linear instability mode is normalised to be $\hat {u}(0)=1$. Following the procedure of a weakly nonlinear analysis (e.g. Drazin Reference Drazin2002), we apply the solvability condition to (4.6) using the adjoint of (4.2). After some simplifications, we arrive at the amplitude equation

(4.8)\begin{equation} \left(C-\frac{F}{D\, Re} \right)\partial_T A- \frac{\lambda}{2} C {\rm \Delta} Ri\,A -\frac{\lambda E}{2 Re} A^2=0, \end{equation}

where $C$, $E$ and $F$ are defined in appendix A. Now, given the quadratic nonlinearity in (4.8), it is clear that all the bifurcations in figure 1(a) are transcritical. We have also computed the slopes of the non-trivial branches in figure 1(a), which are given by $-C\, Re/E$, when they cross each bifurcation point. As shown by the dot-dashed lines in figure 1(a), the computed slopes match with those of the numerically computed nonlinear solutions perfectly at the bifurcation points.

Finally, the neutrally stable eigenmodes (4.4) used for the weakly nonlinear analysis (red line, with appropriate sign) are compared with the numerical solutions (blue line) around each bifurcation point (at $U(0)={\pm }1$) in figure 1. As expected, there is excellent agreement between them.

5. Existence of steady solution

As previously noted by Hwang & Pedley (Reference Hwang and Pedley2014) and Fung et al. (Reference Fung, Bearon and Hwang2020), there are three parameters $Q$, $G$ and $Ri$ that control the bifurcation of (2.1). Here, $Q$ and $G$ are dependent on each other, providing only two degrees of freedom in total. Now, without loss of generality, we shall vary the three parameters $Q$, $G$ and $Ri$ in a controlled manner to explore the existence of the steady solutions reported in § 3. We first prescribe $G$ and continue the steady solutions to (2.1) by changing $Q$ – the related bifurcation diagrams for $G=0$ are also shown in figure 2 (dotted lines). Here, we note that the change of $Q$ for a given $G$ requires a change of $Ri$, given the relation of the three parameters. Figure 3(a) shows how $U(0)$, $Ri$ and $N(1)$ change with $Q$ for several prescribed $G$. It is found that all the solutions blow up at a certain respective threshold of $Q$ (say $Q_{c,G}(G)$): as $Q \rightarrow Q_{c,G}(G)$, $U(0)$ and $Ri$ blow up and $N(1)$ approaches zero (i.e. depletion of the cell number density at the wall).

Figure 3.

Continuations of the steady solution emerged from the first bifurcation point with $U(0)<0$. $(a)$ Plots of $U(0)$ (), $Ri$ () and $N(1)$ () for several $G$ on increasing $Q$ from $0$. $(b)$ The relation between $G$ and $Q$ for several fixed $Ri$. Inset: the maximum achievable flow rate $Q_{max}$ plotted for each $Ri$.

To further explain the existence of a threshold value of $Q$, we also perform continuation by changing $G$ for prescribed $Ri$. Figure 3(b) shows the continuation of the first two steady solution branches reported in figure 1(a) for fixed $Ri\in [200,1400]$. Now, it becomes apparent that there exists a maximum achievable $Q\ (\equiv Q_{max,Ri}(Ri))$ at each $Ri$. Plotting $Q_{max,Ri}(Ri)$ against $Ri$ (inset in figure 3b) shows that $Q_{max,Ri}$ reaches its maximum at $Q_{max} \approx 3.1$ and $Ri \approx 1174$. The maximum downward flow rate $Q_{max}\ (>0)$ is typically achieved with a downward pressure gradient ($G<0$) (i.e. $Q_{max}=\max (Q_{c,G}(G))$ in figure 3a). However, any further decrease of $G$ counter-intuitively decreases, rather than increases, the flow rate of the steady solution.

Lastly, we note that figure 3 is only for the first and second solution branches that emerge from $Q=0$ in figure 1(a). However, the same qualitative behaviours have been found from all the other solution branches: for example, the relation between $N(1)$ and $Q$, the blow-up of $U(0)$ and $Ri$ as $Q\rightarrow Q_{c,G}$, and the existence of $Q_{max}$.