1. Introduction
Liquid flows with dispersed gas phases have been investigated in various fields of scientific research to understand their complex multi-scale flow structures and to develop novel applications of energy and environmental systems (Lohse Reference Lohse2018). Presence of dispersed bubbles alters turbulent flow structures and induces large-scale liquid circulations (Risso Reference Risso2018). This pertains to industrial applications of, for example, bubble column reactors (Sommerfeld Reference Sommerfeld2004) used as efficient gas–liquid contactors for enhancement of mass and heat transfer, in which bubbles are injected into a liquid layer from the bottom. The injected bubbles rise in liquid due to Archimedean buoyancy forces, exchanging mass and heat at gas–liquid interfaces (Mudde Reference Mudde2005). The two-phase flows in reactors show a variety of flow patterns depending on the bubble size, the gas volume fraction, and the geometric shapes of reactors (Mudde Reference Mudde2005; Risso Reference Risso2018). When the gas injection flux is small, bubbles ascend with a uniform distribution (Mudde Reference Mudde2005; Ruzicka Reference Ruzicka2013). At large gas fluxes, the bubble distribution becomes heterogeneous and convective motions occur (Mudde Reference Mudde2005).
The generation of convection in two-phase flow systems with a bubble-rich bottom layer has been modelled as a gravitational instability provoked in an inverse density stratification (Iga & Kimura Reference Iga and Kimura2007). Ruzicka & Thomas (Reference Ruzicka and Thomas2003) elucidated the analogy to the Rayleigh–Bénard (RB) convection. The bubble-induced convection is thus often characterised in the literature by the Rayleigh and Prandtl numbers, 
$Ra$ and 
$Pr$. They are identical to those used for thermal convections but with the bubble residence time 
$T_{R}=d/V_{\infty }$ replacing the thermal diffusion time: 
$Ra= T_{\nu } T_{R} / T_{B}^{2}$, 
$Pr = T_{R} / T_{\nu }$, where 
$d$ is the height of a system, 
$V_{\infty }$ is the terminal velocity of a rising bubble, 
$T_{\nu }=d^{2}/\nu$ is the momentum diffusion time with the kinematic viscosity 
$\nu$ of the liquid, and 
$T_{B} = \sqrt {d/\varepsilon g}$ is the buoyancy time (
$\varepsilon$, the gas volume fraction; 
$g$, gravitational acceleration). Iga & Kimura (Reference Iga and Kimura2007) showed the formation of steady convection rolls at 
$Ra=2.5 \times 10^{5}$ in an experiment with bubbles of approximately 
$10 \ \mathrm {\mu }\textrm {m}$ in diameter. A direct numerical simulation (DNS) based on a point bubble model by Climent & Magnaudet (Reference Climent and Magnaudet1999) showed that steady convection rolls developed for 
$Ra \sim 2.0 \times 10^{5}$ from an initial two-layer state where a bubble-rich fluid layer underlies a pure liquid layer. These works both mentioned that the selection of the horizontal size of convection rolls is not unique and the size changes depending on randomly determined initial conditions.
Recently, we performed a linear stability analysis of bubble-induced convection (Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020), assuming an inverse density profile formed spontaneously from the accelerated ascending motion of bubbles (figure 1
$a$). The analysis is based on a two-fluid model, in which the added-mass, drag, shear-induced lift and buoyancy forces on bubbles are considered in the momentum exchanges between liquid and gas phases. In contrast to the RB convection, the predicted marginal stability curves 
$Ra=Ra_m(k)$, where k is the wavenumber of perturbation flows, have two local minima and their relationship varies as function of 
$Pr$. Convection rolls of the height of the fluid layer (whole-layered mode, see figure 1
$c$-i) are critical for 
$Pr$ smaller than a threshold 
$Pr^{*}$, while superposed convection rolls (multi-layered mode, see figure 1
$c$-ii) are critical otherwise (Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020). Although the values of 
$Pr^{*}$ and 
$Ra_m$ depend on the radius and the injection velocity of bubbles and on the model of drag coefficient, the qualitative behaviour of the marginal curves and the resulting mode selection remain the same.
Figure 1. (
$a$) Geometric configuration of the problem, (
$b$) profiles of bubble rising velocity 
$\bar {w}_G(z)$ and gas volume fraction 
$\bar {\varepsilon }_G(z)$ in the base state for different values of 
$Pr$, and (
$c$) eigenvectors at critical conditions for (i) a whole-layered mode and (ii) a multi-layered mode. The dimensionless bubble injection velocity and bubble radius are set as 
$w_{G, 0} = 1000$ and 
$r_b = 0.01$, respectively. The perturbation fields of liquid velocity 
$\boldsymbol {u}'$ are shown by vectors and streamlines. The perturbation fields of gas volume fraction 
$\varepsilon _G'$ are shown by colours. Bubble-rich and bubble-poor cells correspond to red and white zones.
The present work aims to reveal the nonlinear development of the whole- and multi-layered modes with a particular focus on their bifurcation behaviour. For this aim, nonlinear flows close to the critical states are computed to draw bifurcation diagrams. In thermal convection of non-Newtonian fluids, the bifurcations depend on the fluid rheology: it can be supercritical as in the RB convection of Newtonian fluids but also subcritical, when shear-thinning effects are significant (Bouteraa et al. Reference Bouteraa, Nouar, Plaut, Metivier and Klack2015). The presence of a dispersed solid phase can also alter the type of bifurcation in the RB convection through the variation of the effective viscosity of suspension resulting from the shear-diffusion of suspended particles (Kang, Yoshikawa & Mirbod Reference Kang, Yoshikawa and Mirbod2021). By the present investigation, we show that even in the absence of rheological effects the presence of dispersed phase can affect the bifurcation in two-phase flow systems. The considered flow system is isothermal and convection develops from an intrinsic gravitational instability of bubbly flows. Dilute bubbly flows are assumed so that the constitutive law of the continuous phase remains Newtonian. Following the introduction, we present a brief summary of a mathematical model used in our previous report (Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020) and describe a numerical method for nonlinear flow determination. The properties of the solution are then explained using bifurcation diagrams showing flow patterns. Finally, we discuss the roles of different forces on bubbles in the nonlinear development of convection.
2. Formulation of the problem
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
$$\begin{gather} \frac{\partial \varepsilon_G}{\partial t} +\boldsymbol{\nabla} \boldsymbol{\cdot} (\varepsilon_G \boldsymbol{u}_G)=0, \end{gather}$$
$$\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}$$
$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}$$
$$\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:
\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:
$$\begin{gather} \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} = 0, \quad w=0, \quad \text{at } z=1, \end{gather}$$
$$\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:
\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
\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
\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:
\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
$$\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}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}' = \boldsymbol{0}, \end{gather}$$
$$\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}$$
$$\begin{gather} \displaystyle \frac{\partial u'}{\partial z} + \frac{\partial w'}{\partial x} = w'=0, \quad \text{at }z=1, \end{gather}$$
$$\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:
\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$):
\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. Nonlinear properties of the flows
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:
\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
\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
\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}
\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]$.
4. Discussion
The sensitivity of nonlinear development of convections to the constitutive law of fluid is reported in the literature for thermal convections under different conditions: Parmentier (Reference Parmentier1978), Solomatov (Reference Solomatov2012) and Curbelo & Mancho (Reference Curbelo and Mancho2014) for Newtonian fluids with a temperature-dependent viscosity; Albaalbaki & Khayat (Reference Albaalbaki and Khayat2011), Benouared, Mamou & Messaoudene (Reference Benouared, Mamou and Messaoudene2014) and Bouteraa et al. (Reference Bouteraa, Nouar, Plaut, Metivier and Klack2015) for non-Newtonian fluids with shear-thinning effects; and Kang et al. (Reference Kang, Yoshikawa and Mirbod2021) for suspensions with an effective viscosity modelled by Krieger's law. The results reported in § 3.1 (figure 2
$a{,}b$) suggest that the bifurcation can also be altered by the law of transport of dispersed phase. As we have assumed dilute bubbly flows, the presence of bubbles has no rheological effect in contrast to the above-mentioned works on convections and affects the liquid dynamics only through the buoyancy force (2.1d). The essential difference of the considered system compared with thermal convections arises from the transport equation of bubbles (2.1b). The energy budget analysis in § 3.2 indicates a significant role of the lift force in sustaining convective flows. The power provided by the lift force behaves differently in different bifurcations (figure 3). The details of energy transfer process can be understood from the local distributions of different powers (figure 4). In the case of subcritical bifurcation, the complex coupling of the bubble transport law (2.1d) with developed convective flows intensifies the driving effect of the lift (
$w_L > 0$) close to the free surface and weakens the impeding effect of the lift (
$w_L<0$) in the lower part of convection rolls (figure 4
$a$). It results in the change of net effect of the lift force from impeding (
$W_L < 0$) to driving (
$W_L > 0$) convective flows. The distribution of the power of the inertia (
$w_I$) is similar to 
$w_L$ but with the opposite sign. In the case of supercritical bifurcation, the distribution of 
$w_L$ behaves differently (figure 4
$b$). Though the development of positive and negative zones of 
$w_L$ in the first convection layer at the bottom is similar to the subcritical case, negative zones of 
$w_L$ grow in the second convection layer (figure 4
$b$). The total effect of the lift, thus, remains impeding even after the revolution of the bifurcation curve.
Dispersed phase can thus change the bifurcation through the variations of different forces in its transport law. Laminar-turbulent transition in multi-phase flows would be affected by this mechanism. In experiments on the transition in particle-laden pipe flows, the transition from the laminar flow was observed to be supercritical at large values of particle volume fraction 
$\phi$ in contrast to the subcritical bifurcation of the Newtonian Hagen–Poiseuille flow (Hogendoorn & Poelma Reference Hogendoorn and Poelma2018; Agrawal, Choueiri & Hof Reference Agrawal, Choueiri and Hof2019). To model this transition, a thorough consideration on particle transport might be required in addition to a rheological consideration.
In the present study, only 2-D flows are considered, assuming a fluid layer of infinite horizontal extent. The presence of lateral boundaries can develop three dimensional (3-D) flow structures even at the first bifurcation point as observed in thermal convection (Dijkstra et al. Reference Dijkstra2014). In bubbly flow systems, furthermore, the physics at bubble scales, such as the path instability and bubble–bubble interactions (Risso Reference Risso2018), provides 3-D perturbations to flows and would induce transition to 3-D convections even in laterally non-confined systems. The bifurcation to 3-D flows is to be examined in a future work.
5. Conclusion
In the present paper, we have investigated the nonlinear development of bubble-induced convection from the critical states predicted in our previous work (Nakamura et al. Reference Nakamura, Yoshikawa, Tasaka and Murai2020). Bifurcations determined for the whole- and multi-layered modes, which are critical for the bubble Prandtl number 
$Pr$ smaller and larger than 
$Pr^{*}$, are subcritical and supercritical, respectively.
We have analysed the energy transfer from the base to perturbation flows of the gas phase. The analysis revealed that the bubble transports due to different forces play important roles in the nonlinear development of convection. In particular, the lift force changes its role from impeding to driving convective flows when the convection is intensified with decreasing control parameter. It is responsible for the observed subcritical behaviour of the convection. Local effects of different transport mechanisms showed that the variation in the role of the lift arises from the coupling of gas transport with nonlinear liquid flows.
The observed sensitivity of flow bifurcation to the transport of dispersed phase suggests that, in addition to the rheological effects, a thorough consideration on the transport would be required for modelling laminar-turbulent transitions in multi-phase flow systems.
Acknowledgements
We thank two anonymous reviewers for their careful reading of our manuscript and helpful comments. H.N.Y. would like to thank Professor M. Nagata for his introduction to bifurcation problems and helpful suggestions on the numerical method.
Funding
This work was supported by the Grant-in-Aid for JSPS Fellows from the JSPS KAKENHI grant no. 19J1064809.
Declaration of interests
The authors report no conflict of interest.
