Dispersion of a gyrotactic micro-organism suspension in a vertical pipe: the buoyancy–flow coupling effect

Abstract

Understanding the transport of micro-organisms in pipes is crucial to many fundamental problems, such as bioconvection and biodiesel production. In this work, we investigate the velocity profile and dispersion of a suspension of negatively buoyant, gyrotactic micro-organisms in a vertical pipe. With an imposed flow rate, the non-uniform radial cell concentration typical of gyrotaxis distorts the simple Poiseuille flow through inhomogeneous buoyancy, which in turn affects the cell concentration distribution. By solving the fundamental Smoluchowski equation and the Navier–Stokes equation simultaneously, we account for this bidirectional buoyancy–flow coupling effect. Asymptotic dispersion coefficients, namely, drift velocity and dispersivity, are further calculated with the obtained radial velocity and cell concentration profiles, which are assumed to be steady, symmetric and axially invariant. Using the gyrotactic micro-organism Chlamydomonas augustae as an example, detailed results are given to illustrate the effect of buoyancy–flow coupling. In downwelling flows, the buoyancy–flow coupling effect intensifies with the Richardson number $Ri$ quantifying the mean cell concentration, but is strongest at a moderate flow strength. The buoyancy–flow coupling effect significantly enhances the velocity and cell concentration in the central region, as well as the drift velocity and dispersivity. In contrast, the buoyancy–flow coupling effect is comparatively limited in upwelling flows, due to the dominant influence of the no-slip boundary condition imposed at the wall. Comparisons with predictions of existing approximate models are also presented.

1. Introduction

Motile micro-organisms are ubiquitous in various environments and exhibit fascinating phenomena. A classical phenomenon named ‘gyrotactic focusing’ was illustrated in a series of experiments by Kessler (1984, 1985a,b, 1986) using the suspension of Chlamydomonas augustae (mistakenly classified as Chlamydomonas nivalis Bees 2020, § 2.1). Gyrotactic focusing describes the centreline-accumulating behaviour of gyrotactic micro-organisms (a micro-organism whose mass centre is located at the far end of the buoyancy centre, e.g. C. augustae) in downwelling pipe flows (Denissenko & Lukaschuk 2007; Bees & Croze 2010; Bearon, Bees & Croze 2012; Croze, Bearon & Bees 2017; Jiang & Chen 2020) and channel flows (Croze et al. 2013). On the contrary, gyrotactic micro-organisms scatter towards the wall in upwelling flows (Kessler 1985a,b; Croze et al. 2013; Fung & Hwang 2020; Jiang & Chen 2020). In addition, because micro-organisms are typically heavier than water (Pedley & Kessler 1992), the focused micro-organisms in the central region exert an excessive negative buoyant force, thereby accelerating the local flow. This beam is subjected further to instability: when the background flow is weak and the average background cell concentration is high, axisymmetric blips with increased cell concentration and radius compared with the beam were observed (Kessler 1985a, 1986; Denissenko & Lukaschuk 2007).

Predicting the macrotransport process of micro-organisms in an external flow with a precise continuum model is crucial to relevant research and applications, such as understanding aforementioned instabilities (Hwang & Pedley 2014a,b; Saintillan 2014; Maretvadakethope, Keaveny & Hwang 2019; Bees 2020; Fung, Bearon & Hwang 2020; Fung & Hwang 2020) and optimising the efficiency of artificial cultivation (Chisti 2007). For gyrotactic micro-organisms, Kessler (1985a,b, 1986) introduced a primitive model. He proposed that the swimming directions of micro-organisms are deterministic and shear dependent, and micro-organisms are subject to isotropic diffusion. Kessler (1985a, 1986) subsequently derived an exponential profile of the steady radial cell concentration distribution in a downwelling pipe flow and analysed the instability. Pedley & Kessler (1990) later proposed a pioneering continuum model, which is denoted as the PK model in the remainder of this paper and also referred to as the FP model by some researchers. Compared with the model of Kessler (1986), the PK model incorporates biological random reorientation of micro-organisms and gives a more reasonable approximation for anisotropic diffusion in the position space. The PK model first finds the shear-dependent probability density distribution (p.d.f.) of the orientation by solving a Fokker–Planck equation in the orientation space and then uses a semi-empirical direction correlation time to calculate the diffusivity tensor.

Following the same routine to calculate the p.d.f. in the orientation space (and thus average swimming direction) with a specified shear rate, the generalised Taylor dispersion (GTD) theory (Frankel & Brenner 1989; Brenner & Edwards 1993) was introduced by Hill & Bees (2002) and Manela & Frankel (2003) to derive a more sophisticated but, at the same time, more sensible diffusivity tensor. It should be noted that both the PK model and the GTD model assume that micro-organisms relax in the orientation space with much smaller time and length scales relative to the corresponding scales of micro-organisms swimming in the position space, that is, quasi-steadiness in the orientation space. Bearon, Hazel & Thorn (2011) further extended the GTD model to flow fields with inhomogeneous shear, and a quasi-homogeneous assumption is further employed: the shear is slow varying, with time and length scales larger than those of relaxation in the orientation space.

Extensive efforts have been made to test the performance of the PK model and the GTD model. Croze et al. (2013) investigated the dispersion of gyrotactic micro-algae in the laminar channel flow with the PK model and the GTD model. Brownian dynamics simulation, in which an individual cell performs a random walk according to the stochastic differential equations (the equivalent of the Smoluchowski equation, i.e. the equation accurately describing the cell's rotation in the orientation space, movement in the position space and the effect of white noise), is employed as a benchmark; the dispersion in the turbulent regime was also simulated using Brownian dynamics simulation. They found that the three models agree well in weak flows, while the GTD model outperforms the PK model in strong flows. The difference originates from the GTD model providing a more rigorous diffusivity tensor than the PK model (Bearon et al. 2012; Fung et al. 2020): as the flow becomes stronger, the diffusivity tensor of the PK model approaches a constant scalar matrix asymptotically (in Cartesian coordinates), whereas the diffusivity tensor of the GTD model approaches a zero matrix. Many studies have indicated the superiority of the GTD model over the PK model, not only because of its better accuracy in investigating the dynamics of gyrotactic micro-organisms (Bearon et al. 2012; Croze et al. 2013; Fung et al. 2020; Fung, Bearon & Hwang 2022), but also due to its versatility in studying other kinds of micro-organisms, such as run–tumble chemotactic bacteria (Bearon 2003). Fung et al. (2020) recently analysed the bifurcation and stability of the suspension of C. augustae in downwelling pipe flows, employing both the GTD model and the PK model for comparison. With the GTD model, a trend consistent with the experimental observation of Kessler (1986) was found in the variation of stability as a function of the average cell concentration and flow rate; in contrast, the PK model predicted inconsistently.

Despite the foregoing merits, the GTD model has its inherent restriction. This restriction stems from simplifying the fundamental Smoluchowski equation to an advection–diffusion equation in position space. It is important to note that the GTD theory by Frankel & Brenner (1989) was not initially devised to facilitate a position–orientation separation to approximate the Smoluchowski equation. Instead, one should identify all bounded coordinates as the local space and all unbounded coordinates as the global space, and the obtained drift velocity and dispersivity only characterise the dispersion in the global space (Brenner & Edwards 1993; Jiang & Chen 2019, 2020). Rigorously, the GTD model extended by Hill & Bees (2002) as well as Manela & Frankel (2003) should be only applied to unbounded steady flows with homogeneous shear, and the obtained dispersion coefficients are valid only for long times. When the GTD model facilitating a position–orientation separation is applied to other situations, quasi-steadiness of the orientation dynamics of the micro-organisms, quasi-homogeneity of the flow shear and negligible boundary effect in the position–orientation space must be assumed to obtain an advection–diffusion equation in the position space. For example, Bearon et al. (2011) showed the discrepancies between the results of the GTD model and the Smoluchowski-equation-based simulation when the assumptions are not met. In another example, where the distribution of non-biased slender bacteria was considered in a channel flow, Bearon & Hazel (2015) showed how the GTD model predicted a uniform cell concentration profile across the width and failed to reproduce the experimentally observed shear-induced trapping (Rusconi, Guasto & Stocker 2014), whereas the results of the full numerical simulation of the Smoluchowski equation are in qualitative agreement with the experimental observations. Bearon & Hazel (2015) then showed a multi-time-scale analysis of the Smoluchowski equation where an effective drift velocity can be derived to explain the shear-induced trapping phenomenon. This multi-time-scale analysis is further formalised by Vennamneni, Nambiar & Subramanian (2020) into an advection–diffusion model, which can capture shear-induced migration of bacteria in a channel flow (Rusconi et al. 2014; Barry et al. 2015).

Although the model of Vennamneni et al. (2020) succeeded in capturing shear-induced trapping of micro-organisms, it only works when the motility is not biased, and therefore cannot be extended to case of cells with an orientation bias, such as gyrotactic micro-organisms. This limitation is overcome by the recent work of Fung et al. (2022), as they first devised an exact transformation of the Smoluchowski equation into an advection–diffusion equation in the position space, before applying a similar multi-time-scale method. In effect, Fung et al. (2022) has extended the work of Bearon & Hazel (2015) and Vennamneni et al. (2020) by including the effect of biased motility in the multi-time-scale asymptotics. This new model, however, still cannot elude the defects inevitably involved in the position–orientation separation.

As noted previously, the Smoluchowski equation can be referred to as a benchmark for other up-scaled models. Solving the full continuum Smoluchowski equation has therefore also attracted some researchers (Chen & Jiang 1999; Saintillan & Shelley 2008a,b; Ezhilan & Saintillan 2015; Jiang & Chen 2021). However, it is not always possible to resolve the full Smoluchowski equation in a three-dimensional flow field due to numerical difficulty. Only in some simplified circumstances, such as in unidirectional flows, can the dimension of the Smoluchowski equation be reduced. Recently, Jiang & Chen (2019) devised the basis functions and solved the Smoluchowski equation in the cross-section of a two-dimensional channel with a Galerkin method. The cross-sectional average drift and diffusivity were accurately derived by directly applying the standard GTD theory (Frankel & Brenner 1989; Brenner & Edwards 1993) without a position–orientation separation invoking assumptions of quasi-steadiness, quasi-homogeneity and negligible boundary effect. Jiang & Chen (2020) further investigated the dispersion of gyrotactic micro-organisms in pipe flows: a detailed comparison with the PK model and the GTD model was presented, explicitly demonstrating the inaccuracy of the PK model and the GTD model when the micro-organisms are highly motile or when the background flow is strong. Alternatively, as a compromised substitute of the Smoluchowski equation, Brownian dynamics simulations are mostly employed for computational convenience (Thorn & Bearon 2010; Bearon et al. 2011; Croze et al. 2013; Peng & Brady 2020). Although straightforward, if Brownian dynamics simulations are used, one has the problem to include the forces and stresses on the flow induced by the micro-organisms, and these effects are non-negligible in semi-dilute and dense suspensions.

The objective of the present study is to address the buoyancy–flow coupling effect with the transport of micro-organisms directly described by the Smoluchowski equation, hence extending the results in Jiang & Chen (2020). In the dilute regime, buoyancy–flow coupling is the dominant effect that negatively buoyant micro-organisms exert on the flow (Pedley & Kessler 1990), and is therefore crucial to the transport and stability of the suspension (Bees 2020). For the fundamental case within a vertical pipe, there are only a handful of works considering the buoyancy–flow coupling effect. Bees & Croze (2010) formulated a dispersion scheme resembling the classical Taylor–Aris dispersion theory (Taylor 1953; Aris 1956), using the asymptotic results from the PK model. Croze et al. (2017) presented an experimental study and directly compared the measured cross-sectional average drift velocity with predictions from the PK model and the GTD model. The bifurcation and stability in such a configuration were analysed by Fung et al. (2020) and Fung & Hwang (2020) for downwelling and upwelling flows using the PK model and the GTD model. Fung (2021, § 5.1) further performed a bifurcation analysis using an analytical solution of the simplified Smoluchowski equation in a vertical pipe considering the buoyancy–flow coupling effect. Although this solution is still coupled to the Navier–Stokes (N–S) equation (because it is related to the flow profile) and only applies to the steady state of spherical micro-organisms, it offers a more accurate basic state compared with those calculated by the PK model and the GTD model.

In this work, we consider a suspension flowing either downwards or upwards on average. The velocity and cell concentration profiles in the radial direction are assumed to be axially invariant and axisymmetric, as in Bees & Croze (2010) and Croze et al. (2017). The stability analysis of the solution, such as Fung et al. (2020) and Fung & Hwang (2020), is beyond the scope of this work. Due to the buoyancy–flow coupling effect, the flow velocity profile will deviate from a simple Poiseuille flow and must be solved simultaneously with the Smoluchowski equation. Compared with cases only considering the unidirectional effect of the flow on the cells (Jiang & Chen 2020), the bidirectional buoyancy–flow interaction couples the velocity and the cell transport, which is the major contribution of this work. In addition, the effect of buoyancy–flow coupling on the dispersion coefficients, namely the drift velocity and the dispersivity, will be discussed. At last, we shall present comparisons between the results of the Smoluchowski equation, the PK model and the GTD model, thus assessing their performances in the presence of the buoyancy–flow coupling effect. This assessment is difficult to realise in the Brownian dynamics simulations.

2. Formulation

2.1. Governing equations

As shown in figure 1, a cylindrical coordinate $(r^{\ast }, \psi, z^{\ast })$ identified by unit vectors $\boldsymbol {e}_r$, $\boldsymbol {e}_{\psi }$ and $\boldsymbol {e}_z$ characterises the dynamics in the position space and two angles $\theta$ and $\phi$ with the corresponding unit vectors $\boldsymbol {e}_{\theta }$ and $\boldsymbol {e}_{\phi }$ characterise the orientation of a micro-organism. Under the Boussinesq approximation, the N–S equation of the steady, unidirectional, axially invariant and axisymmetric flow velocity $\boldsymbol {u}^{\ast } = U^{\ast }(r^{\ast }) \boldsymbol {e}_z$ is

(2.1)\begin{equation} 0 ={-}\frac{1}{\rho^{{\ast}}} \frac{\mathrm{d} p^{{\ast}}}{\mathrm{d} z^{{\ast}}} + \nu^{{\ast}} \frac{1}{r^{{\ast}}} \frac{\mathrm{d}}{\mathrm{d} r^{{\ast}}} \left(r^{{\ast}} \frac{\mathrm{d} U^{{\ast}}}{\mathrm{d} r^{{\ast}}} \right) + n^{{\ast}} v^{{\ast}} \frac{\Delta \rho^{{\ast}}}{\rho^{{\ast}}} g^{{\ast}}.\end{equation}

On the right-hand-side of (2.1), the first term is the imposed pressure gradient ($\rho ^{\ast }$ is the density of water, and $\mathrm {d} p^{\ast }/\mathrm {d} z^{\ast }$ is the constant pressure gradient along the axial direction $\boldsymbol {e}_z$), the second term is the viscous stress ($\nu ^{\ast }$ is the kinematic viscosity of water) and the third term is the negative buoyancy of the cells ($n^{\ast }$ is the cell number density in the position space $(r^{\ast }, \psi, z^{\ast })$, $v^{\ast}$ is the volume of a cell, $\Delta \rho ^{\ast }$ is the density difference between the cell and water and $g^{\ast }$ is the gravitational acceleration).

Figure 1. Sketch of the problem. (a) In an exemplary downwelling flow, the cells accumulate in the central region and thus accelerate the local flow via their negative buoyancy. (b) The top view of the pipe. Unit vectors $\boldsymbol {e}_r$, $\boldsymbol {e}_{\psi }$ and $\boldsymbol {e}_z$ identify the cylindrical coordinate in the position space. (c) The spherical coordinate characterising the orientation of the micro-organism is featured with a polar angle $\theta$ between $\boldsymbol {p}$ and $\boldsymbol {e}_{\psi }$ and an azimuthal angle $\phi$ between $-\boldsymbol {e}_z$ and the projection of $\boldsymbol {p}$ on the $r$–$z$ plane.

In the N–S equation (2.1), the only effect of the micro-organisms on the flow is buoyancy. Here, we further consider the effect of active hydrodynamic stresslet (Saintillan 2018)

(2.2)\begin{equation} \boldsymbol{\varSigma}_p^{{\ast}} = n^{{\ast}} T^{{\ast}} \left(\int \boldsymbol{pp} \frac{f^{{\ast}}}{n^{{\ast}}} \, \mathrm{d}^2 \boldsymbol{p}- \frac{1}{3} {\boldsymbol{\mathsf{I}}}\right), \end{equation}

where $T^{\ast } = 10^{-10}\ \mathrm {g} \ \mathrm {cm}^2 \ \mathrm {s}^{-2}$ for C. augustae (Pedley 2010; Fung et al. 2020), $\boldsymbol{p}$ is the vector for the orientation of the micro-organism, and $\boldsymbol{\mathsf{I}}$ is the identity tensor. The $z$-component of $(1/\rho ^{\ast })\boldsymbol {\nabla }^{\ast }\boldsymbol {\cdot }\boldsymbol {\varSigma }_p^{\ast }$ is now added to the right-hand side of (2.1). A scaling analysis for the relative strength of the stresslet and buoyancy can be readily achieved by examining $T^{\ast }/(\nu ^{\ast } \Delta \rho ^{\ast } g^{\ast }R^{\ast })$. Substituting the parameters used in this work (shown later in table 1), we obtain a relative strength of $0.01$. Therefore, the stresslet can be safely neglected.

Table 1. Values and ranges of relevant parameters of C. augustae, most of which are taken from Pedley & Kessler (1990) and Hwang & Pedley (2014a,b).

Assuming that the total number of cells is conserved, the governing equation of the steady-state cell number density $f^{\ast }$ in the full position–orientation space $(r^{\ast }, \psi, z^{\ast }, \theta, \phi )$, namely the Smoluchowski equation, is

(2.3)\begin{equation} \boldsymbol{\nabla}^{{\ast}} \boldsymbol{\cdot} [(\boldsymbol{u}^{{\ast}} + V_s^{{\ast}} \boldsymbol{p} + V_d^{{\ast}} \boldsymbol{e}_z) f^{{\ast}} - D_t^{{\ast}} \boldsymbol{\nabla}^{{\ast}} f^{{\ast}}] + \boldsymbol{\nabla}_p \boldsymbol{\cdot} (\boldsymbol{\dot{p}}^{{\ast}} f^{{\ast}} - D_r^{{\ast}} \boldsymbol{\nabla}_p f^{{\ast}}) = 0. \end{equation}

Here, $\boldsymbol {\nabla }^{\ast }=\boldsymbol {e}_r (\partial /\partial r^{\ast }) + \boldsymbol {e}_{\psi } (1/r^{\ast })(\partial /\partial \psi ) + \boldsymbol {e}_z (\partial /\partial z^{\ast })$ and $\boldsymbol {\nabla }_p=\boldsymbol {e}_{\theta } (\partial /\partial \theta ) + \boldsymbol {e}_{\phi } (1/\sin \theta ) (\partial /\partial \phi )$ denote the gradient operators in the position space and in the orientation space, respectively; $V_s^{\ast }$ is the swimming velocity of the micro-organisms, with the orientation denoted by $\boldsymbol {p}$; $V_d^{\ast }$ is the settling velocity; $D_t^{\ast }$ and $D_r^{\ast }$ are respectively the intrinsic translational and orientational diffusivities (since we do not consider the hydrodynamic interaction between cells (Ishikawa & Pedley 2007), these two diffusivities are purely due to the physiology of the micro-organisms); $\boldsymbol {\dot {p}}^{\ast }$ is the time derivative of $\boldsymbol {p}$, and is related to the angular velocity of the micro-organism (denoted by $\boldsymbol {\varOmega }^{\ast }$) via

(2.4)\begin{align} \boldsymbol{\dot{p}}^{{\ast}} &= \boldsymbol{\varOmega}^{{\ast}} \times \boldsymbol{p},\nonumber\\ &=\left[\frac{1}{2} \boldsymbol{\omega}^{{\ast}} + \alpha_0 \boldsymbol{p} \times ({\boldsymbol{\mathsf{E}}}^{{\ast}} \boldsymbol{\cdot} \boldsymbol{p}) - \frac{1}{2B^{{\ast}}} (\boldsymbol{p} \times \boldsymbol{e}_z) \right] \times \boldsymbol{p},\nonumber\\ &= \frac{1}{2} \boldsymbol{\omega}^{{\ast}} \times \boldsymbol{p} + \alpha_0 \boldsymbol{p} \boldsymbol{\cdot} {\boldsymbol{\mathsf{E}}}^{{\ast}} \boldsymbol{\cdot} ({\boldsymbol{\mathsf{I}}}- \boldsymbol{pp}) + \frac{1}{2B^{{\ast}}}\boldsymbol{p} (\boldsymbol{p} \boldsymbol{\cdot} \boldsymbol{e}_z) - \frac{1}{2B^{{\ast}}}\boldsymbol{e}_z, \end{align}

with

(2.5)\begin{equation} \boldsymbol{\omega}^{{\ast}} = \boldsymbol{\nabla}^{{\ast}} \times \boldsymbol{u}^{{\ast}}, \end{equation}

being the vorticity of the flow, and

(2.6)\begin{equation} {\boldsymbol{\mathsf{E}}}^{{\ast}} = \tfrac{1}{2} [\boldsymbol{\nabla}^{{\ast}} \boldsymbol{u}^{{\ast}} + (\boldsymbol{\nabla}^{{\ast}} \boldsymbol{u}^{{\ast}})^{{T}}], \end{equation}

being the rate of strain; $B^{\ast }$ is the gravitactic reorientation time (smaller for a strong gyrotaxis); $\alpha _0=(\mathcal {AR}^2-1)/(\mathcal {AR}^2+1)$ characterises the shape of the cell ($\mathcal {AR} \geq 1$ is the aspect ratio of the spheroid). It is important to note that, since the micro-organisms considered in this work are spherical, i.e. $\alpha _0=0$, we have neglected the settling velocity component not parallel to $\boldsymbol {e}_z$ in (2.3), which is, however, non-zero for ellipsoids.

Integrating $f^{\ast }(r^{\ast }, \theta, \phi )$ over the orientation space yields the cell concentration in the position space, which is denoted by $n^{\ast }$

(2.7)\begin{equation} n^{{\ast}}(r^{{\ast}}) \triangleq \int_0^{\rm \pi} \int_0^{2{\rm \pi}} f^{{\ast}} \sin \theta \, \mathrm{d} \phi \, \mathrm{d} \theta. \end{equation}

From here on, we will eliminate the dependences of $f^{\ast }$ and $n^{\ast }$ on the position coordinates $\psi$ and $z^{\ast }$ using the symmetric and axially invariant assumption.

Dimensionless variables are introduced as

(2.8) \begin{equation} \left. \begin{array}{c@{}} \displaystyle t = \dfrac{t^{{\ast}} V_s^{{\ast}}}{R^{{\ast}}}, \quad r = \dfrac{r^{{\ast}}}{R^{{\ast}}}, \quad z = \dfrac{z^{{\ast}}}{R^{{\ast}}}, \quad D_t = \dfrac{D_t^{{\ast}}}{V_s^{{\ast}} R^{{\ast}}},\\ \displaystyle D_r=\dfrac{D_r^{{\ast}} R^{{\ast}}}{V_s^{{\ast}}},\quad \boldsymbol{u} = \dfrac{\boldsymbol{u}^{{\ast}}}{V_s^{{\ast}}}, \quad f=\dfrac{f^{{\ast}}}{N^{{\ast}}}, \quad n= \dfrac{n^{{\ast}}}{N^{{\ast}}},\quad p = \dfrac{p^{{\ast}}}{{\rho^{{\ast}} V_s^{{\ast}}}^2}, \\ \displaystyle \lambda = \dfrac{1}{2B^{{\ast}} D_r^{{\ast}}}, \quad Re = \dfrac{R^{{\ast}} V_s^{{\ast}}}{\nu^{{\ast}}}, \quad Ri = \dfrac{N^{{\ast}}v^{{\ast}} \Delta \rho^{{\ast}} g^{{\ast}}R^{{\ast}}}{\rho^{{\ast}}{V_s^{{\ast}}}^2}, \quad Se = \dfrac{V_d^{{\ast}}}{V_s^{{\ast}}}. \end{array}\right\} \end{equation}

Here, $D_t$ and $D_r$ are the dimensionless translational and rotational diffusivities, respectively; $\lambda$ is the dimensionless bias parameter; $f$ and $n$ are respectively the cell p.d.f. in the radius–orientation space $(r,\phi,\theta )$ and the dimensionless cell concentration in the radial direction ($r$) normalised by the mean cell concentration $N^{\ast }$, thus we have the normalisation condition

(2.9)\begin{equation} \int_0^1\int_0^{\rm \pi} \int_0^{2{\rm \pi}} 2 rf\sin \theta \, \mathrm{d} \phi \,\mathrm{d} \theta \,\mathrm{d}r = \int_0^1 2 r n \, \mathrm{d}r = 1; \end{equation}

$Re$ is the Reynolds number quantifying the viscous effect; $Ri$ is the Richardson number quantifying the mean cell concentration; and $Se$ is the dimensionless settling velocity.

The dimensionless form of (2.1) is

(2.10)\begin{equation} -\frac{\mathrm{d} p}{\mathrm{d} z} + \frac{1}{Re} \frac{1}{r} \frac{\mathrm{d}}{\mathrm{d} r} \left(r \frac{\mathrm{d} U}{\mathrm{d} r} \right)+ Ri \, n= 0, \end{equation}

and the dimensionless form of (2.3) becomes

(2.11)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} [(\boldsymbol{u}+\boldsymbol{p} + Se \boldsymbol{e}_z)f - D_t \boldsymbol{\nabla}f] + \boldsymbol{\nabla}_p \boldsymbol{\cdot} (\boldsymbol{\dot{p}}_r f - D_r \boldsymbol{\nabla}_p f)= 0. \end{equation}

Note that

(2.12)\begin{equation} \boldsymbol{\dot{p}}_r = \boldsymbol{\dot{p}} - \boldsymbol{\dot{p}}_c,\end{equation}

is the dimensionless time derivative of $\boldsymbol {p}$ expressed with basis vectors in the orientation space, namely $\boldsymbol {e}_{\theta }$ and $\boldsymbol {e}_{\phi }$, thus a term

(2.13)\begin{equation} \boldsymbol{\dot{p}}_c = \frac{\cos\theta}{r} \boldsymbol{e}_z \times \boldsymbol{p}, \end{equation}

representing the Coriolis effect must be subtracted from $\boldsymbol {\dot {p}}$ defined in the inertial coordinate. Specifically, $\boldsymbol {\dot {p}}_c$ emerges due to the swimming of micro-organism in the $\boldsymbol {e}_{\psi }$-direction and the dependence of spherical coordinate characterising the orientation of the micro-organism on the cylindrical coordinate fixed in the position space. The angular velocity of the rotating spherical coordinate is $\cos \theta /r \boldsymbol {e}_z$. The detailed expression of $\boldsymbol {\dot {p}}_r$ is

(2.14)\begin{equation} \boldsymbol{\dot{p}}_r = \dot{\theta} \boldsymbol{e}_{\theta} + \dot{\phi} \sin\theta \boldsymbol{e}_{\phi}, \end{equation}

where

(2.15)$$\begin{gather} \dot{\theta} = \cos\phi \cos\theta \left(D_r \lambda + \alpha_0 \frac{\mathrm{d} U}{\mathrm{d} r} \sin\phi \sin\theta \right) - \frac{\cos\theta \sin\phi}{r}, \end{gather}$$

(2.16)$$\begin{gather}\dot{\phi} ={-}\frac{1}{2} \frac{\mathrm{d} U}{\mathrm{d} r} + \frac{1}{2} \alpha_0 \frac{\mathrm{d} U}{\mathrm{d} r} \cos 2\phi - D_r \lambda \csc\theta \sin\phi - \frac{\cos\phi \cos\theta \cot\theta}{r}. \end{gather}$$

The steady-state Smoluchowski equation (2.11) for the cell p.d.f. $f(r, \theta, \phi )$ is expressed as

(2.17)\begin{align} \mathcal{L} f &\triangleq -\sin\theta \sin\phi \frac{1}{r} \frac{\partial (r f)}{\partial r} - D_t \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial f}{\partial r}\right) +\frac{1}{\sin\theta} \frac{\partial (\dot{\theta} \sin\theta f)}{\partial \theta} + \frac{\partial (\dot{\phi}f)}{\partial \phi}\nonumber\\ &\quad - D_r \left[ \frac{1}{\sin\theta} \frac{\partial}{\partial \theta} \left( \sin\theta \frac{\partial f}{\partial \theta} \right) + \frac{1}{\sin^2\theta} \frac{\partial^2 f}{\partial \phi^2} \right] = 0, \end{align}

where $\mathcal {L}$ denotes the operator characterising the transport in the radius–orientation space $(r, \theta, \phi )$.

2.2. Solutions for velocity and cell p.d.f.

The total flow rate is calculated by

(2.18)\begin{equation} Q = \int_0^1 2 {\rm \pi}rU(r) \, \mathrm{d}r ={\rm \pi} U_m, \end{equation}

where $U_m$ denotes the mean flow velocity. We shall prescribe the mean flow velocity $U_m$, then there are two variables $U(r)$ and $f(r, \theta, \phi )$ along with (2.10) and (2.17). Theoretically, the coupled equations can be solved with appropriate constraints and boundary conditions.

The boundary conditions for $U(r)$ are a no-slip condition at the wall and a boundedness condition at the centreline

(2.19)$$\begin{gather} U(1) = 0, \end{gather}$$

(2.20)$$\begin{gather}U(0) \neq \infty. \end{gather}$$

For $f(r, \theta, \phi )$, following Jiang & Chen (2020), we use specular reflection boundary conditions to ensure the no-flux requirement of cell concentration at the wall

(2.21)$$\begin{gather} f(r, \theta, \phi, t)|_{r=1} = f(r, \theta, -\phi, t)|_{r=1}, \end{gather}$$

(2.22)$$\begin{gather}\left. \frac{\partial f(r, \theta, \phi, t)}{\partial r}\right|_{r=1} =\left. -\frac{\partial f(r, \theta, -\phi, t)}{\partial r}\right|_{r=1}. \end{gather}$$

We note that the reflection boundary conditions may not be appropriate for micro-organisms featuring wall adhesion, for example, non-biased bacteria accumulate at walls (Berke et al. 2008; Li & Tang 2009). However, we suppose that the reflection boundary conditions are suitable for a vertical pipe configuration, because most of the suspended gyrotactic micro-organisms accumulate in the central region in a downwelling flow and the high shear rate at the wall dominates the reorientation of cells in an upwelling flow. Furthermore, as noted previously, the PK model and the GTD model are unable to incorporate boundary conditions in the position–orientation space, and no-flux boundary conditions in the position space are imposed; therefore, using the reflection boundary conditions does not affect the comparisons between the PK model, the GTD model and the Smoluchowski equation (Bearon et al. 2011; Croze et al. 2013; Bearon & Hazel 2015; Jiang & Chen 2020).

The boundedness condition for $\theta$ is

(2.23a,b)\begin{equation} f|_{\theta=0} \neq \infty, \quad f|_{\theta={\rm \pi}} \neq \infty, \end{equation}

and the periodic boundary conditions for $\phi$ are

(2.24)$$\begin{gather} f_{\phi=2 {\rm \pi}} = f|_{\phi=0}, \end{gather}$$

(2.25)$$\begin{gather}\frac{\partial f}{\partial \phi}\bigg|_{\phi=2 {\rm \pi}} = \frac{\partial f}{\partial \phi}\bigg|_{\phi=0}. \end{gather}$$

In this work, we shall use a Galerkin method for the Smoluchowski equation and an eigenfunction expansion method for the N–S equation. To resolve the buoyancy–flow coupling, an iterative method is employed (Hwang & Pedley 2014b). The details of the numerical method are given in Appendix A.

2.3. Solutions for dispersion coefficients

The drift velocity $U_d$ and the dispersivity $D_{{T}}$ constitute two key parameters in the drift–dispersion equation for the cross-sectional average cell concentration at the asymptotically long time. However, care must be taken in understanding these two coefficients in the current situation: they only describe the dispersion of a small number of labelled cells in an existing background flow formed by the imposed flow and a large number of unlabelled cells. Specifically, compared with the unlabelled cells, the distortion on the flow induced by the labelled cells is negligible. This assumption decouples the dispersion of labelled cells from the background buoyancy–flow coupling problem and ensures the background flow is always steady (Bees & Croze 2010; Croze et al. 2017), which are required for the application of Taylor–Aris theory (Taylor 1953; Aris 1956; Guan et al. 2022; Jiang et al. 2022) and the GTD theory (Frankel & Brenner 1989; Brenner & Edwards 1993; Jiang & Chen 2019, 2020).

With the solved velocity profile $U(r)$ and cell p.d.f. profile $f(r, \theta, \phi )$, one can apply the standard GTD theory (Frankel & Brenner 1989; Jiang & Chen 2019, 2020) to evaluate the cross-sectional average drift velocity $U_d$ and dispersivity $D_{{T}}$. Let $\langle ({\cdot }) \rangle$ denote the integration $2 \int _0^1 \int _0^{{\rm \pi} } \int _0^{2{\rm \pi} } ({\cdot }) r \sin \theta \, \mathrm {d} \phi \,\mathrm {d} \theta \,\mathrm {d} r$, the cross-sectional average drift velocity $U_d$ and dispersivity $D_{{T}}$ are calculated by Jiang & Chen (2020, (3.34) and (3.35))

(2.26)$$\begin{gather} U_d = \langle V_z(r, \theta, \phi) f(r, \theta, \phi) \rangle, \end{gather}$$

(2.27)$$\begin{gather}D_{{T}} = D_t + \langle V_z(r, \theta, \phi) b(r, \theta, \phi) \rangle, \end{gather}$$

where

(2.28)\begin{equation} V_z(r, \theta, \phi) = U(r) + p_z(\theta, \phi) + Se = U(r) - \sin\theta \cos\phi + Se, \end{equation}

is the local effective velocity, and $b(r, \theta, \phi )$ is a scalar field constituting the solution of the equation (Frankel & Brenner 1989; Brenner & Edwards 1993; Haugerud, Linga & Flekkøy 2022)

(2.29)\begin{equation} \mathcal{L} b = f(V_z - U_d).\end{equation}

An additional normalisation condition for $b(r, \theta, \phi )$ is imposed to ensure its uniqueness

(2.30)\begin{equation} \langle b(r, \theta, \phi) \rangle =0,\end{equation}

note that this normalisation condition does not affect the dispersivity. Then $b(r, \theta, \phi )$ is also solved by a Galerkin method (see Appendix A for details).

3. Results and discussion

Following Jiang & Chen (2020), among others (Pedley & Kessler 1992; Bees & Hill 1998; Hwang & Pedley 2014a,b; Zeng & Pedley 2018), parameters of the gyrotactic species C. augustae are used for calculation to illustrate the buoyancy–flow coupling effect, as listed in table 1.

It is noteworthy that, at high flow rate, the Richardson number $Ri$ must be below a threshold $Ri_c$ to avoid blow-up of the focused beam, i.e. to ensure the existence of a steady solution. We adopt a formula from Fung et al. (2020, § 3.3) to estimate this threshold, that is, $Ri_c=16/(\lambda Re)$. This formula is deduced using the analytical solution of the simplified Smoluchowski equation (Fung 2021, § 5.1), and can be also obtained from a simplified model resembling the GTD model. Substituting the values of $\lambda$ and $Re$ in our case, we obtain $Ri_c=122.85$. If $Ri$ is larger than $Ri_c$ at a high flow rate, a steady solution may not exist and the numerical method in this work becomes invalid. It is acknowledged that this threshold numerically estimated by the full Smoluchowski equation may slightly differ from the one given by the formula which assumes that all cells are concentrated in a narrow core around the axis; however, since the bifurcation structure is not the focus of this work, we do not try to find out the exact $Ri_c$ and an upper limit of $Ri=100$ is tested to guarantee a steady solution.

3.1. Velocity and cell concentration profiles

We first report the results in downwelling flows, as shown in figure 2 (the profiles of the velocity deviation above the mean) and figure 3 (the cell concentration profiles) for three Richardson numbers ($Ri=0$, $36.6$ and $73.1$) and three mean flow velocities ($U_m=1$, $5$ and $10$). The velocity deviation above the mean is defined as

(3.1)\begin{equation} \chi(r) = \frac{U(r)}{U_m}-1. \end{equation}

Figure 2. The velocity deviation profiles $\chi (r)$ in the downwelling flows with $Ri=0$, $36.6$ and $73.1$. The cell transport is calculated using the Smoluchowski equation; (a) $U_m=1$, (b) $U_m=5$, (c) $U_m=10$.

Figure 3. The cell concentration profiles $n(r)$ in the downwelling flows with $Ri=0$, $36.6$ and $73.1$. The cell transport is calculated using the Smoluchowski equation; (a) $U_m=1$, (b) $U_m=5$, (c) $U_m=10$.

As seen in figure 2, when $Ri=0$, i.e. the buoyancy–flow coupling effect vanishes and the velocity profile is completely decoupled from the cell concentration profile, the flow is exactly a Poiseuille flow characterised by $\chi (0)=1$ and $\chi (1)=-1$ with $r=0$ and $1$ corresponding to the centreline and wall, respectively. As $Ri$ increases, we see that the flow profile deviates from the Poiseuille flow, as evidenced by an increased velocity at the centreline ($\chi (0)>1$). In addition, we observe that $\chi (0)$ varies non-monotonically with $U_m$ when $Ri$ is set to positive (see also figure 5c), indicating a competition between the buoyancy–flow coupling effect and the background flow. We shall investigate this non-monotonicity quantitatively later. The cell concentration profiles typical of gyrotactic focusing at the central region in the downwelling flows are depicted in figure 3. The increases of both $U_m$ and $Ri$ have a profound enhancement on the gyrotactic focusing, as reflected by the increased cell concentration at the centreline.

Now we turn to the upwelling flows, the velocity deviation profiles with $U_m=-10$ are shown in figure 4(a). We find no discernible difference between the profiles obtained with different $Ri$ at this selected $U_m$, which is relatively large, i.e. profiles almost identical to an upwelling Poiseuille flow are predicted. This observation is the key result in the upwelling cases, and is also in deep contrast to the downwelling cases where the buoyancy–flow coupling effect has a significant influence on the velocity deviation profile even when $|U_m|=10$. In the upwelling flows, the gyrotactic cells scatter towards the wall where they were expected to distort the flow via their negative buoyancy. Nonetheless, near the wall, the no-slip condition (2.19) dominates the buoyancy–flow coupling effect, and thus the velocity deviation profiles do not vary as much as they do in the downwelling cases when $Ri$ is increased. Although we see in figure 5(c) that the buoyancy–flow coupling effect still accelerates the peripheral flow (in the $\boldsymbol {e}_x$ direction) and consequently leads to an increased velocity deviation at the centreline ($\chi (0)>1$), this effect vanishes faster than that in the downwelling flows as $|U_m|$ increases. We present the cell concentration profiles in the upwelling flows with $U_m=-10$ in figure 4(b), and the almost identical cell concentration profiles at different $Ri$ can be understood by the almost identical velocity deviation profiles. Meanwhile, since micro-organisms are mostly located in the near-wall region where the buoyancy–flow coupling effect is insignificant, we do not see a notable difference between the cell concentration at the wall $n(1)$ at different $Ri$ as shown in figure 5(b).

Figure 4. (a) The velocity deviation profiles $\chi (r)$. (b) The cell concentration profiles $n(r)$ in the upwelling flows. The cell transport is calculated using the Smoluchowski equation. Note that the velocity deviation profiles $\chi (r)$ at all sampled mean flow velocity ($U_m=-1$, $-5$ and $-10$) are very close, and therefore only profiles with $U_m=-10$ are shown here.

Figure 5. (a) Variations of the cell concentration at the centreline $n(0)$ as a function of the mean velocity $U_m$ in the downwelling flows. (b) Variations of the cell concentration at the wall $n(1)$ as a function of the mean velocity $U_m$ in the upwelling flows. (c) Variations of the velocity deviation at the centreline $\chi (0)$ as a function of the mean velocity $U_m$. Here, ‘S’, ‘G’ and ‘P’ denote results obtained from the Smoluchowski equation, the GTD model and the PK model, respectively.

As noted previously, there exists a competition between the buoyancy–flow coupling effect and the background flow, as evidenced by the non-monotonic variation of $\chi (0)$ with increasing $U_m$ (see e.g. figures 3 and 5a): $\chi (0)$ first experiences an increase and then a decrease to unity as $U_m$ increases when $Ri>0$. An intuitive impression on the strength of the buoyancy–flow coupling effect in the $U_m$–$Ri$ plane is given in figure 6, where the contours of the velocity deviation at the centreline $\chi (0)$ and the cell concentration at the centreline $n(0)$ are plotted. We observe that, while $\chi (0)$ increases monotonically with $Ri$ regardless of the value of $U_m$, $\chi (0)$ is maximally enhanced at a moderate $U_m$ (approximately equals $2$). In contrast, $n(0)$ increases monotonically with both $U_m$ and $Ri$.

Figure 6. Contours of (a) the velocity deviation at the centreline $\chi (0)$ and (b) the cell concentration at the centreline $n(0)$ in the $U_m$–$Ri$ plane.

3.2. Drift velocity and dispersivity

Figure 7 illustrates the variations of the drift velocity $U_d$ and the dispersivity $D_{{T}}$ as functions of the mean flow speed $U_m$ for three Richardson numbers ($Ri=0$, $36.6$ and $73.1$). In the downwelling flows, $U_d$ increases monotonically with $U_m$ and $Ri$. This trend can be well understood by the accelerated flow velocity and concentrated cell concentration in the central region, both of which contribute to a larger $U_d$. For $D_{{T}}$, we observe a boost with increasing $Ri$, which can be probably explained by the enhanced shear in the central region benefiting dispersion; however, since $D_{{T}}$ is a quantity related to the complex scalar field $b(r, \theta, \phi )$, an intuitive understanding of its variation with the flow strength and the cell concentration profile is potentially difficult. For example, as shown in figure 7(a), in the downwelling flows, $D_{{T}}$ exhibits a non-monotonic variation with $U_m$: $D_{{T}}$ first increases and then decreases as $U_m$ increases, which is also reported in non-coupling cases (Croze et al. 2013; Jiang & Chen 2020) and coupling cases (Croze et al. 2017). This non-monotonic variation is due to the effect of the enhanced gyrotactic focusing on the dispersion. We shall consider two extreme cases to understand this variation: first, there is no flow, and thus dispersivity $D_{{T}}$ is solely contributed by swimming, which is very small due to biased swimming induced by gyrotaxis; second, there is an extremely strong flow, and thus nearly all the cells are located in the centre, meaning that differential advection again contributes nothing to dispersivity; between these two extremes, the dispersivity is increased by both the swimming-induced dispersion and differential-advection-induced dispersion.

Figure 7. Variations of drift velocity $U_d$ and dispersivity $D_{{T}}$ as functions of the mean velocity $U_m$ (a,b) in the downwelling flows and (c,d) in the upwelling flows. Here, ‘S’, ‘G’ and ‘P’ denote results obtained from the Smoluchowski equation, the GTD model and the PK model, respectively.

In the upwelling flows, $U_d$ is negative because both the advective velocity of the background flow and the average swimming velocity of the gyrotactic micro-organisms are directed upwards. In contrast to the downwelling cases, $U_d$ finally saturates in the upwelling flows. A sensible explanation for this phenomenon is that, in the limit of extremely strong upwelling flow, all micro-organisms scatter towards the wall, thus $U_d$ is independent of $U_m$ and is solely determined by the near-wall orientation distribution of the micro-organisms (which determines the swimming-induced drift). The non-monotonic variation of $D_{{T}}$ with $U_m$ in the upwelling flows is analogous to that observed in the downwelling flows, and is also due to an initial increment and a following reduction of the differential-advection-induced dispersion as the flow becomes strong. Although the buoyancy–flow coupling effect is relatively weak in the upwelling flows because of the dominance of the no-slip wall, we are still able to observe slight differences of $U_d$ and $D_{{T}}$ between cases with different $Ri$ at a moderate range of $U_m$ (around $2$). It is seen that $U_d$ is increased with a larger $Ri$, this is due to the fact that the micro-organisms now accelerate the flow near the wall in the $\boldsymbol {e}_x$ direction, while the background flow is in the $-\boldsymbol {e}_x$ direction; this counteracting effect finally leads to the increase of $U_d$ with $Ri$. The decrease of $D_{{T}}$ with $Ri$ can be probably also understood by the modification on the flow profile: the shear in the near-wall region where the micro-organisms accumulate is decreased, and thus the contribution of differential advection to the dispersion is weakened.

3.3. Comparisons with PK model and GTD model

For dispersion of gyrotactic micro-organisms in a pipe, we can employ a simplified method for the calculation using the PK model and the GTD model, see e.g. Bearon et al. (2012, Appendix E). Now using the parameters of C. augustae and the calculation method presented in Appendix B, we shall compare the results obtained from the Smoluchowski equation with those from the PK model and the GTD model (the comparisons are already included in figure 5 and figure 7). It is also noteworthy that Croze et al. (2017) used a different species called Dunaliella salina in their experiment. We note that $\lambda =0.21$ for D. salina, which means that their gyrotactic focusing will be weak, and thus the effect of buoyancy–flow coupling will be less significant. In Appendix C, we present the comparisons between the results of D. salina calculated using the PK model, the GTD model, the Smoluchowski equation and the results measured in the experiments.

We first discuss the comparisons in the downwelling flows. As shown in figure 5(a), when $Ri=0$, the Smoluchowski equation and the GTD model predict almost the same $n(0)$ profiles, typical of gyrotactic focusing, while $n(0)$ predicted by the PK model is very close to unity. As $Ri$ increases, we see that the flow deviates from a Poiseuille flow for all models, as evidenced by an increased velocity at the centreline $\chi (0)>1$ (figure 5c). However, the buoyancy–flow coupling effect is less profound predicted by the PK model, while the GTD model slightly overestimates the effect at a moderate $U_m$ and underestimates it at a large $U_m$. The unsatisfactory performance of the PK model in high-shear flows has long been recognised (Croze et al. 2013, 2017; Bees 2020; Fung et al. 2020; Jiang & Chen 2020). This weakness is because the diffusion tensor calculated by the PK model cannot reflect the shear-induced dispersion in the position space (Croze et al. 2013). In a Poiseuille flow, if the flow rate continues increasing, the cell concentration profile predicted by the PK model will tend to be uniform, e.g. Bearon et al. (2012, figure 2) and Croze et al. (2017, figure 1a), as a result of the mean swimming velocity approaching to zero and the diffusivity tensor approaching to a constant scalar matrix (in Cartesian coordinate). In other words, gyrotactic micro-organisms will behave in a similar manner to passive tracers if one applies the PK model in a high-shear flow. For the buoyancy–flow coupling effect, although the velocity deviation at the centreline $\chi (0)$ does reflect the buoyancy–flow coupling effect to some extent, we do not see a noteworthy variation of $n(0)$ as $Ri$ increases with the PK model.

The velocity deviation at the centreline $\chi (0)$ predicted by the Smoluchowski equation and the GTD model are close (figure 5c). Nonetheless, when it turns to the cell concentration at the centreline $n(0)$, an apparent difference is seen when $Ri > 0$, as is seen in figure 5(a). The GTD model underestimates the gyrotactic focusing in the presence of the buoyancy–flow coupling effect. This result is due to these models sharing the same governing equation for the flow but fundamentally different governing equations for the cell transport (Bees 2020). In the analytical solution of the Smoluchowski equation (Fung 2021, § 5.1) and in models assuming a position–orientation separation, the radial cell concentration has an exponential scaling with the flow profile (see e.g. Fung 2021, (5.1)) or the ratio between the radial components of the local drift velocity and the diffusivity tensor (see e.g. (B6)), which also depends on the flow profile. Thus, any perturbation on the flow profile will result in a remarkable change in the cell concentration, causing the cell concentration to be more sensitive to the models.

Now we turn to cases where the flow is upwards, the comparisons of the velocity deviation at the centreline $\chi (0)$ and the cell concentration at the wall $n(1)$ with varying $U_m$ are shown in figures 5(c) and 5(b). We find that all models predict almost identical $\chi (0)$. Furthermore, as the flow becomes strong, $\chi (0)$ values obtained with all models approach unity. This coincidence is again in deep contrast with the downwelling cases where the cell transport model has a significant influence on $\chi (0)$, and can be also attributed to the no-slip boundary condition imposed at the wall. For $n(1)$, we see relatively good agreement between the predictions of the Smoluchowski equation and the GTD model at weak and moderate flow strengths, deviation is, however, observed in strong flow. In contrast, the PK model yields a much more uniform concentration profile and a decreasing trend of $n(1)$ as the flow becomes strong when $U_m$ is smaller than $-2$.

Figure 7(a) compares the variations of the drift velocity $U_d$ as functions of the mean flow speed $U_m$ in the downwelling flows between different models. We see that $U_d$ predicted by the GTD model and the Smoluchowski equation are nearly twice as large as the one predicted by the PK model with a large $U_m$. This is because both the GTD model and the Smoluchowski equation can effectively capture the gyrotactic focusing phenomenon in a strong flow. Specifically, cells accumulate at the centreline where the local advective speed is twice the mean flow velocity in non-coupling cases. In contrast, $U_d$ predicted by the PK model gradually approaches $U_m$, as a consequence of the nearly uniform radial cell concentration profile resembling that of solute tracers. For the effect of buoyancy–flow coupling, we see that the GTD model and the Smoluchowski equation still give almost identical $U_d$ (figure 7a). Although we have seen that the cell concentration at the centreline $n(0)$ of the Smoluchowski equation is significantly larger than that of the GTD model (figure 5a), this concentration enhancement in the central region does not lead to a discernible increase of $U_d$, since the area of the narrow central core scales with $r^2$. In addition, the curves of the PK model at all sampled $Ri$ nearly overlap except for a very small $U_m$.

The variation of dispersivity $D_{{T}}$ as a function of the mean flow speed $U_m$ in the downwelling flows between different models is presented in figure 7(b). Different from the results of the Smoluchowski equation, the GTD model predicts a saturated instead of a decreasing $D_{{T}}$. In contrast, although the PK model captures the variation of $D_{{T}}$ with increasing $Ri$, it fails to reflect neither the decrease nor the saturation of $D_{{T}}$ when the flow becomes strong. In summary, for the cases in the downwelling flows, the GTD model outperforms the PK model, but still overestimates $D_{{T}}$.

The dispersion predictions in the upwelling flows by different models are plotted in figures 7(c) and 7(d). First, although $U_d$ of the GTD model finally becomes nearly a constant with an initial drop as $|U_m|$ increases, an extra turning point is observed compared with the curve predicted by the Smoluchowski equation. In comparison, the curves of the PK model still behave very differently from the other two models: a tracer-like dispersion behaviour is predicted again as $|U_m|$ increases. For the dispersivity $D_{{T}}$, although the GTD model also predicts a non-monotonic trend as the flow becomes strong, the agreement with the results of the Smoluchowski equation remains qualitative.

4. Concluding remarks

In this paper, by solving the N–S equation and the Smoluchowski equation simultaneously, we investigate the effect of buoyancy–flow coupling on the velocity profile, cell concentration profile and Taylor dispersion of a suspension of gyrotactic micro-organisms in a vertical pipe. The buoyancy–flow coupling effect is more significant when the average flow is directed downwards, as a consequence of gyrotactic focusing of negatively buoyant micro-organisms accelerating the flow in the central region. In contrast, the buoyancy–flow coupling effect has a minor influence in the upwelling flows because the negative buoyancy of the wall-accumulated micro-organisms is inferior to the effect of the no-slip condition imposed at the wall. In the downwelling flows, the buoyancy–flow coupling effect leads to enhanced flow velocity and gyrotactic focusing in the central region, as well as enhanced drift velocity and dispersivity. The strength of the buoyancy–flow coupling effect can be quantified by the Richardson number $Ri$; however, its modification to the flow, which we use the velocity deviation at the centreline $\chi (0)$ to quantify, is largest at a moderate mean flow velocity $U_m$. This phenomenon is understood by considering that a small $U_m$ does not arouse obvious gyrotactic focusing and thus the cell concentration is too uniform to result in acceleration of the flow in the central region, while a large $U_m$ due to the imposed pressure gradient overwhelms the buoyancy–flow coupling effect.

Another contribution of this paper is giving a benchmark test for the PK model and the GTD model in the presence of the buoyancy–flow coupling effect. It has been shown in previous studies that the PK model performs unsatisfactorily in strong shear flows (Croze et al. 2013, 2017; Jiang & Chen 2020). Here, we further confirm its defect in the presence of the buoyancy–flow coupling effect. In comparison, although the GTD model captures some trends in the variations, its predictions of the dispersivity still deviate a lot from the exact results of the Smoluchowski equation.

One may inevitably question the extent to which the Smoluchowski equation can explain the pioneering experimental measurements in Croze et al. (2017). We show the comparisons in Appendix C; however, no obvious improvement on the GTD model is seen, since $\lambda =0.21$ for D. salina and thus assumptions of the GTD model basically hold. Although somewhat disappointed, we still advocate using the Smoluchowski equation in relevant problems, at least as a benchmark test, because we have shown the inaccuracy of the GTD model when the assumptions associated with orientation–position separation are not met. The fact that the experimental results do not agree satisfactorily with the GTD model and the Smoluchowski equation is open to many possible explanations, such as heterogeneity of micro-organisms (Bees, Hill & Pedley 1998), cell–cell interactions in the central core (Fung et al. 2020) and the microscopic kinematic model of the cell (Omori et al. 2022; Walker et al. 2022; Zeng, Jiang & Pedley 2022), which deserve in-depth research.

Another important issue not covered by the current study is the possible bifurcation and stability of the steady-state solutions. Since we have limited the Richardson number $Ri$ to less than $100$, we can safely evade the regimes with bifurcation and instability. It is of great interest to extend the current work to linear hydrodynamic stability analysis of the solutions; however, great difficulty arises in converting the linearised equation to an eigenvalue problem due to the non-parallel perturbed flow field and complex axial transport of cells. Fung (2021, § 5.1) has given a preliminary analysis and this problem should be explored in future work.