2.1. Governing equations

We consider two-dimensional (2-D) bubbly flows in a horizontal fluid layer. Small spherical bubbles are injected into the liquid layer from the bottom wall (figure 1$a$). The injection is assumed uniform in space and constant in time. Injected bubbles rise in the layer due to buoyancy forces and are ejected from the top free surface. We model the dynamics of this immiscible two-phase flow system using the Euler–Euler approach, in which the liquid and gas phases are regarded as two dynamical continua exchanging momentum and energy with each other. Mass and momentum conservations for the two phases are modelled as

(2.1a)$$\begin{gather} \frac{\partial \varepsilon_G}{\partial t} +\boldsymbol{\nabla} \boldsymbol{\cdot} (\varepsilon_G \boldsymbol{u}_G)=0, \end{gather}$$

(2.1b)$$\begin{gather}\frac{\textrm{D} \boldsymbol{u}_G}{\textrm{D}t} = 3\frac{\textrm{D} \boldsymbol{u}}{\textrm{D} t}-\frac{18}{r_b^{2}}(\boldsymbol{u}_G-\boldsymbol{u})-(\boldsymbol{u}_G-\boldsymbol{u}) \times (\boldsymbol{\nabla} \times \boldsymbol{u})+\frac{18}{Pr r_b^{2}} \boldsymbol{e}_z, \end{gather}$$

(2.1c)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$

(2.1d)$$\begin{gather}\frac{\textrm{D} \boldsymbol{u}}{\textrm{D} t} ={-}\boldsymbol{\nabla} p + \Delta \boldsymbol{u}-\frac{9}{Pr r_b^{2}} \boldsymbol{e}_z+\varepsilon_G\frac{Ra r_b^{2}}{9Pr} \left(\frac{\textrm{D} \boldsymbol{u}}{\textrm{D} t}+\frac{9}{Pr r_b^{2}} \boldsymbol{e}_z \right), \end{gather}$$

where $\varepsilon _G$ is gas volume fraction, $\boldsymbol {u}_G = u_G\boldsymbol {e}_x + w_G\boldsymbol {e}_z$ and $\boldsymbol {u} = u\boldsymbol {e}_x + w\boldsymbol {e}_z$ are the velocity fields of gas and liquid, and $p$ is the pressure field (see Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020, and references therein for details). The unit vectors in the $x$ and $z$ directions are denoted by $\boldsymbol {e}_x$ and $\boldsymbol {e}_z$. The differential operator $\textrm {D}/\textrm {D}t$ represents $\textrm {D}{\boldsymbol u}_G/\textrm {D}t = \partial _t {\boldsymbol u}_G + {\boldsymbol u}_G \boldsymbol {\cdot } \boldsymbol {\nabla } {\boldsymbol u}_G$ and $\textrm {D}{\boldsymbol u}/\textrm {D}t = \partial _t {\boldsymbol u} + {\boldsymbol u} \boldsymbol {\cdot } \boldsymbol {\nabla }{\boldsymbol u}$. These equations have been non-dimensionalised with scales $d$ of length, $\nu / d$ of velocity and $d^{2}/\nu$ of time. The gas volume fraction has been normalised with a scale $J_{G,0} d /\nu$, where $J_{G,0}$ is the gas injection flux. We have introduced the following four dimensionless parameters:

(2.2a–d)\begin{align} Pr = \frac{\nu}{V_\infty d}= \frac{9\nu^{2}}{g R_b^{2} d}, \quad Ra =\frac{gd^{2} J_{G,0}}{V_\infty^{2} \nu}= \frac{81\nu d^{2} J_{G,0}}{gR_b^{4}}, \quad w_{G,0}=\frac{W_{G,0}\, d}{\nu}, \quad r_b=\frac{R_b}{d}, \end{align}

where $R_b$ and $W_{G,0}$ are the radius and the injection velocity of bubbles, respectively. For bubbles of radius 0.5 mm injected into a water layer of thickness 100 mm at room temperature with $W_{G,0}=20\ \ \textrm {mm}\,\textrm {s}^{-1}$ and $J_{G,0}=1\ \textrm {mm}\,\textrm {s}^{-1}$, these parameters take values of $Pr=7.34\times 10^{-5}$, $Ra = 330$, $w_{G,0}=1000$, $r_b = 0.01$. Throughout the present work, the dimensionless bubble injection velocity $w_{G,0}$ and bubble radius $r_b$ are fixed at these values. In the gas momentum equation (2.1b), the hydrodynamic diffusion has been omitted, assuming dilute bubbly flows. We have also assumed spherical and non-deformable bubbles in a pure liquid without contamination by surfactants to model drag and shear-induced lift forces and adopt the drag, added-mass and lift coefficients $C_D = 48/{Re}_b$ (Levich Reference Levich1962), where $Re_b=2R_b \|{\bf u}_G - {\bf u}\| / \nu$ is the bubble Reynolds number, and $C_A=C_L = 1/2$ (Magnaudet & Eames Reference Magnaudet and Eames2000). This model of $C_D$ describes well bubble motions for ${Re}_b \sim 30 - 400$ (Clift, Grace & Weber Reference Clift, Grace and Weber2005), e.g. bubbles of $R_b = 0.5\ \textrm {mm}$ in a steady ascending motion in water at room temperature. In the liquid mass conservation equation (2.1c) and the liquid momentum equation (2.1d), liquid flows are considered effectively incompressible, as we have assumed dilute bubbly flows. We have also restricted our attention to flows at scales larger than the mean distance between bubbles. The dynamics of the liquid phase is thus affected by bubbles only through the mesoscale reaction force, $\varepsilon _G(\textrm {D}{\boldsymbol u}/\textrm {D}t-\boldsymbol {g})$ (Druzhinin & Elghobashi Reference Druzhinin and Elghobashi1998).

The upper surface of the liquid layer is assumed as flat and shear-free for simplicity. At the bottom wall, no-slip conditions on the liquid velocity, constant gas velocity and constant gas flux are imposed:

(2.3a)$$\begin{gather} \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} = 0, \quad w=0, \quad \text{at } z=1, \end{gather}$$

(2.3b)$$\begin{gather}{\boldsymbol u}={\boldsymbol 0}, \quad w_{G}=w_{G,0},\quad \varepsilon_{G}=1, \quad \text{at } z={-}1. \end{gather}$$

2.2. Base state

Two-dimensional flows would respect the translational symmetry of the system along the $x$ direction when the gas flux $J_{G,0}$, or the Rayleigh number $Ra$ in the dimensionless description, is small. We thus assume a steady bubbly flow in the homogeneous regime, where the flow fields are laterally uniform:

(2.4a–c)\begin{equation} \boldsymbol{u} = \boldsymbol{0}, \quad \boldsymbol{u}_G = \bar{w}_G(z)\, \boldsymbol{e}_z, \quad \varepsilon_G = \bar{\varepsilon}_G(z). \end{equation}

For these fields, (2.1a) and (2.1b) read

(2.5a,b)\begin{equation} \frac{\textrm{d}}{\textrm{d}z}(\bar{\varepsilon}_G \bar{w}_G)=0, \quad \bar{w}_G \frac{\textrm{d} \bar{w}_G}{\textrm{d}z}={-}\frac{18 \bar{w}_G}{r_b^{2}} + \frac{18}{Pr r_b^{2}}, \end{equation}

and the boundary conditions (2.3) are reduced to

(2.6a,b)\begin{equation} \bar{w}_G=w_{G,0},\quad \bar{\varepsilon}_G = 1, \quad \text{at }z={-}1. \end{equation}

Some profiles of velocity $\bar {w}_G$ and gas fraction $\bar {\varepsilon }_G$ calculated from (2.5a,b) and (2.6a,b) are shown in figure 1$(b)$. The gas velocity increases toward the terminal velocity $V_{\infty }$ in a fluid sublayer attached to the wall. Inside the sublayer, the fraction $\bar {\varepsilon }_G$ decreases sharply so that the stratification in the mean density of the liquid–gas mixture is potentially unstable to the gravity.

The profiles of $\bar {w}_G$ and $\bar {\varepsilon }_G$ depend on $Pr$. At small $Pr$, the density gradient is small and the unstable sublayer extends over a thick zone of the fluid layer. Our linear stability analysis (Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020) showed that the whole-layered convection rolls develop when $Ra$ exceeds a critical value (figure 1$c$-i) for $Pr$ smaller than a threshold $Pr^{*}$. For $Pr$ larger than $Pr^{*}$, the unstable sublayer is thin, and the multi-layered mode (figure 1$c$-ii) is selected as a critical state. The value of $Pr^{*}$ is $3.25 \times 10^{-5}$ for $(r_b, w_{G,0})=(0.01, 1000)$. At this Prandtl number, both whole- and multi-layered modes are critical at $Ra=645$. Similar formation of whole- and multi-layered convection rolls sensitive to density profiles is reported for thermal convections (Ogura & Kondo Reference Ogura and Kondo1970; Nield Reference Nield1975).

2.3. Perturbation analysis

We consider perturbations around the base state:

(2.7a–c)\begin{equation} \boldsymbol{u}=\boldsymbol{u}', \quad \boldsymbol{u}_G=\bar{w}_G \boldsymbol{e}_z+\boldsymbol{u}_G', \quad \varepsilon_G=\bar{\varepsilon}_G + \varepsilon_G', \end{equation}

where perturbation components are indicated by primes. The velocity fields $\boldsymbol {u}'$ and $\boldsymbol {u}_G'$ are assumed 2-D: $\boldsymbol {u}' = u' \boldsymbol {e}_x + w' \boldsymbol {e}_z$, $\boldsymbol {u}_G' = u'_G \boldsymbol {e}_x + w'_G \boldsymbol {e}_z$. Substituting (2.7a–c) in (2.1) and (2.3), we have

(2.8a)$$\begin{gather} \frac{\partial \varepsilon_G'}{\partial t} +\left(\bar{w}_G \frac{\partial}{\partial z}+\frac{\textrm{d}\bar{w}_G}{\textrm{d}z} \right)\varepsilon_G' +\bar{\varepsilon}_G \frac{\partial u_G'}{\partial x}+\left(\bar{\varepsilon}_G \frac{\partial}{\partial z}+\frac{\textrm{d} \bar{\varepsilon}_G}{\textrm{d}z} \right)w_G'+\boldsymbol{\nabla} \boldsymbol{\cdot} (\varepsilon_G' \boldsymbol{u}_G')=0, \\ \hskip-16.3pc \frac{\partial \boldsymbol{u}_G'}{\partial t}+\bar{w}_G \frac{\partial \boldsymbol{u}_G'}{\partial z} +\frac{\textrm{d}\bar{w}_G}{\textrm{d}z}w_G' \boldsymbol{e}_z + \boldsymbol{u}_G' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_G' \nonumber\\ =3\left(\frac{\partial \boldsymbol{u}'}{\partial t}+\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'\right) -\frac{18}{r_b^{2}}(\boldsymbol{u}_G' - \boldsymbol{u}') +\,(\bar{w}_G+w_G'-w') \zeta' \boldsymbol{e}_x -(u_G'-u')\zeta' \boldsymbol{e}_z,\end{gather}$$

(2.8b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}' = \boldsymbol{0}, \end{gather}$$

(2.8c)$$\begin{gather}\frac{\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}' ={-}\boldsymbol{\nabla} {\rm \pi}'+\Delta \boldsymbol{u}' + \frac{Ra r_b^{2}}{9Pr} (\bar{\varepsilon}_G+\varepsilon_G') \left(\frac{\partial \boldsymbol{u}'}{\partial t} + \boldsymbol{u}'\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}'\right)+\varepsilon_G' \frac{Ra}{Pr^{2}}\boldsymbol{e}_z, \end{gather}$$

(2.9a)$$\begin{gather} \displaystyle \frac{\partial u'}{\partial z} + \frac{\partial w'}{\partial x} = w'=0, \quad \text{at }z=1, \end{gather}$$

(2.9b)$$\begin{gather}\displaystyle u'=w'=w_{G}'=\varepsilon_G'=0, \quad \text{at }z={-}1, \end{gather}$$

where $\zeta'= \partial_z u' - \partial_x w'$ is the liquid vorticity perturbation and ${\rm \pi}'$ is the perturbation of a reduced pressure. To solve the set of (2.8) and (2.9), we assume travelling wave modes and expand the perturbations $(\boldsymbol {u}_G', \boldsymbol {u}', \varepsilon _G')$ into the Fourier–Chebyshev series:

(2.10)\begin{equation} \left( \boldsymbol{u}_G',\, \boldsymbol{u}',\, \varepsilon_G' \right) = \sum_{m={-}\infty}^{\infty}\sum_{l=0}^{\infty} \left( \widehat{\boldsymbol{u}}_{G,m,l},\, \widehat{\boldsymbol{u}}_{m,l},\, \widehat{\varepsilon}_{G,m,l} \right) T_l(z) \exp({\textrm{i}m\alpha (x - c t)}), \end{equation}

where $T_l(z)$ is the Chebyshev polynomial of degree $l$, $\alpha$ is the wavenumber of fundamental mode and $c$ is the phase velocity. The case of $c=0$ corresponds to a steady convection. Substituting (2.10) for perturbation fields in (2.8a)–(2.9b) and truncating the series at $m=\pm M$ and $l=L$, we obtain a set of coupled nonlinear algebraic equations for $\widehat {\boldsymbol {u}}_{G,m,l} = (\hat {u}'_{G,m,l}, \hat {w}'_{G,m,l})$, $\widehat {\boldsymbol {u}}_{m,l} = (\hat {u}'_{m,l}, \hat {w}'_{m,l})$, and $\widehat {\varepsilon }_{G,m,l}$ ($m=0, \pm 1, \pm 2, \dots, \pm M$; $l=0, 1, \dots, L$):

(2.11)\begin{equation} F_ i = L_{ij}X_j + N^{(2)}_{ijk}X_jX_k + N^{(3)}_{ijkl}X_jX_kX_l + i \alpha c\, C_{ij} X_j = 0, \end{equation}

where $\{X_j\} =\{ \hat {u}'_{G,m,l}, \hat {v}'_{G,m,l}, \hat {u}'_{m,l}, \hat {w}'_{m,l}, \widehat {\varepsilon }_{G,m,l}\, |\, m=-M, \dots , M; l = 0, 1, \dots , L \}$ and $L_{ij}$, $N^{(2)}_{ijk}$, $N^{(3)}_{ijkl}$, and $C_{ij}$ are coefficients. We solve these equations by the numerical continuation with varying $Ra$ from the critical value (Dijkstra et al. Reference Dijkstra2014). To determine the solution at a given $Ra$, the Newton–Raphson iterative scheme is invoked with a convergence criterion ${\varDelta } =\max _j ({\varDelta }_j) < 10^{-6}$, where ${\varDelta }_j$ is the relative improvement to $(I-1)$th estimate $X_j^{I-1}$ to $X_j$ by the $I$th iteration (Deguchi & Nagata Reference Deguchi and Nagata2011). The improvement ${\varDelta }_j$ is defined as ${\varDelta }_j = | (X_j^{I} - X_j^{I-1}) / X_j^{I-1}|$ when $|X_j^{I}|, |X_{j}^{I-1}| > 10^{-6}$; ${\varDelta }_j=0$ otherwise. For the solutions computed under this criterion, (2.11) are satisfied with a tolerance of $O(10^{-8})$. The truncation parameters are fixed at $M = 25$ and $L = 35$ for convergence.

3.1. Bifurcation diagram

We compute the flows developing from critical modes with varying $Ra$ for different values of $Pr$ in the range $0.56 \leq Pr / Pr^{*} \leq 2.3$, which was considered in Nakamura et al. (Reference Nakamura, Yoshikawa, Tasaka and Murai2020) for the linear stability analysis. The obtained bifurcation diagrams for a case of $Pr < Pr^{*}$ and another case of $Pr > Pr^{*}$ are shown in figure 2, where the behaviour of the solution is monitored through the variation of $\max {(u')}$ as a function of the bifurcation parameter $\delta = (Ra - Ra_c)/Ra_c$. Determined flows are all stationary ($c=0$). Flow fields at different points on the bifurcation curves are also shown. For the whole-layered mode ($Pr < Pr^{*}$), the convection develops through a subcritical bifurcation (figure 2a). The curve revolves once $\delta$ goes down to $-0.3$ and, then, $\delta$ starts increasing with $\max {(u')}$. For the multi-layered mode ($Pr > Pr^{*}$), in contrast, the flow grows through a supercritical bifurcation (figure 2$b$). A revolution occurs at approximately $\delta = 0.02$ and, then, $\max {(u')}$ continues to increase but with decreasing $\delta$. In both cases, pairs of equally sized convection rolls develop at weak nonlinearity, i.e. at small $\max {(u')}$, as shown in figure 2($a1$,$b1$,$b2$). Sharp ascending plumes are observed at large $\max {(u')}$, see figure 2($a2$,$a3$,$b3$). These plumes are coincident with bubble-concentrated zones and delimited laterally by bubble-poor zones.

Figure 2.

Bifurcation diagrams for ($a$) the whole-layered mode and ($b$) the multi-layered mode: ($a$1–$a$3; $b$1–$b$3) flow patterns at different points in the diagrams: the bifurcation parameter $\delta$ is defined by $\delta =(Ra-Ra_c)/Ra_c$ with $Ra_c=337$ and $764$ for ($a$) and ($b$), respectively. The ordinates $\max {u'}$ are the maximum horizontal velocity of the liquid phase. In panels ($a$1–$a$3) and ($b$1–$b$3), the liquid velocity fields $\boldsymbol {u}'$ are shown by arrows and streamlines. The gas volume fraction fields $\varepsilon _G'$ are shown by colours.

We have examined the effects of the drag force model on bifurcations by determining bifurcation diagrams with $C_D=16/Re_b$ (Hadamard Reference Hadamard1911; Rybczynski Reference Rybczynski1911) and $C_D=24/Re_b$ (Stokes Reference Stokes1851). The former and latter models are valid in the low-Reynolds-number limit for bubbles in a pure liquid and for bubbles in a liquid contaminated by surfactants, respectively. In all cases, the bifurcations of whole- and multi-layered critical modes are found to be subcritical and supercritical, respectively.

3.2. Energy budget

Nakamura et al. (Reference Nakamura, Yoshikawa, Tasaka and Murai2020) showed essential differences between the instabilities of thermal and bubble-induced convections by an energy budget analysis. Both convections result from heterogeneous distribution of buoyancy sources, i.e. the temperature and the gas fraction in thermal and bubble-induced convections, respectively. However, the transport processes of these sources and resulting effects are distinct from each other. In thermal convection, the advection and diffusion of thermal energy bring about destabilising and stabilising effects, respectively. In bubble-induced convection, the lift and drag forces on bubbles tend to homogenise bubble distributions and, as a consequence, to stabilise the base state, while the inertial effect of liquid brings about destabilising effects.

We perform a similar energy budget analysis to consider the energy transfer from the base to perturbation flows and to reveal the driving mechanism of convective flows. Since (2.8b) and (2.8c) governing the liquid phase dynamics are analogous to the corresponding equations for thermal convection (Drazin & Reid Reference Drazin and Reid2010), the energy transfer to liquid perturbation flows is similar to that in the RB convection: the flows gain and lose energy due to buoyancy and viscous energy dissipation, respectively. The energy transfer to perturbation flows of the gas phase is more insightful. Taking the inner product of (2.8b) with $\boldsymbol {u}_G'$ and averaging the resulting equation over the whole fluid domain, we obtain the evolution equation of the kinetic energy of the gas phase:

(3.1)\begin{equation} \frac{\textrm{d}K_G}{\textrm{d}t}=W_I+W_D+W_L, \end{equation}

where the kinetic energy $K_G$ and the powers of the inertia $W_I$, of the drag $W_D$ and of the lift $W_L$ are defined by

(3.2a)\begin{equation} K_G = \left\langle \frac{u_G^{\prime 2} + w_G^{\prime 2}}{2} \right\rangle, \quad W_I = \left\langle w_I\right\rangle, \quad W_D = \left\langle w_D\right\rangle, \quad W_L = \left\langle w_L\right\rangle,\end{equation}

with

(3.2b) \begin{align} w_I&={-}\left[\bar{w}_G \frac{\partial}{\partial z} \frac{\|\boldsymbol{u}_G'\|^{2}}{2}+\frac{d \bar{w}_G}{\textrm{d}z}{w_G'}^{2}+(\boldsymbol{u}_G' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_G') \boldsymbol{\cdot} \boldsymbol{u}_G' \right] \nonumber\\ &\quad+3\left[\frac{\partial \boldsymbol{u}'}{\partial t} \boldsymbol{\cdot} \boldsymbol{u}_G' + (\boldsymbol{u}' \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}') \boldsymbol{\cdot} \boldsymbol{u}_G' \right], \end{align}

(3.2c)\begin{equation} w_D={-}\frac{18}{r_b^{2}}(\boldsymbol{u}_G'-\boldsymbol{u}')\boldsymbol{\cdot} \boldsymbol{u}_G' , \quad w_L=\bar{w}_G \zeta' u_G' + (u' w_G' -w' u_G')\zeta'. \end{equation}

The angle brackets denote the following integral operation: $\left \langle \, \bullet \, \right \rangle =({\alpha }/{2{\rm \pi} })\int _{-1}^{1} \int _0^{2{\rm \pi} /\alpha } \bullet {\textrm {d}\kern0.06em x}\, \textrm {d}z$. The integrands $w_I$, $w_D$ and $w_L$ are local densities of the inertia, drag and lift powers.

For both whole- and multi-layered convective flows, the drag is impeding ($W_D < 0$) convective flows (figure 3). The roles of the inertia and lift, however, depend on the mode type and can vary with the flow nonlinearity. In the case of the whole-layered flows (figure 3$a$), the inertia and lift forces are, respectively, driving ($W_I > 0$) and impeding ($W_L < 0$) convections close to the critical state as observed in the linear stability analysis (Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020). However, these effects relative to the drag force, i.e. the effects measured by $W_I/|W_D|$ and $W_L/|W_D|$, diminish and increase with the flow development, as shown in the inset of figure 3$(a)$. This implies that the lift plays an important role to drive convection at $Ra$ smaller than the critical value. In the case of the multi-layered modes, the powers of inertia and lift forces remain, respectively, positive and negative during the flow development along the bifurcation curve (figure 3$b$). Around the revolution of the curve, however, the effect of the lift relative to the drag decreases and switches finally to a driving one.

Figure 3.

Energy transfer from the base to perturbation flows of gas phase for (a) whole-layered modes and (b) multi-layered modes. The powers $W_I$, $W_D$ and $W_L$ of the inertia, drag and lift forces are shown in function of the bifurcation parameter $\delta$. These powers are normalised by twice the kinetic energy $K_G$ of gas perturbation flows. In panel (a), the powers $W_I$ and $W_L$ normalised by $|W_D|$ are also shown as inset for $\delta \in [-0.15, \, 0]$.