Bioconvective Linear Stability of Gravitactic Microorganisms

Ildebrando Pérez-Reyes and Luis Antonio Dávalos-Orozco

## Abstract

Interesting results on the linear bioconvective instability of a suspension of gravitactic microorganisms have been calculated. The hydrodynamic stability is characterized by dimensionless parameters such as the bioconvection Rayleigh number R, the gyrotaxis number G, the motility of microorganisms d, and the wavenumber k of the perturbation. Analytical and numerical solutions are calculated. The analytical one is an asymptotic solution for small wavenumbers (and for any motility number) which agrees very well with the numerical solutions. Two numerical methods are used for the sake of comparison. They are a shooting method and a Galerkin method. Marginal curves of R against k for fixed values of d and G are presented along with curves corresponding to the variation of the critical values of Rc and kc. Moreover, those critical values are compared with the experimental data reported in the literature, where the gyrotactic algae Chlamydomonas nivalis is the suspended microorganism. It is shown that the agreement between the present theoretical results and the experiments is very good.

Keywords: bioconvection, hydrodynamic stability, Galerkin method

#### 1. Introduction

Since many years ago, efforts in the experimental and theoretical investigation of the bioconvection phenomenon have been made. These efforts, which lead to the understanding of bioconvective instability, have produced novel and interesting applications. For example, Noever and Matsos [1] proposed a biosensor for monitoring the heavy metal Cadmium based in bioconvective patterns as redundant technique for analysis, a number of researchers [2–6] have been working on the control of bioconvection by applying electrical fields (as in galvanotaxis) to use it as a live micromechanical system to handle small objects immersing in suspensions, Itoh et al. [7, 8] use some ideas of bioconvection in a study for the motion control of microorganism groups like Euglena gracilis to manipulate objects by using its phototactic orientation (as in phototaxis), and more recently possibly bioconvection seeded the investigation of Kim et al. [9, 10] for using a feedback control strategy to manipulate the motions of Tetrahymena pyriformis as a microbiorobot, among others. Perhaps, further applications on biomimetics [11–13] at the nano- and microscale could be driven by this contribution.

The term bioconvection was first coined by Platt [14] as the spontaneous pattern formation in suspensions of swimming microorganisms. This phenomenon has some similarity with Rayleigh-Benard convection but originates solely from

diffusion and the swimming of the organisms. Reviews about this topic have been published by Pedley and Kessler [15] and Hill and Pedley [16]. Ideas and theories on cellular motility can be found in the book of Murase [17], and the effect of gravity on the behavior of microorganisms is widely explained in the book of Hader et al. [18]. In 1975, Childress et al. [19] presented a model for bioconvection of purely gravitactic microorganisms and their results of a linear theoretical study, and later Harashima et al. [20] studied the nonlinear equations of this model. According to the model of Childress et al. [19], the critical wavenumber at the onset of the instability is always zero. In ordinary particles and colloidal suspensions, the internal degrees of freedom like the internal rotation or spin are important under some geometrical and physical conditions [21, 22]. The case of a suspension of microorganisms is not an exception. For this case, Pedley et al. [23] proposed a gyrotactic model for a suspension of infinite depth. Their model includes the displacement of the gravity from the geometric center in the organisms along their axis of symmetry. Hill et al. [24] performed an analysis of the linear instability of a suspension of gyrotactic microorganisms of finite depth using the model of Pedley et al. [23]. Hill et al. [24] found finite wavenumbers at the onset of the instability, but agreement with the experiment was not good. Later, Pedley and Kessler [25] reported a model for suspensions of gyrotactic microorganisms where account was taken of randomness in the swimming direction of the cells. In a study of the linear instability of the system based on the model of Pedley and Kessler [25], Bees and Hill [26] found disagreement between their theoretical results and the experimental data reported by Bees and Hill [27]. Several experimental investigations of bioconvection have been reported by Loeffer and Mefferd [28] and Fornshell [29], by Kessler [30] and Bees and Hill [27] who take into account the gyrotaxis, by Dombrowski et al. [31] and Tuval et al. [32] who take into account the oxitaxis, and more recently by Akiyama et al. [33] who observed a pattern alteration response characterized by a rapid decrease in the bioconvective patterns. Pattern formation has been observed in cultures of different microorganisms such as Chlamydomonas nivalis, Chlamydomonas reinhardtii, Euglena gracilis, Bacillus subtilis, Paramecium tetraurelia, and Tetrahymena pyriformis.

More recently, investigations have been reported for a semi-dilute suspension of swimming microorganisms where cell–cell interactions are considered [34–38]. On the other hand, Kitsunezaki et al. [39] investigated the effect of oxygen and depth on bioconvective patterns in suspensions with high concentrations of Paramecium tetraurelia. Bioconvection is also studied from other points of view in gravitational biology. Interesting results are also available in Refs. [40–42] about the pattern formation in suspensions of Tetrahymena and Chlamydomonas subject to different gravity conditions. Further results are due to Sawai et al. [43] who investigate the proliferation of Paramecium under simulated microgravity, to Mogami et al. [44] who report an investigation of the formed patterns by Tetrahymena and Chlamydomonas as well as a physiological comparison, to Takeda et al. [45] who give an explanation of the gravitactic behavior of single cells of Paramecium in terms of the swimming velocity and swimming direction, to Mogami et al. [46] who present theory and experiments of two mechanisms of gravitactic behavior for microorganisms, and to Itoh et al. [47] who investigate the modification of bioconvective patterns under strong gravitational fields.

This chapter presents interesting results about the bioconvective linear stability of a suspension of swimming microorganisms. Use is made of the equations presented by Ghorai and Hill [48, 49] some years ago. In their approach, Ghorai and Hill [48, 49] used a different dimensionalization scale for the concentration microorganisms which gives distinct meaning to the basic state for the concentration of microorganisms and a bioconvective Rayleigh number defined in terms of the mean cell

Bioconvective Linear Stability of Gravitactic Microorganisms DOI: http://dx.doi.org/10.5772/intechopen.83724

concentration. To the authors best knowledge, those equations along with the change in the basic state and Rayleigh number definitions have not been used to determine the linear bioconvective instability in an infinite horizontal fluid layer and to compare the results with experiment. These results were obtained by means of both numerical and analytical techniques. The critical values of the Rayleigh number Rc and the wavenumber kc, for fixed values of the gyrotaxis number G and the motility of microorganisms d, that characterize the hydrodynamic stability of the system are compared with the experimental data presented in Table I of Bees and Hill [27] and Table II of Bees and Hill [26] where the gyrotactic biflagellate alga Chlamydomonas nivalis is used as suspended microorganism. Below, it is shown for the first time that the numerical results have a very good agreement with the experimental data.

The chapter is organized as follows. The governing equations and boundary conditions [48, 49] as well as the basic state can be found in Section 2. Nondimensionalization and linearization of the system of equations is outlined in Section 3. In Section 4, use is made of an asymptotic expansion [50–53] method and a Galerkin method [54] to find limiting cases and predict critical values of R and k for the instability. The numerical calculations done by means of the shooting method along with the graphics corresponding to the marginal curves are given in Section 5. In Section 6, the experimental data [27] are compared with the numerical results. A discussion is given in the final section.

## 2. Equations of motion

We consider an infinite horizontal layer of a suspension of gyrotactic microor- <sup>∗</sup> ganisms. The fluid layer is bounded at <sup>z</sup> ¼ �H, 0. The fluid where the cellular microorganisms swim is water with density ρ. Each cell has a volume υ and density ρ þ Δρ, where Δρ ≪ ρ. The suspension is considered dilute and incompressible. Density fluctuations in the suspension are small enough such that the Boussinesq approximation is valid and the corresponding governing equations are

$$
\rho \frac{D\mathbf{u}^\*}{Dt^\*} = -\nabla p^\* - n^\* \text{ug} \Delta \rho \mathbf{k} + \mu \nabla^2 \mathbf{u}^\* \tag{1}
$$

$$\frac{\partial \mathbf{n}^\*}{\partial t^\*} = -\nabla \cdot \mathbf{J}^\* \tag{2}$$

$$\nabla \cdot \mathbf{u}^\* = \mathbf{0} \tag{3}$$

where t <sup>∗</sup> is the time, u<sup>∗</sup> is the suspension velocity, p<sup>∗</sup> is the pressure, gk is the acceleration due to gravity, k is the vertical unit vector, μ is the viscosity, n<sup>∗</sup> is the concentration of microorganisms, and J <sup>∗</sup> is the flux density of organisms through the fluid defined as

$$\mathbf{J}^\* = n^\* \mathbf{u}^\* + n^\* V\_c \mathbf{p}^\* - D\_c \nabla n^\* \tag{4}$$

where Vc is the cell swimming speed, p<sup>∗</sup> is a unit vector representing the average orientation of cells, and Dc is a scalar microorganism mass diffusion coefficient independent of the other parameters of the problem. Use is made of Cartesian coordinates with the z-axis in the vertical direction. The walls at z<sup>∗</sup> ¼ �H, 0 are considered to be rigid. As pointed out by Hill et al. [24] although the top boundary is open to the air, algal cells tend to collect at the surface forming what appears to be a packed layer, and it is unlikely that the boundary is ever fully stress-free. Then the boundary conditions are

Heat and Mass Transfer - Advances in Science and Technology Applications

$$\mathbf{u}^\* = \mathbf{0} \,\text{at} \,\text{z}^\* = -H, \mathbf{0} \tag{5}$$

$$\mathbf{J}^\* \cdot \mathbf{k} = \mathbf{0} \,\mathrm{at} \,\mathrm{z}^\* = -H, \,\mathbf{0} \tag{6}$$

<sup>∗</sup> In the basic state, the fluid velocity is zero and n<sup>∗</sup> <sup>¼</sup> <sup>n</sup>0ð Þ<sup>z</sup> and p0 <sup>¼</sup> <sup>k</sup>. Thus for <sup>∗</sup> <sup>n</sup>0ð Þ<sup>z</sup> from Eq. (2) with the boundary conditions (6), we have

$$m\_0^\*(\mathbf{z}) = \frac{\overline{n}V\_cH \exp\left[V\_c\mathbf{z}^\*/D\_c\right]}{D\_c(\mathbf{1} - \exp\left[-V\_cH/D\_c\right])}\tag{7}$$

where n represents the average concentration of organisms. Eq. (7) is the same basic state as presented by Ghorai and Hill [48, 49] whose linear stability will be investigated. It differs from that of Childress et al. [19] and Hill et al. [24] by the coefficient

$$\frac{V\_cH}{D\_c(1 - \exp\left[-V\_cH/D\_c\right])}$$

#### 3. Linear stability

We make the governing Eqs. (1–3) nondimensional by scaling all lengths with H, the time with H<sup>2</sup>=Dc, the fluid velocity with Dc=H, the pressure with νDcρ=H<sup>2</sup> , and the cell concentration with nHVc=Dc. Now the dimensionless variables are expressed without star. The boundaries are located at z ¼ �1, 0 and the basic state is

$$\mathbf{u\_0} = \mathbf{0}, \quad \mathbf{p\_0} = \mathbf{k}, \quad n\_0(z) = \frac{e^{dx}}{1 - e^{-d}}$$

where the nondimensional quantity d ¼ VcH=Dc is the ratio of swimming speed of microorganisms and their representative mass diffusion velocity. Here, d is called the motility of the microorganisms. In order to investigate the linear stability of the system, small perturbations have to be considered. They are

$$\mathbf{u} = \mathbf{u}\_0 + \delta \mathbf{u}\_1, \ p = p\_0 + \delta \mathbf{p}\_1, \ \mathbf{p} = \mathbf{p}\_0 + \delta \mathbf{p}\_1, \ n = n\_0 + \delta n\_1$$

where <sup>δ</sup> <sup>≪</sup> 1. The components of u1 are ðu1; <sup>v</sup>1; <sup>w</sup>1Þ. In this way, the nondimensional governing Eqs. (1–3) are linearized to order Oð Þδ . Then, we have the following linear equations

$$\mathbf{S}\mathbf{r}^{-1}\frac{\partial \mathbf{u}\_1}{\partial t} = -\nabla p - R n\_1 \mathbf{k} + \nabla^2 \mathbf{u}\_1 \tag{8}$$

$$\frac{\partial n\_1}{\partial t} = -w\_1 \frac{dn\_0}{dz} - d \frac{\partial n\_1}{\partial z} - dn\_0 \nabla \cdot \mathbf{p\_1} + \nabla^2 n\_1 \tag{9}$$

$$\nabla \cdot \mathbf{u\_1} = \mathbf{0} \tag{10}$$

where

$$\mathcal{S}\_{\mathfrak{c}} = \frac{\nu}{D\_{\mathfrak{c}}}, \quad R = \frac{\overline{n} \text{ug} \Delta \rho H^3 d}{D\_{\mathfrak{c}} \nu \rho}$$

are the Schmidt and bioconvection Rayleigh numbers, respectively. Pedley and Kessler [55] give a definition of the vector p1 for swimming microorganisms with spheroidal shape. They determine p1 in terms of u1 that in nondimensional form is Bioconvective Linear Stability of Gravitactic Microorganisms DOI: http://dx.doi.org/10.5772/intechopen.83724

$$\mathbf{p\_1} = G \left[ (\mathbf{1} + \alpha\_0) \frac{\partial \mathbf{u\_{1\perp}}}{\partial \mathbf{z}} - (\mathbf{1} - \alpha\_0) \nabla\_\perp w\_1 \right] \tag{11}$$

$$G = \frac{BD\_c}{H^2} \tag{12}$$

where the subscript ⊥ denotes the horizontal component, α0 is the cell eccentricity, and G is the nondimensional form of the gyrotactic orientation parameter B. Finally after substitution of p1 and n0, the governing equations become

$$\mathbf{Sc}^{-1}\frac{\partial \mathbf{u}}{\partial t} = -\nabla p - R n \mathbf{k} + \nabla^2 \mathbf{u} \tag{13}$$

$$\frac{\partial n}{\partial t} = \frac{d e^{dx}}{1 - e^{-d}} \left( -w + G \left[ (\mathbf{1} + \mathbf{a}\_0) \frac{\partial^2 w}{\partial x^2} + (\mathbf{1} - \mathbf{a}\_0) \nabla\_\perp^2 w \right] \right) - d \frac{\partial n}{\partial x} + \nabla^2 n \tag{14}$$

$$\nabla \cdot \mathbf{u} = \mathbf{0} \tag{15}$$

with boundary conditions

$$\mathbf{u} = \mathbf{0}\,\mathrm{at}\,\mathrm{z} = -\mathbf{1}, \,\mathbf{0} \quad \text{and} \quad \frac{\partial \mathbf{u}}{\partial \mathbf{z}} - d\boldsymbol{n} = \mathbf{0}\,\mathrm{at}\,\mathbf{z} = -\mathbf{1}, \,\mathbf{0} \tag{16}$$

where the superscripts have been deleted. Notice that the adimensionalization of the equations is different from that of Hill et al. [24]. Here, an application of a more general asymptotic analysis for any magnitude of d is used. An analytic Galerkin method and a shooting numerical method for the solution of the proper value problem allowed us to have an interesting perspective of the stability of the present problem under research. The results are used here to compare with the experimental data of the flagellated alga Chlamydomonas nivalis.

By elimination of the pressure from Eqs. (13–15), it is possible to obtain a coupled system of two equations, for w and n, to describe the instability of the system. The perturbations of the variables will be analyzed in terms of normal modes of the form

$$\begin{aligned} w &= W(z) \exp\left[ (k\_{\ge} \mathfrak{x} + k\_{\ge} \mathfrak{y})i + \mathfrak{c}t \right], \\ n &= \Phi(z) \exp\left[ (k\_{\ge} \mathfrak{x} + k\_{\ge} \mathfrak{y})i + \mathfrak{c}t \right] \end{aligned}$$

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where <sup>k</sup> <sup>¼</sup> <sup>k</sup><sup>2</sup><sup>þ</sup> <sup>k</sup>2 is the wavenumber of the disturbance and <sup>σ</sup> is the growth x y rate. The wavenumber is scaled as <sup>k</sup> <sup>¼</sup> <sup>k</sup><sup>∗</sup> H corresponding to a nondimensional wavelength λ ¼ 2π=k. Thus, the governing equations become

$$-Rk^2\Phi = \left(\frac{\sigma}{S\_c} + k^2 - D^2\right) \left(k^2 - D^2\right) W \tag{17}$$

$$\left(\sigma + dD + k^2 - D^2\right)\Phi = \frac{d\dot{\mathcal{C}}^x}{1 - e^{-d}} \left(-\mathbf{1} + G\left[ (\mathbf{1} + \alpha\_0)D^2 - (\mathbf{1} - \alpha\_0)k^2 \right] \right) W \tag{18}$$

subject to the boundary conditions

$$DW = DW = D\Phi - d\Phi = 0 \quad \text{at} \ z = -1, 0 \tag{19}$$

where D ¼ d=dz. The variables of the above problem can be changed in order to simplify the analysis. The change of dependent variable is

$$\Phi = F(z)e^{dx}$$

Then, Eqs. (17) and (18) and the boundary conditions Eq. (19) become

$$-Rk^2Fe^{dx} = \left(\frac{\sigma}{S\_c} + k^2 - D^2\right)\left(k^2 - D^2\right)W\tag{20}$$

$$\left(\sigma - dD + k^2 - D^2\right)F = \frac{d}{1 - e^{-d}} \left(-\mathbf{1} + G\left[ (\mathbf{1} + \mathbf{a}\_0)D^2 - (\mathbf{1} - \mathbf{a}\_0)k^2 \right] \right)W \tag{21}$$

subject to the new boundary conditions

$$DW = DW = DF = \mathbf{0} \text{ at } z = -\mathbf{1} \text{ } \mathbf{0} \tag{22}$$

In this form, the equations are very similar to those of the well-known problem of thermal convection in an infinite horizontal fluid layer between nonconducting boundaries [50–53, 56]. The familiar fixed heat flux boundary condition is the main characteristic of those thermal convection problems and is analogous to that presented in Eq. (22). The equations derived by Childress et al. [19] can also be analyzed from the present view point of this change of variable. In the theory of thermal convection as in that of Childress et al. [19], a zero critical wavenumber is found as a result of the fixed flux boundary condition. In more recent models, which include the effects of gyrotaxis, the similarity with the thermal convection problem is not valid unless G ¼ 0.

#### 4. Asymptotic analysis

In this section, the eigenvalue problem stated in the system of Eqs. (13–15) with boundary condition Eq. (16) is investigated by means of two analytic methods. The magnitude of the marginal value of R is a function of all the other parameters. The way in which the solution of the stability problem is to be carried out is as follows. For a given value of d and G, we must determine the lowest value for R with respect to the wavenumber k. The values obtained are the critical Rayleigh numbers Rc at which instability will first occur.

#### 4.1 Asymptotic analysis

We conducted a general asymptotic analysis in comparison with those used before [19, 24, 26] which included the restrictions of the limits d ≪ 1 for shallow layers and d≫1 for deep layers along with different restrictions for G. In a similar way, as in other problems of convection, we follow the steps of Chapman and Proctor [51], Dávalos-Orozco [52], and Dávalos-Orozco and Manero [53]. Under the above conditions, the analysis is very complex, the reason why use has been made of the Maple algebra package. Thus, we look for a solution to Eqs. (20) and (21) using the following expansions:

$$\mathcal{W} = \mathcal{W}\_0 + \varepsilon \mathcal{W}\_1 + \dots,\tag{23}$$

$$
\Phi = \Phi\_0 + \varepsilon \Phi\_1 + \dots \tag{24}
$$

$$R = R\_0 + \varepsilon R\_1 + \dots \tag{25}$$

$$
\sigma = \varepsilon \sigma\_0 + \varepsilon^2 \sigma\_1 + \dots \tag{26}
$$

where ε ≪ 1. We also consider no restrictions for d and G and rescale the wavenumber as <sup>k</sup> <sup>¼</sup> <sup>ε</sup><sup>1</sup>=2<sup>~</sup> k. Thus, after substitution of expansions Eqs. (23–26) and the mentioned scalings in Eqs. (20) and (21) with boundary condition Eq. (22), we obtain the following systems of equations at different orders.

At order Oð Þ1

$$\left(D^4W\_0 + \tilde{k}^2 R\_0 e^{dx} F\_0 = \mathbf{0}, \tag{27}$$

$$(\mathbf{J}^2 + d)\mathbf{J}'\_0 = \mathbf{0} \tag{28}$$

subject to

$$DW\_0 = DW\_0 = DF\_0 = \mathbf{0} \text{ at } z = -\mathbf{1} \text{ } \mathbf{0}.\tag{29}$$

At order Oð Þε

$$D^4 W\_1 - \left(\tilde{\mathcal{M}}^2 + \frac{\sigma\_0}{\text{Sc}}\right) \not{p}^2 W\_0 + \tilde{k}^2 e^{\text{d}x} (R\_0 F\_1 + R\_1 F\_0) = \mathbf{0},\tag{30}$$

$$[G(\mathbf{1} + \alpha\_1)\mathbf{D}^2 - \mathbf{1}] \mathcal{W}\_c - \left(\not{\mathbf{1}}\_c + \tilde{k}^2\right) \not{p}^\prime\_c + (D^2 + d) \mathcal{H}\_1 = \mathbf{0}.\tag{31}$$

$$\frac{d}{\mathbf{1} - \mathbf{e}^{-d}} \left[ G(\mathbf{1} + \mathbf{a}\_0) D^2 - \mathbf{1} \right] \mathbf{W}\_0 - \left( \mathbf{e}\_0^\prime + \tilde{\mathbf{k}}^2 \right) \mathbf{f}\_0^\prime + (D^2 + d) \mathbf{f}\_1^\prime = \mathbf{0} \tag{31}$$

subject to

$$DW\_1 = DW\_1 = DF\_1 = 0 \quad \text{at} \ z = -1, 0. \tag{32}$$

At order <sup>O</sup> <sup>ε</sup><sup>2</sup> ð Þ

$$D^4 \mathcal{W}\_2 - \left(\tilde{\mathcal{K}}^2 + \frac{\sigma\_0}{\mathcal{Sc}}\right) \tilde{\mathcal{C}}^2 \mathcal{W}\_1 + \left.\tilde{k}\right|^4 + \frac{\sigma\_0 \tilde{k}^2}{\mathcal{Sc}} - \frac{\sigma\_1}{\mathcal{Sc}} \tilde{D}^2 \Big/ \left(\mathcal{W}\_0 + \tilde{k}^2 e^{\mathrm{d}x} (R\_0 F\_2 + R\_1 F\_1 + R\_2 F\_0)\right) = \mathbf{0},\tag{33}$$

$$\frac{d}{1 - e^{-d}} \left\{ \left[ G(1 + a\_0) D^2 - 1 \right] \mathcal{W}\_1 - G(1 + a\_0) \bar{k}^2 \mathcal{W}\_0 \right\} - \left( \sigma\_0 + \bar{k}^2 \right) F\_1 - \sigma\_1 F\_0 + \left( D^2 + d \right) \mathcal{Y}\_2' = 0 \tag{34}$$

subject to

$$DW\_2 = DW\_2 = DF\_2 = 0 \quad \text{at} \ z = -1, 0. \tag{35}$$

The systems of equations at order Oð Þε and higher are inhomogeneous and must satisfy their corresponding solvability conditions allowing to compute the Rayleigh number R as an eigenvalue in terms of the other parameters of the problem. Solvability conditions are found as usual [57]: each inhomogeneous system is multiplied by the solution to the adjoint of the homogeneous system and integrated over the range of the independent variable. The resulting integral must vanish.

Thus, the solvability conditions at Oð Þ<sup>ε</sup> and <sup>O</sup> <sup>ε</sup><sup>2</sup> ð Þ are, respectively,

$$\begin{split} \mathbf{0} &= \frac{d}{\mathbf{1} - \mathbf{c}^{-d}} \int\_{-1}^{0} e^{d\mathbf{z}} \left[ G(\mathbf{1} + \mathbf{a}\_{0}) D^{2} - \mathbf{1} \right] \mathbf{W}\_{0} d\mathbf{z} - \left( \mathbf{q}\_{0} + \bar{\mathbf{k}}^{2} \right) \int\_{-1}^{0} e^{d\mathbf{z}} F\_{0} d\mathbf{z} \\ \mathbf{0} &= \frac{d}{\mathbf{1} - \mathbf{c}^{-d}} \int\_{-1}^{0} e^{d\mathbf{z}} \left\{ \left[ G(\mathbf{1} + \mathbf{a}\_{0}) D^{2} - \mathbf{1} \right] \mathbf{W}\_{1} - G(\mathbf{1} + \mathbf{a}\_{0}) \bar{\mathbf{k}}^{2} W\_{0} \right\} d\mathbf{z} \\ &- \left( \mathbf{q}\_{0} + \bar{\mathbf{k}}^{2} \right) \int\_{-1}^{0} e^{d\mathbf{z}} F\_{1} d\mathbf{z} - \sigma\_{1} \int\_{-1}^{0} e^{d\mathbf{z}} F\_{0} d\mathbf{z} \end{split} \tag{37}$$

The solutions of the system of equations at leading order are

Heat and Mass Transfer - Advances in Science and Technology Applications

$$F\_0 = \mathbf{1}, \ \ W\_0 = f\_1(z, d) \mathcal{R}\_0 \tilde{k}^2 \tag{38}$$

where the function f <sup>1</sup>ðz; dÞ can be obtained from the authors upon request. For convenience, the solution F0 has been normalized to 1. The next step is to evaluate the solvability condition Eq. (36) at Oð Þε and obtain an expression for σ<sup>0</sup>

$$
\sigma\_0 = \tilde{k}^2 \left\{ \left[ \mathbf{1} - d^2 G (\mathbf{1} + \alpha\_0) \right] A\_0 R\_0 - \mathbf{1} \right\} \tag{39}
$$

The constant A0 is large and can be obtained from the authors upon request. At order Oð Þε similar steps as those for solving the system of equations at Oð Þ1 are followed to find F1 and W1. Then, algebraically F1 is

$$F\_1 = \left[ G(\mathbf{1} + \mathbf{a}\_0) f\_2(z, d) + f\_3(z, d) \right] R\_0 \tilde{k}^2 \tag{40}$$

After substitution of F0, W0, F1, and σ0 into Eq. (30), the velocity W1 can be calculated subject to its corresponding boundary condition Eq. (32). Because of the term edz appearing in the system of equations at Oð Þ1 , the expression of W1 is very large and complicated and will not be given here. The evaluation of the solvability condition at order <sup>O</sup> <sup>ε</sup><sup>2</sup> ð Þ given in Eq. (37) yields

$$\begin{aligned} \sigma\_1 &= \left[1 - d^2 G(\mathbf{1} + \mathbf{a}\_0)\right] \bar{k}^2 A\_0 \mathbf{R}\_1 - \left\{ \left[ G^2 (\mathbf{1} + \mathbf{a}\_0)^2 A\_1 + G(\mathbf{1} + \mathbf{a}\_0) A\_2 + A\_3 \right] \mathbf{R}\_0^2 \right. \\\\ &+ G(\mathbf{1} - \mathbf{a}\_0) A\_4 \mathbf{R}\_0 + \mathbf{R}\_0 \left[ \mathbf{1} - d^2 G(\mathbf{1} + \mathbf{a}\_0) \right] \left\{ \left[ G(\mathbf{1} + \mathbf{a}\_0) A\_5 + A\_6 \right] \mathbf{R}\_0 \right. \\\\ &+ A\_7 + \left( \left[ \mathbf{1} - d^2 G(\mathbf{1} + \mathbf{a}\_0) \right] A\_8 \mathbf{R}\_0 + A\_7/2 \right) \bar{k}^4 \end{aligned} \tag{41}$$

The growth rate can now be obtained by substitution of σ0 and σ1 into the expansion for σ given in Eq. (26). However σ is omitted to save space but can be obtained from the authors upon request. Finally, use is made of the expansion Eq. (25) for R.

Now, the transition from stability to instability via a stationary state is investigated by setting σ ¼ 0. Thus, the corresponding value of R for the marginal state is

$$R = \frac{1}{\left[1 - d^2 G(1 + a\_0)\right] A\_0} + \frac{k^2}{\left[1 - d^2 G(1 + a\_0)\right] A\_0}$$

$$\left\{ \frac{G(1 - a\_0)}{\left[1 - d^2 G(1 + a\_0)\right]} + \frac{G(1 + a\_0) A\_5 + A\_6}{\left[1 - d^2 G(1 + a\_0)\right] A\_0^2} \right. \\ \left. + \frac{G^2(1 + a\_0)^2 A\_1 + G(1 + a\_0) A\_2 + A\_3}{\left[1 - d^2 G(1 + a\_0)\right]^2 A\_0^2} + \frac{A\_7}{A\_0} \right\} \right\} \tag{42}$$

where some simplifications have been made with the use of R0 obtained from Eq. (39). The functions f <sup>2</sup>ðz; dÞ and f <sup>3</sup>ðz; dÞ and the coefficients A1 to A8 appearing in the above expressions are functions of d and can be obtained from the authors upon request. The result for R0 is

$$R\_0 = \frac{1}{\left[1 - d^2 G (1 + a\_0)\right] A\_0} \tag{43}$$

From the expression for the Rayleigh number given in Eq. (42), it is possible to calculate the limit for d ≪ 1. In this case, we consider d and k to be of the same order in such a way that kd ¼ k=d and kd � Oð Þ1 . Then, under these assumptions, the approximation of R kð ; d; G; α0Þ is

Bioconvective Linear Stability of Gravitactic Microorganisms DOI: http://dx.doi.org/10.5772/intechopen.83724

$$R = 720\left\{ 1 + d^2 \left[ \frac{17}{420} + G(1 + a\_0) + k\_d^2 \left( \frac{17}{462} - \frac{2}{7} G(5 - 2a\_0) \right) \right] \right\} + O\left(d^3\right) \tag{44}$$

Here we point out that in the present chapter, our definition of the Rayleigh number differs from that defined by Hill et al. [24]. If our approximation given in ˙ <sup>ˆ</sup> <sup>d</sup> 2 4 Eq. (44) is multiplied by �<sup>d</sup> <sup>¼</sup> <sup>1</sup><sup>þ</sup> <sup>d</sup>=2 <sup>þ</sup> <sup>d</sup> <sup>=</sup>12 � <sup>d</sup> <sup>=</sup>720… , the same approxi- <sup>1</sup>�<sup>e</sup> mation by Hill et al. [24] is obtained. Moreover, if G ¼ 0 this becomes that given by Childress et al. [19]. In the more general expression of R for a small wavenumber approximation, Eq. (42) has a special characteristic due to its dependence on the square of the wavenumber k. The first coefficient at zeroth order in k corresponds to R0, and that at the second order in k is R1. Even though in the experiments on bioconvection [27], only finite critical wavenumbers kc > 0 have found these coefficients are very useful. For example, it can be shown from R0 that if A<sup>0</sup> > 0 and

$$1 - d^2 G(1 + a\_0) > 0\tag{45}$$

then R0 < 0. This corresponds to a stable stratification, which is not the case here. The second coefficient R1 in Eq. (42) is a very important result, because it provides information about the shape of the marginal curve with respect the critical wavenumber. That information can be obtained by making zero the coefficient R<sup>1</sup> and calculating the following critical value of the gyrotaxis number Gc

$$\begin{aligned} \mathbf{0} &= G^2 (\mathbf{1} + \mathbf{a}\_0)^2 \left[ A\_1 - d^2 A\_5 - d^2 A\_0^2 \frac{\mathbf{1} - \mathbf{a}\_0}{\mathbf{1} + \mathbf{a}\_0} \right] + G(\mathbf{1} + \mathbf{a}\_0) \\ &\quad \left[ A\_2 + A\_5 - d^2 (A\_6 + A\_7 A\_0) + A\_0^2 \frac{\mathbf{1} - \mathbf{a}\_0}{\mathbf{1} + \mathbf{a}\_0} \right] + A\_3 + A\_6 + A\_7 A\_0 \end{aligned} \tag{46}$$

From this equation, two admissible cases are possible when Eq. (45) is satisfied. First, for fixed values of d, α0, and G > Gc, the marginal curve starts at k ¼ 0 and then decreases monotonically. However, according to the numerical analysis presented below, the marginal curves in fact first decrease and then start to grow monotonically after a minimum is attained, at the critical wavenumber. In the second case, for fixed values of d, α0, and G < Gc, the marginal curves start at k ¼ 0 and then grow monotonically. Here, the critical wavenumber is always zero. The importance of these results is that the magnitude of Gc agrees very well with the results of the marginal curves found in the numerical analysis given below. This critical value determines the magnitude for which the curves have a finite critical wave number (G > Gc) or a zero critical wavenumber (G < Gc). It is of interest to note that some of the magnitudes of the gyrotaxis parameter G calculated from the data in the literature are very near but above of Gc. This is the reason why some of the curves found from numerical analysis are almost flat in a range of wavenumbers near to zero. Because the experimental critical wavenumbers found for gyrotactic bioconvection are always finite, we conclude that Gc is important to point out where to find the theoretical limitations of the model.

#### 4.2 Analytic Galerkin method

Here use is made of the analytical Galerkin method to study the eigenvalue problem of Eqs. (17)–(18) with the boundary condition Eq. (19). This method has been used before by Pellew and Soutwell [58], Chandrasekhar [54], and Gershuni and Zhukovitskii [59]. Even though this is an approximate method, it has a very high precision. The advantage of the method is that it is possible to obtain an explicit expression of the Rayleigh number R. Here, it is supposed that σ ¼ 0.

Briefly, the method consists in assuming a trial function which satisfies the boundary conditions for each of the dependent variables. Let that variable be Φ which after substitution in one of the equations of the problem allows for the exact solution of the other variables, let us say W. Both trial functions are now substituted into the other coupled equation. Then, use is made of the orthogonality properties of the solutions in this equation to obtain the proper value of the Rayleigh number as a function of the other parameters [60].

In this way, the proposed expansions of Φ and W are

$$\Phi = \sum\_{n=0}^{\infty} B\_n \exp\left(dz\right) \cos \pi nz \quad \text{and} \quad W = \sum\_{n=0}^{\infty} B\_n W\_n \tag{47}$$

then, after substitution of Φ of Eq. (47) into Eq. (20), Wn is the solution of the following differential equation:

$$\left(\left(\text{D}^2 - k^2\right)^2 W\_n = -Rk^2 \exp\left(dz\right) \cos \pi nz \tag{48}$$

which is subjected to the conditions

$$W\_n = DW\_n = \mathbf{0} \text{ at } z = -\mathbf{1}.\text{0}.\tag{49}$$

The solution is

$$\mathcal{W}\_n = (a\_1 \cos \pi nz + a\_2 \sin \pi nz)e^{kz}k^2R + c\_1e^{kz} + c\_2e^{-kz} + c\_3ze^{kz} + c\_4ze^{-kz} \tag{50}$$

where c1 to c4 are constants of integration which can be found by evaluating in the boundary condition Eq. (49).

Next, Eq. (18) is multiplied by Φ and is integrated in the range z ¼ �1 to z ¼ 0, to get

$$\begin{aligned} \frac{d}{1 - e^{-d}} \int\_{-1}^{0} \Phi \left( -\mathbf{1} + G \left[ (\mathbf{1} + \mathbf{a}\_0) D^2 - (\mathbf{1} - \mathbf{a}\_0) k^2 \right] \right) W dz \\ - \int\_{-1}^{0} \Phi \left[ D \left( e^{-dx} D \right) - k^2 e^{-dx} \right] \Phi dz = \mathbf{0} \end{aligned} \tag{51}$$

After substitution of Φ and W, given in Eqs. (47) and (50) and some simplifications, we obtain

$$\begin{aligned} & \left| \frac{d}{\mathbf{1} - e^{-d}} \right|\_{-1}^{0} \Phi\_m \left( -\mathbf{1} + G \left[ (\mathbf{1} + \mathbf{a}\_0) D^2 - (\mathbf{1} - \mathbf{a}\_0) k^2 \right] \right) \mathcal{W}\_n dz \\ & - \int\_{-1}^{0} \Phi\_m \left[ D \left( e^{-dx} D \right) - k^2 e^{-dx} \right] \Phi\_n dz \right| = \mathbf{0} \end{aligned} \tag{52}$$

This determinant, calculated with the help of the software Maple, is the solvability condition from which the eigenvalue R is calculated. The resulting algebraic expression of the integrals in this equation is very complex and will not be presented here. However, the first approximation of R, corresponding to the element (0,0) of the matrix in the determinant Eq. (52), is

$$R = \frac{(\mathbf{1} - e^{-d})A\_{10}}{de^{-d} \left[\mathbf{1} + G(\mathbf{1} - \mathbf{a}\_0)k^2 - d^2G(\mathbf{1} + \mathbf{a}\_0)\right]A\_9} \tag{53}$$

where A9 and A10 are functions of the wavenumber k and d and can be obtained from authors upon request. This result is new because it includes, for the first time,

#### Bioconvective Linear Stability of Gravitactic Microorganisms DOI: http://dx.doi.org/10.5772/intechopen.83724

all the parameters of the problem without any approximation. In the limit of d, k, G ! 0, R reduces to the well-known value of 720. Higher-order estimates of R can be obtained from Eq. (52), which provides a useful check on numerical calculations. The comparison of R given in Eq. (42), obtained from the asymptotic analysis, and that of Eq. (52) shows that in the limit k ! 0, the agreement was very good.
