Self-similar solution for laminar bubbly flow evolving from a vertical plate

Abstract

The development of a bubble plume from a vertical gas-evolving electrode is driven by buoyancy and hydrodynamic bubble dispersion. This canonical fluid mechanics problem is relevant for both thermal and electrochemical processes. We adopt a mixture model formulation for the two-phase flow, considering variable density (beyond Boussinesq), viscosity and hydrodynamic bubble dispersion. Introducing a new change of coordinates, inspired by the Lees–Dorodnitsyn transformation, we obtain a new self-similar solution for the laminar boundary layer equations. The results predict a wall gas fraction and gas plume thickness that increase with height to the power of 1/5 before asymptotically reaching unity and scaling with height to the power 2/5, respectively. The vertical velocity scales with height to the power of 3/5. Our analysis shows that self-similarity is only possible if gas conservation is entirely formulated in terms of the gas specific volume instead of the gas fraction.

1. Introduction

Gas evolution from a vertical plate arises in boiling as well as various electrochemical processes, including water electrolysis for the production of green hydrogen (Le Bideau et al. 2020; Lee et al. 2022; Khalighi et al. 2023), the production of chlorine and chlorate (Hedenstedt et al. 2017) and the production of aluminium (Suzdaltsev, Nikolaev & Zaikov 2021). Buoyant forces set the electrolyte in an inexpensive, convective motion, which is advantageous for the mass transport of reactants and products and the removal of heat and bubbles (Haverkort 2024a). Similar to heat, bubbles diffuse away from high concentrations. Besides buoyancy, this is the dominant contribution to the bubble motion at low gas fractions (Dahlkild 2001; Schillings, Doche & Deseure 2015; Rajora & Haverkort 2023). While such two-phase systems resemble thermally driven natural convection, the analogy of gas fraction with temperature does not work exactly at higher gas fractions due to the increasing importance of other sources of bubble slip. Also, the usual Boussinesq approximation no longer holds due to the strong dependence of density and viscosity with void fraction. This makes reusing well-known self-similarity results for thermal convection (Ostrach 1953; Sparrow & Gregg 1956; Sparrow, Eichhorn & Gregg 1959) less accurate unless the gas fraction is very low. Instead, we need to employ a two-fluid formulation to account for the effects of variable density and viscosity (Ishii & Hibiki 2011). The mixture model is an affordable and common choice for this kind of problem, which has already produced good results in water electrolysis (Dahlkild 2001; Schillings et al. 2015, 2017). An alternative analogy is with particle-laden flows (Osiptsov 1981), which is only present for forced flows and diluted suspensions. Due to the nonlinearity of the Navier–Stokes equations, buoyant convective motion can develop instabilities – and eventually turbulence – along a flat plate (Osiptsov 1981; Boronin & Osiptsov 2008). However, in this work we limit ourselves to the study of the laminar regime.

The mathematical formulation of the mixture model resembles that of high-speed compressible flow, for which extensive literature on compressible boundary layers exists (Anderson 2006; Schlichting & Gersten 2017) for aerodynamic applications. Whereas gravity is typically disregarded in high-speed flows, it will be essential for gas evolution applications. Furthermore, whereas often a constant wall temperature is assumed, for gas evolution a constant gas flux is usually a more relevant boundary condition. Including buoyancy, taking a constant gas flux at the wall and considering variable physical properties constitute the main novelties of the analysis presented in this paper. They allow us to obtain new analytical solutions for the industrially important configuration of gas evolution from vertical walls or electrodes.

2. Formulation

We consider a steady-state two-phase flow of dispersed bubbles as depicted in figure 1, and use a two-fluid formulation (Ishii & Hibiki 2011) to describe the average flow field and phase distribution along the vertical plate. Assuming no acceleration between liquid ($l$) and gas ($g$) phases, typically valid for bubbly flows, we adopt the mixture model (Taivassalo & Kallio 1996; Ishii & Hibiki 2011). This model considers the gas–liquid mixture as a single fluid and includes additional closure relations for the unresolved flow features. The mixture density ($\rho$) is the volume-averaged density of the liquid and gas phases, which in terms of the gas volume fraction ($\varepsilon$) and liquid density ($\rho _{l}$) reads

(2.1)\begin{equation} \rho = \rho_{{l}} (1-\varepsilon), \end{equation}

where we neglected the contribution of the gas phase to the mixture density due to its typically much lower density (i.e. $\rho _{g} \ll \rho _{l}$). We then introduce the volume ($\boldsymbol {U}$)- and mass ($\boldsymbol {u}$)-averaged velocities

(2.2)$$\begin{gather} \boldsymbol{U} \equiv \boldsymbol{U}_{{l}} + \boldsymbol{U}_{{g}}, \end{gather}$$

(2.3)$$\begin{gather}\rho \boldsymbol{u} \equiv \rho_{{l}} \boldsymbol{U}_{{l}} + \rho_{{g}} \boldsymbol{U}_{{g}}, \end{gather}$$

where $\boldsymbol {U}_{{l}} \equiv (1-\varepsilon )\boldsymbol {u}_{{l}}$ and $\boldsymbol {U}_{{g}} \equiv \varepsilon \boldsymbol {u}_{{g}}$ are the superficial liquid and gas velocities expressed in terms of the interstitial velocities $\boldsymbol {u}_{{l}}$ and $\boldsymbol {u}_{{g}}$. Note that introducing the mixture formulation amounts to averaging the governing equations over a sufficiently large space/time scale. However, the description of the multiphase flow system requires capturing subscale features (e.g. those related to the motion of bubbles with respect to the mixture) which need to be modelled instead via so-called closure models. This includes the superficial slip velocity,

(2.4)\begin{equation} \boldsymbol{U}_{{slip}} \equiv \boldsymbol{U}_{{g}} - \varepsilon \boldsymbol{U} = \varepsilon (\boldsymbol{u}_{{g}} - \boldsymbol{U} ), \end{equation}

which is sometimes referred to as the drift flux or slip flux. The slip velocity is the difference between the real flux of gas (i.e. $\varepsilon \boldsymbol {u}_{{g}}$) and the flux predicted by the mixture (i.e. $\varepsilon \boldsymbol {U}$).

Figure 1. Schematic of the set-up for a gas-evolving vertical plate. Typical profiles of vertical velocity $w$ and the gas-to-liquid volume ratio $\theta \propto \varepsilon /(1-\varepsilon )$ at an arbitrary location.

Since neither $\boldsymbol {u}_{{g}}$ nor $\boldsymbol {u}_{{l}}$ is resolved, central to the mixture model are the semi-empirical relations that describe the slip velocities. Slip velocities result from the balance between drag and bubble forces, like lift, buoyancy or bubble–bubble interactions. Additional terms, such as bubble coalescence, break-up or turbulent interaction, are neglected in this study. In defining the slip velocities, we assume that the vertical buoyant rise velocity of bubbles relative to the liquid can be neglected with respect to the liquid velocity itself. In the horizontal direction, we neglect the lift forces and consider only the bubble–bubble interactions that lead to hydrodynamic dispersion:

(2.5)\begin{equation} \boldsymbol{U}_{{slip}}=-D_{eff}\boldsymbol{\nabla} \varepsilon, \end{equation}

where we introduce a new, effective, hydrodynamic dispersion coefficient:

(2.6)\begin{equation} D_{{eff}} = \frac{D_{{b}}}{(1-\varepsilon)^2},\end{equation}

which includes a $1/( 1-\varepsilon )^{2}$ factor to approximately account for the additional repulsive effects that arise as the gas fraction increases. Details concerning the derivation of this particular form can be found in Appendix C. This can be seen as one of the mechanisms keeping volume fractions below $1$. For particle suspensions, this is often modelled by adding a ‘solid pressure’ (Johnson & Jackson 1987), which mathematically has a similar effect.

The magnitude of the diffusivity is obtained by analogy with solid particle suspensions, which was found to be approximately given by the particle radius times the Stokes settling velocity (Ham & Homsy 1988; Nicolai et al. 1995), so

(2.7)\begin{equation} D_{{b}} = \frac{d_{{b}} w_{St}}{2},\end{equation}

with $d_{{b}}$ the bubble diameter and $w_{St} = { g d_{{b}}^2}/{18\nu _{{l}}}$ the Stokes rise velocity in terms of the kinematic liquid viscosity $\nu _{{l}}=\mu _{{l}}/\rho _{{l}}$. Note that we replace here the bubble size distribution with a single representative bubble diameter.

While this form of the slip velocity expression is arguably the least validated assumption of our model, it does reduce to the correct empirically validated limit at low gas fractions and allows us to describe the essence of solid pressure effects in a simplified way. In addition, only this particular form allows for a self-similarity solution of the boundary layer equations, as is demonstrated in Appendix A. As shown in § 5, the predictions following from this assumption are in good agreement with experiments.

For flows developing a laminar boundary layer along a vertical plate, we modify the compressible version of Prandtl's boundary layer theory (Anderson 2006) to include gravity:

(2.8)$$\begin{gather} \partial_x ( \rho u ) + \partial_z ( \rho w )= 0 , \end{gather}$$

(2.9)$$\begin{gather}\rho u w_x + \rho w w_z = \partial_x (\mu w_x ) + ( \rho_{l} - \rho ) g , \end{gather}$$

(2.10)$$\begin{gather}\rho u \partial_x \frac{1}{\rho} + \rho w \partial_z \frac{1}{\rho} = \partial_x \left(\rho_{l}\frac{D_{{b}} }{1-\varepsilon} \partial_x \frac{1}{\rho} \right), \end{gather}$$

where $u$ and $w$ are the components of $\boldsymbol {u}$ in the horizontal $x$ direction and vertical $z$ direction, respectively. Far away from the wall, where $(\rho _{{l}} - \rho )g=0$, a hydrostatic pressure gradient is assumed to cancel the gravitation force. Close to the wall, this term is positive, due to the upwards force that buoyancy exerts on the mixture. The gas transport equation (2.10) follows from taking the divergence of $u = U - U_{slip}/ (1-\varepsilon)$, $\boldsymbol {\nabla }\boldsymbol {\cdot}\boldsymbol {U}=0$ (since gas and liquid are conserved, i.e. $\boldsymbol {\nabla }\boldsymbol {\cdot}\boldsymbol {U}_{{g}}=\boldsymbol {\nabla }\boldsymbol {\cdot}\boldsymbol {U}_{{l}}=0$), inserting (2.1) and (2.5) and finally invoking $1/( 1-\varepsilon )^2 \partial _x \varepsilon = \partial _x ( 1 / ( 1-\varepsilon ) )$. Note that additional $\rho$ and $1/\rho$ terms are included in the convective term to render the final form in terms of the mass-averaged mixture velocity $(u,w)$ and the specific volume $1/\rho$.

The resulting equations (2.8)–(2.10) resemble those of high-speed aerodynamics, inspiring our solution approach. It is worth mentioning that the last equation corresponds to the conservation of the gas phase expressed in terms of specific volume ($1/\rho$).

We assume the following rheological relation for mixture viscosity (Ishii & Zuber 1979):

(2.11)\begin{equation} \mu = \frac{\mu_{{l}}}{1-\varepsilon}, \end{equation}

and define the following boundary conditions:

(2.12a,b)$$\begin{align} u\rvert_{x=0} = w\rvert_{x=0} &= 0 , \lim_{x\to\infty}w = 0 , \end{align}$$

(2.13a,b)$$\begin{align}\left.-\frac{D_{{b}}}{( 1-\varepsilon )^2} \frac{\partial \varepsilon}{\partial x}\right\vert_{x=0} &= U_{{w}} (1-\varepsilon_{{w}}) ,\quad \lim_{x\to \infty} \varepsilon = 0, \end{align}$$

which completes the formulation of the problem. The Neumann boundary condition for the gas fraction follows from (2.4) and (2.5), which give $- ({D_{{b}}}/{( 1-\varepsilon )^{2}}) \boldsymbol {\nabla } \varepsilon = \boldsymbol {U}_{{g}}-\boldsymbol {U}\varepsilon$, and assuming $\boldsymbol {U}_{{l}}\approx 0$ owing to the very low specific volume of the liquid with respect to the gas phase. In typical electrochemical systems, the superficial gas velocity at the wall $U_{{w}} = j\mathcal {V}_{{m}}/n{F}$ can be related to the current density $j$, the number of gas molecules $n$ produced per electron converted in the reaction, Faraday's constant ${F}$ and the molar volume of the gas $\mathcal {V}_{m}$ – typically given by the ideal gas law as $\mathcal {V}_{{m}} = {R}T/p$, where ${R}$ is the ideal gas constant, $T$ the absolute temperature and $p$ the pressure.

3. Methodology

Next, we summarize the steps leading to the development of the self-similarity solution. The procedure is similar to that employed in obtaining a self-similar solution for high-speed compressible flows (Anderson 2006), with variations to include the dependency on gravity.

Stream function formulation. We first introduce a stream function $\varPsi$ for the mass flux:

(3.1a,b)\begin{equation} \rho u =-\varPsi_z , \quad \rho w = \varPsi_x, \end{equation}

which makes (2.8) automatically satisfied. We rewrite the other governing equations as

(3.2)$$\begin{gather} - \varPsi_z \partial_x \frac{\varPsi_x}{\rho} + \varPsi_x \partial_z \frac{\varPsi_x}{\rho} = \partial_x \left( \frac{\rho_{{l}}\mu_{{l}}}{\rho} \partial_x \frac{\varPsi_x}{\rho} \right) + ( \rho_{{l}}-\rho )g, \end{gather}$$

(3.3)$$\begin{gather}- \varPsi_z \partial_x \frac{1}{\rho} + \varPsi_x \partial_z \frac{1}{\rho} = \partial_x \left(D_{{b}}\frac{ \rho_{l}^2}{\rho} \partial_x \frac{1}{\rho} \right), \end{gather}$$

where subscripts $x$ and $z$ denote partial derivatives.

Dorodnitsyn transformation. Now, we introduce the following self-similarity variable $\eta$:

(3.4)\begin{equation} \eta = \frac{Ar_{*z}^{1/5}}{z}\int_0^x \frac{\rho}{\rho_{{l}}}\, {{\rm d}\kern0.06em x}', \end{equation}

where we introduce the modified Archimedes number $Ar_{*z} = gU_{w}\,z^4/\nu _{l}^2 D_{b}$. This self-similarity transformation can be seen as combining the transformation proposed by Sparrow & Gregg (1956) for Boussinesq flows developing along a vertical flat plate subject to constant heat flux with the original Dorodnitsyn transformation (Anderson 2006) for compressible flows. Note that in this particular type of transformation, the new coordinate $\eta$ explicitly depends on the vertical coordinate $z$, while the dependency on the horizontal $x$ coordinate is introduced via the density-weighted integral. The intent of such transformation is to simplify the calculation of $w$ in Prandtl's boundary layer equations. In particular,

(3.5)\begin{equation} w = \frac{\varPsi_x}{\rho} = \frac{\eta_x}{\rho}\varPsi_\eta = \frac{Ar_{*z}^{1/5}}{z \rho_{l}} \varPsi_\eta,\end{equation}

which removes the explicit dependency on $\rho$ of the terms involving $w$, arguably simplifying the treatment of the boundary layer equations despite the variable density. Step-by-step details on the application of such transformation can be found in Appendix A.

Function transformation. Next, we adopt the following transformation for the unknowns $\varPsi$ and $1/(1-\varepsilon )$ in terms of the self-similar functions $f(\eta )$ and $\theta (\eta )$:

(3.6)$$\begin{gather} \varPsi = \mu_{{l}} Ar_{*z}^{1/5} f(\eta), \end{gather}$$

(3.7)$$\begin{gather}\frac{1}{1-\varepsilon} = 1 +Pe_{{b}z} Ar_{*z}^{-1/5} \theta(\eta), \end{gather}$$

where we introduce the Péclet number $Pe_{{b}z} = U_{w} z / D_{b}$. These transform equations (3.2) and (3.3) into the system

(3.8)$$\begin{gather} f''' = \frac{3}{5} \,f'^2 - \frac{4}{5}\, f\, f'' - \theta, \end{gather}$$

(3.9)$$\begin{gather}\theta'' = \frac{Pr_{{b}}}{5} (\,f'\theta - 4\, f\theta'), \end{gather}$$

where the bubble Prandtl number $Pr_{{b}} = \nu _{{l}} / D_{{b}}$. The system is subject to the boundary conditions

(3.10)$$\begin{align} f\rvert_{\eta=0} = f'\rvert_{\eta=0} &=0 \quad \lim_{\eta\to \infty}f' =0, \end{align}$$

(3.11)$$\begin{align}\theta'\rvert_{\eta=0} &=-1 \lim_{\eta\to \infty}\theta =0, \end{align}$$

where a prime ($'$) denotes a derivative with respect to $\eta$ to shorten the expressions. These transformations now allow us to derive a similarity solution. The system of (3.9) and (3.10) is the equivalent of the classical result for Boussinesq natural convection with constant heat flux (Sparrow & Gregg 1956). However, the most remarkable change is the shape of the buoyant scalar (i.e. $\varepsilon$). In particular, (2.10) gives

(3.12)\begin{equation} \varepsilon = \frac{Pe_{{b}z} Ar_{*z}^{-1/5}\theta}{1+Pe_{{b}z} Ar_{*z}^{-1/5}\theta}, \end{equation}

whereas in the thermal convection case (Sparrow & Gregg 1956) $\varepsilon = Pe_{{b}z} Ar_{*z}^{-1/5}\theta$, which agrees with (3.12) when $Pe_{{b}z} Ar_{*z}^{-1/5}\theta {\ll 1}$. This is expected since for very low gas fractions, we should recover the Boussinesq hypothesis. The denominator in (3.12) ensures that the gas fraction does not exceed unity, a phenomenon in bubbly flow with no equivalent in thermal natural convection.

The vertical velocity field can be recovered from (3.5) and (3.6) as

(3.13)\begin{equation} w = \frac{\nu_{l}}{z} Ar_{*z}^{2/5} f', \end{equation}

which gives a scaling with height proportional to $z^{3/5}$, equal to the thermal convection solution (Sparrow & Gregg 1956).

The mass flow rate per unit width [${\rm kg}\ {\rm m}^{-1}\ {\rm s}^{-1}$] is then equal to

(3.14)\begin{equation} \int_0^\infty \rho w \,{{\rm d}\kern0.06em x} = \int_0^\infty \varPsi_x \,{{\rm d}\kern0.06em x} = \varPsi_\infty = \mu_{{l}}Ar_{*z}^{1/5}f_\infty, \end{equation}

which scales with $z^{4/5}$, accounting for the suction in the horizontal direction from outside the plume.

The wall shear stress $\tau _{w}$ [Pa] is given by

(3.15)\begin{equation} \tau_{w} \equiv \left. { \mu \frac{\partial w}{\partial x}}\right\vert_{w} = \mu_{l} \frac{\nu_{l}}{z^2} Ar_{*z}^{3/5} f''_{w}, \end{equation}

where we used (2.11) for the mixture viscosity $\mu = {\mu _{{l}}}/({1- \varepsilon })$ as well as (2.1), (3.4) and (3.13).

The wall shear stress shows the same scaling with height and current density irrespective of the local gas fraction, as in the Boussinesq hypothesis case (Sparrow & Gregg 1956). This is due to the variations of $\mu$ and $\rho$ with $\varepsilon$ cancelling each other. However, the wall shear rate $w'_{{w}} = {\tau _{{w}} }/{ \mu } = ({\tau _{{w}} }/{ \mu _{{l}}}) ( 1- \varepsilon _{{w}} )$ [s$^{-1}$] does show a transition between $z^{2/5}$ and $z^{3/5}$. Using (3.12) this gives

(3.16)\begin{equation} w'_{{w}} = \frac{\nu_{l}}{z^2} \frac{f''_{w} Ar_{*z}^{3/5}}{1+Pe_{{b}z} Ar_{*z}^{-1/5}\theta_{{w}}}. \end{equation}

Finally, we may define the dimensionless gas plume thickness as ${\delta _{{g}}}/{z} = - ({z^{-1}}/\eta_x) ({\theta _{{w}} }/{\theta '_{{w}} })$. It follows from the boundary conditions in (3.11) and (3.4) and (2.1) that

(3.17)\begin{equation} \frac{\delta_{{g}}}{z} = Ar_{*z}^{-1/5} ( 1 + Pe_{{b}z}Ar_{*z}^{-1/5} \theta_{{w}} )\theta_{{w}}, \end{equation}

which shows a transitional scaling between $z^{-4/5}$ and $z^{-3/5}$. Note, however, that as $z \to \infty$, the model loses accuracy as turbulence may set in.

4. Results

The final system of (3.8)–(3.11) can be solved by numerical integration for different values of $Pr_{{b}}$, which we choose corresponding to a wide range of bubble diameters.

The numerical method consists of a classical shooting algorithm, which adjusts the value for $\theta$ at $\eta =0$ iteratively until the boundary conditions are satisfied. The domain of integration ranges from $\eta =0$ to $\eta = 50$. This range was chosen such that parameters in the vicinity of the wall do not change substantially upon a further increase. For $f_\infty$, convergence was achieved for a substantially larger domain up to $\eta =800$. Results for $f'$ and $\theta$ are shown in figure 2. Through (3.12) and (3.13) these dimensionless quantities are related to the gas fraction $\varepsilon$ and the vertical velocity $w$, respectively.

Figure 2. The numerical solution of (3.8) and (3.9) for $f'$ (a) and $\theta$ (b) near the wall for several $Pr_{{b}}$ numbers as indicated in the figure. The computational domain is $\eta \in [0,50]$, but for visualization purposes, a shorter interval is shown.

From these profiles, we can obtain relevant parameters such as wall gas fraction $\varepsilon _{w}$ and wall shear stress $\tau _{w}$, which depend on $\theta _{w}$ and $f''_{w}$; and also on vertical mass flow rate, which depends on $f_\infty$. Values of $\theta _{w}$, $f''_{w}$ and $f_\infty$ are reported in figure 3 for different $Pr_{{b}}$ numbers. All profiles show a transition at $Pr_{{b}}\approx 1$, corresponding to the switch in the predominance of the momentum diffusivity $\nu$ over the bubble diffusivity $D_{{b}}$. For convenience, figure 3 has as an additional axis showing typical bubble sizes assuming Stokes drag and an aqueous-like viscosity. Results show different asymptotic behaviours for small bubbles ($d_{{b}}<100\ {\rm \mu}\textrm {m}$, typical for electrolytic bubbles) and larger bubbles ($d_{{b}}>100\ \mathrm {\mu }\textrm {m}$, typical for bubbles due to boiling). A power-law scaling is obtained for the two regions, shown in the figure, resulting in approximate asymptotic relations for both parameters:

Figure 3. The wall values $\theta _{w}$ and $f''_{{w}}$ and $f_{\infty }$ as a function of the bubble Prandtl number $Pr_{{b}} = \nu / D_{{b}}$. The upper $x$ axis shows the corresponding bubble diameter using (2.7), assuming typical values of $g=9.81\ {\rm m}\ {\rm s}^{-2}$ and $\nu _{{l}}=10^{-6}\ {\rm m}^2\ {\rm s}^{-1}$. We note that for bubbles with $d_{{b}} \gtrsim 100\ \mathrm {\mu }$m, $w_{{St}}$ used in (2.7) no longer corresponds to the buoyant slip velocity as Stokes drag becomes invalid.

(4.1)$$\begin{gather} \theta_{w} \approx \begin{cases} 1.63Pr_{{b}}^{-1/5} & Pr_{{b}} \gtrsim 1 \\ 1.63 Pr_{{b}}^{-1/3} & Pr_{{b}} \lesssim 1, \end{cases} \end{gather}$$

(4.2)$$\begin{gather}f''_{w} \approx \begin{cases} 1.50Pr_{{b}}^{-2/5} & Pr_{{b}} \gtrsim 1 \\ 1.50 Pr_{{b}}^{-1/3} & Pr_{{b}} \lesssim 1, \end{cases} \end{gather}$$

(4.3)$$\begin{gather}f_\infty \approx \begin{cases} 1.40Pr_{{b}}^{-3/5} & Pr_{{b}} \gtrsim 1 \\ 1.25 Pr_{{b}}^{-3/10} & Pr_{{b}} \lesssim 1. \end{cases} \end{gather}$$

Approximations combining both asymptotic limits can be found in Appendix A. These values can be used to approximate wall gas fraction $\varepsilon _{w}$, wall shear stress $\tau _{w}$ and wall shear rate $w'_{{w}}$ through (3.12), (3.15) and (3.16), respectively, to give

(4.4)$$\begin{gather} \frac{ \varepsilon_{w} }{1 - \varepsilon_{w}} \approx 1.63 \times \begin{cases} \left( \dfrac{ \nu_{l} U_{{w}}^4\, z }{g D_{b}^3 } \right)^{1/5} & d_{{b}}\lesssim 100 \ \mathrm{\mu}\textrm{m} \\ \left( \dfrac{\nu_{l}^{1/3} U_{{w}}^4\, z}{g D_{b}^{7/3} } \right)^{1/5} & d_{{b}}\gtrsim 100\ \mathrm{\mu}\textrm{m}, \end{cases} \end{gather}$$

(4.5)$$\begin{gather}\tau_{w} \approx 1.5 \rho_{{l}} \times \begin{cases} \left( \dfrac{g^3 \nu_{l}^2 U_{w}^3\, z^2 }{ D_{{b}} }\right)^{{1}/{5}} & d_{{b}}\lesssim100\ \mathrm{\mu}\textrm{m} \\ ( g^3 (\nu_{l}^{7}/D_{{b}}^{4} )^{{1}/{3}} U_{w}^3\, z^2 )^{{1}/{5}} & d_{{b}}\gtrsim100\ \mathrm{\mu}\textrm{m} \\ \end{cases} \end{gather}$$

and

(4.6)\begin{equation} w'_{w} \approx 1.5\times \begin{cases} \dfrac{ \left( \dfrac{g^3 U_{w}^3\, z^2 }{ D_{{b}} \nu_{l}^3 }\right)^{1/5} }{1+1.63 \left( \dfrac{ \nu_{l} U_{{w}}^4\, z }{g D_{b}^3 } \right)^{1/5} } & d_{{b}}\lesssim100\ \mathrm{\mu}\textrm{m} \\ \dfrac{\left( \dfrac{ g^3 U_{w}^3\, z^2 }{ (\nu_{l}^{8} D_{{b}}^{4} )^{{1}/{3}} } \right)^{1/5}}{1+1.63 \left( \dfrac{\nu_{l}^{1/3} U_{{w}}^4\, z }{g D_{b}^{7/3} } \right)^{1/5} } & d_{{b}}\gtrsim100 \ \mathrm{\mu}\textrm{m}. \end{cases} \end{equation}

This last relation shows a transition from a proportionality with $U_{{w}}^{3/5}z^{2/5}$ at low wall gas fractions to $U_{{w}}^{-1/5}z^{1/5}$ at wall gas fractions close to one. Interestingly, the viscosity increases faster with increasing $U_{{w}}$ than the wall shear stress, resulting in a decreasing wall shear rate with increasing gas flux in this regime. This limiting result strongly depends on the exact relation between effective viscosity and gas fraction, i.e. (2.11), and the assumed maximum gas fraction of unity.

The mass flux per unit width, from (3.14), reads

(4.7)\begin{equation} \int_0^\infty \rho w \,{{\rm d}\kern0.06em x} = \begin{cases} 1.4\rho_{{l}} ( g D_{b}^2 U_{w}\,z^4 )^{1/5} & d_{{b}}\lesssim100 \ \mathrm{\mu}\textrm{m} \\ 1.25\rho_{{l}} ( g \nu_{l}^{3/2} D_{b}^{1/2} U_{w}\,z^4)^{1/5} & d_{{b}}\gtrsim 100 \ \mathrm{\mu}\textrm{m}. \end{cases} \end{equation}

This shows an almost linear increase with $z$ and a much weaker dependence on the wall gas flux $U_{{w}}$. Rajora & Haverkort (2023) found a linear dependence on both $z$ and $U_{{w}}$, for a plume in the shape of a step function. When $Pe_{{b}z} Ar_{*z}^{-1/5} \gg 1$ we have $\varepsilon _{{w}}\approx 1$ and the gas plume described by (3.12) similarly becomes a smoothed step function. However, in this limit, there is no liquid inside the plateau region of the gas plume and the analysis of Rajora & Haverkort (2023) fails in this case.

Finally, using (3.17), (4.1) and (4.4) the gas plume thickness becomes

(4.8)\begin{equation} \delta_{{g}} \approx 1.63 \left(\frac{z}{ g U_{w} } \right)^{{1}/{5}} \times \begin{cases} \nu_{l}^{{1}/{5}} D_{b}^{{2}/{5}}\left( 1+\left( \frac{ \nu_{l} U_{{w}}^4\, z }{g D_{b}^3 } \right)^{1/5} \right) & d_{{b}}\lesssim100\ \mathrm{\mu}\textrm{m} \\ \nu_{l}^{{1}/{15}} D_{b}^{{8}/{15}}\left( 1+ \left( \frac{\nu_{l}^{1/3} U_{{w}}^4\, z }{g D_{b}^{7/3} } \right)^{1/5} \right) & d_{{b}}\gtrsim100 \ \mathrm{\mu}\textrm{m}. \end{cases} \end{equation}

This shows a transition from a proportionality with $(z / U_{{w}} )^{1/5}$ at low gas fraction to a proportionality with $U_{{w}}^{3/5}z^{2/5}$ at high wall gas fractions. The latter positive dependence on gas flux is in agreement with most experimental findings in which the plume thickness increases with increasing current density. A similar transition was found by Rajora & Haverkort (2023) using approximate methods, where at high gas fractions, depending on $Pr_{{b}}$, a proportionality of between $( U_{{w}}z)^{1/3}$ and $U_{{w}}^{0.73}z^{0.43}$ was obtained. The result obtained here is virtually exact, valid for all wall gas fractions, and arguably more elegant.

5. Validation

Here we compare our results with experimental data for gas-evolving electrodes. Two different sets of experiments are considered: wall shear stress measurements of Hiraoka et al. (1986) and gas plume thickness data of Fukunaka et al. (1989). Both datasets concern hydrogen-evolving electrodes.

The experiments of Hiraoka et al. (1986) report direct measurements of the average wall shear stress. Their set-up consisted of a planar vertical working electrode immersed in a KOH solution between two counter electrodes. The measured weight of the working electrode was recorded against the current density and subtracted from the weight that was measured when no gas evolution was present. This measurement allowed to determine $\tau _{w}$.

A comparison with our model predictions is shown in figure 4. We averaged wall shear stress for (3.15) as

(5.1)\begin{equation} \langle {\tau_w} \rangle = \frac{1}{h}\int_0^h \tau_w\, {\rm d}z = \frac{5}{7} \tau_w \vert_h. \end{equation}

Figure 4. Experimental values of the wall shear stress $\tau _{{w}}$ as a function of the gas-evolving velocity $U_{w}$. The corresponding current density for H$_2$ evolution is reported on the top $x$ axis. Data taken from Fukunaka et al. (1989). Bubble diameter is set to $d_{{b}} = 200\ \mathrm {\mu }\textrm {m}$.

The average wall shear stress at the electrode surface for electrodes of different heights in the range of $h \in [4\unicode{x2013}10]\ \textrm {cm}$ is shown. The predicted results are overall in good agreement with the experimental data except for very low gas-evolving velocities (i.e. low current densities).

The gas plume thickness is compared with the experimental data of Fukunaka et al. (1989). Their experimental set-up consists of a segmented working electrode immersed in a transparent beaker with 0.05 M CuSO$_4$–1.85 M H$_2$SO$_4$ solution and a single counter electrode. The gas plume thickness was measured from pictures taken from the side of the beaker. The bubble detachment diameter was found to be in the range $d_{{b}} \in [60\unicode{x2013}200]\ \mathrm {\mu }\textrm {m}$.

The experimental data of the gas plume thickness are compared with our model predictions in figure 5. The bubble diameter, $d_{b}$, was used as a fitting parameter, which provided results within the observed range. The obtained increase in bubble size with increasing current density is in good agreement with the data reported for H$_2$SO$_4$ by Janssen & Hoogland (1973). Figure 5 includes a comparison of the present theory with the Boussinesq hypotheses using a corresponding bubble diameter of $d_{{b,B}}$ (see figure 5). Our result shows a slope that is marginally closer to that of the data, but both results can explain the data reasonably well and with reasonable bubble diameters. The reason is that not even a decade of different heights is included in the data, and (4.8) depends rather weakly on height.

Figure 5. Experimental values of the gas plume thickness $\delta _{g}$ as a function of height h for the new model (solid line) and with the classical Boussinesq assumption (dashed line). Data taken from Hiraoka et al. (1986).

A much larger difference between the predictions of the present model and the Boussinesq approximation was found for the mass transfer coefficient, as recently reported in Valle & Haverkort (2024) and Haverkort (2024b). The reason is that the mass transfer coefficient is strongly dependent on the shear at the wall $w'$, which shows a much more acute difference. Here, the predictions from the present model are in much better agreement with the data than those of the Boussinesq approximation.

6. Conclusion

We have developed a new self-similarity solution for laminar bubbly flow evolving from a vertical plate, which considers variable density, viscosity and hydrodynamic dispersion, all depending on the local gas fraction. Results for both small and large bubbles produce simple power-law relations for the wall shear stress and the gas fraction at the wall, which are critical parameters in process technology. The wall strain rate and gas plume thickness show a transition between low to moderate gas fractions and gas fractions near one.

The new results go beyond the classical Boussinesq theory, which is only valid for very low gas fractions. Owing to the strong dependence of density and viscosity with gas fraction, the new results expand the validity of this analysis to general void fractions, as far as the flow remains laminar.

Our theoretical analysis reveals that a self-similarity solution is only possible if the bubble diffusivity is proportional to $(1-\varepsilon )^{-2}$, which is equivalent to formulating the gas fraction convection–diffusion equation in terms of specific volumes and using a constant diffusivity. While this is a convenient assumption that allows for a neat analytical treatment, it can also be physically motivated to approximately model increased diffusion at high gas fractions, preventing the system from reaching non-physical gas fractions $\varepsilon > 1$.

By introducing a novel change of coordinates and variables we observe that the solutions show an asymptotic transition between the classical Boussinesq result for relatively low gas flow rates and short heights and a new asymptotic solution for which the wall gas fraction tends to one.

Comparison with experimental results shows good agreement with existing experimental literature data for wall shear stress and bubble plume thickness. Considering the simplified nature of the model and the number of assumptions made, the agreement in the observed trends is highly satisfactory.

Further work in this area requires an improvement of the closure models for hydrodynamic dispersion of bubbles. An immediate step would be to perform a linear stability analysis of the base flow presented above.