3.1. Degassing-dominated PBE

Recalling (2.6), $S^+_f(a)$ is similar in magnitude to $S^-_f(a)$, so $\varLambda (a)\gg \varOmega (a)$ allows us to remove both the fragmentation terms from (3.1). This gives the degassing-dominated PBE, which, after rearranging terms, reads

(3.2)\begin{equation} I(a)= \varLambda(a)\,N(a) + \partial N(a) / \partial t {.} \end{equation}

We see that entrainment is balanced by two terms: degassing, and increasing bubble size distribution $\partial N(a) / \partial t$. We define the ratio of the first term to the second term:

(3.3)\begin{equation} \gamma(a) \equiv [\varLambda(a)\,N(a)]\big/[I(a)-\varLambda(a)\,N(a)] {.} \end{equation}

As a measure of this ratio across radii, we integrate (weighted by bubble volume $a^3$) the numerator and denominator to define $\varGamma \equiv Q_D/(Q_I-Q_D)$. This integral measure can be related directly to $\mathcal {D}$ as $\varGamma =\mathcal {D}/(1-\mathcal {D})$. For a given radius, $\gamma (a) \ll 1$ means that $\varLambda (a)\,N(a) \ll \partial N(a) / \partial t$, and (3.2) becomes $\partial N(a) / \partial t \approx I(a)$, giving $N(a)\sim t\,I(a)$. Because the degassing term grows linearly with $N(a)$, $\varLambda (a)\,N(a) \ll \partial N(a) / \partial t$ is true only for time scales less than order $1/\varLambda (a)$. Apart from these short time scales, we expect $\varLambda (a)\,N(a) > \partial N(a) / \partial t$, corresponding to $\gamma (a) > 1$, for degassing-dominated populations. In this case, the primary balance is between entrainment and degassing, and (3.2) becomes

(3.4)\begin{equation} N(a)\approx I(a)/\varLambda(a) {.} \end{equation}

Given an entrainment distribution $I(a)\propto a^{-10/3}$ (Yu et al. Reference Yu, Hendrickson and Yue2020), and assuming that degassing follows a power law $\varLambda (a)\propto a^{\alpha }$, the power-law slope of $N(a)$ is

(3.5)\begin{equation} \beta=-10/3-\alpha {.} \end{equation}

3.2. Power-law scaling of degassing in free-surface turbulence

Our objective is to find $\beta$ for degassing-dominated populations using (3.5). To predict $\alpha$, we derive a characteristic bubble depth $L_\varLambda$ and bubble rise velocity $U_\varLambda$ such that $\varLambda (a)\propto U_\varLambda /L_\varLambda$. The scaling of $U_\varLambda$ with bubble radius is a widely studied problem not unique to air-generating flows. A challenge is the characteristic depth $L_\varLambda$, which should describe some initial depth of an entrained bubble. For FST, turbulent eddies drive entrainment, and to derive $L_\varLambda$, we start by estimating how deep a single eddy will push a bubble. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) show that air-generating FST creates a region of nearly homogeneous isotropic turbulence (HIT) near the free surface. For HIT, an eddy of size $\ell$ has velocity $u_\ell \sim \varepsilon ^{1/3} \ell ^{1/3}$ and is coherent over $t_\ell \sim \varepsilon ^{-1/3} \ell ^{2/3}$. A bubble of radius $a$ has an (added) mass $m\sim \rho _w a^3$ and (at large ${{Re}}$) feels a force from the eddy $F\sim \rho _w a^2 u_\ell ^2$, giving an acceleration $\dot {w}\sim u_\ell ^2/a$. This moves the bubble to a depth $z\sim {t_\ell }^2 \dot {w}=\ell ^2/a$. To apply this single-eddy model to HIT where a range of eddy sizes are present, we assume that the energy-containing eddies $\ell \sim u_{rms}^3 \varepsilon ^{-1}$ are primarily responsible for $L_\varLambda$, giving $L_\varLambda \sim u_{rms}^6 \varepsilon ^{-2} a^{-1}$.

For characteristic velocity $U_\varLambda$, buoyant rise or turbulent advection can be relevant. The terminal buoyant rise velocity of a bubble in quiescent flow, $W_t$, is characterized by the bubble Reynolds number ${{Re}}_b \equiv 2 a W_t/\nu _w$. The bubbles in our DNS are in the inertial regime $1<{{Re}}_b<100$, where

(3.6)\begin{equation} W_{t} = 0.144\,{g}^{5/6}\,\nu_w^{2/3}\,{(2a)}^{3/2} {,} \end{equation}

with the prefactor coming from experiments (Wallis Reference Wallis1974; Park et al. Reference Park, Park, Lee and Lee2017). While this $W_{t}$ is appropriate for our bubbles, we note that there are different ${{Re}}_b$ regimes, and that corrections for surface effects are available (Park et al. Reference Park, Park, Lee and Lee2017). For small bubbles (small $W_t$), we expect the movement of the bubbles to be dominated by turbulent advection rather than buoyant rise. The dominance of turbulent advection implies a characteristic velocity for small bubbles $U_\varLambda \sim u_{rms}$ rather than $U_\varLambda \sim W_t$ for large bubbles. Assuming that $u_{rms}=W_{t}$ defines the transition between these two regimes, we solve (3.6) for $a$ to define the critical radius

(3.7)\begin{equation} a_{\varLambda} = 1.82\,u_{rms}^{2/3}\,\nu_w^{4/9}\,{g}^{-5/9} {.} \end{equation}

Although more advanced models consider the interaction between buoyancy and turbulence and incorporate surface tension effects (e.g. Salibindla et al. Reference Salibindla, Masuk, Tan and Ni2020; Ruth et al. Reference Ruth, Vernet, Perrard and Deike2021), we find that this simple two-regime model is sufficient to explain the power-law scaling of degassing at large ${{We}}$.

Using these $L_\varLambda$ and $U_\varLambda$ with $\varLambda (a)\propto U_\varLambda /L_\varLambda$, we obtain a split power law for $\alpha$. For large buoyancy-driven bubbles, $\alpha =5/2$, and for small turbulence-driven bubbles, $\alpha =1$:

(3.8)\begin{equation} \varLambda(a) \propto \varLambda_0 \begin{cases} (a/a_\varLambda), & a < a_\varLambda, \\ (a/a_\varLambda)^{5/2}, & a > a_\varLambda, \end{cases} \end{equation}

where $\varLambda _0 = 1.82\,u_{rms}^{-13/3}\,\varepsilon ^{2}\,\nu _w^{4/9}\,{g}^{-5/9}$. Thus, using (3.5) for degassing-dominated populations, we obtain

(3.9)\begin{equation} \beta = \begin{cases} -13/3, & a < a_\varLambda, \\ -35/6, & a > a_\varLambda. \end{cases} \end{equation}

4.1. Methodology

We study FST generated by an initial shear profile,

(4.1)\begin{equation} u(z, t=0)/U = 1 - 0.9988 \operatorname{sech}(0.88137 \,z/L) {,} \end{equation}

characterized by the shear velocity $U$ and shear length $L$. At sufficient Reynolds number ${{Re}}= U L/\nu _w$, a small, random initial perturbation causes the flow to become turbulent (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999). Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019, Reference Yu, Hendrickson and Yue2020) show that for Froude number squared (${{Fr}}^2= U^2/L {g}$) greater than a critical ${{Fr}}^2_{cr}\approx 5$, this FST entrains bubbles. There have also been extensive DNS characterizing this FST at low ${{Fr}}^2$ (Shen et al. Reference Shen, Zhang, Yue and Triantafyllou1999; Shen, Triantafyllou & Yue Reference Shen, Triantafyllou and Yue2000). Following all previous work, we choose ${{Re}}=1000$, and set the domain size to $(10.472 L)^2 \times 6 L$, where initially $z/L\in [-4,0]$ is water and $z/L\in [0,2]$ is air (density ratio $\rho _a/\rho _w=0.00123$, viscosity ratio $\mu _a/\mu _w=0.0159$). This gives an initial free-surface area $A/L^2=10.472^2$. For conciseness, we set $U$ and $L$ to unity. Thus all otherwise dimensional values that follow should be interpreted as normalized by $U$ and $L$ such that they are non-dimensional. The rendering in figure 1 comes from a simulation of ${{Fr}}^2=15$ at time $t=60$, with the shear velocity going from left to right.

For DNS, we use the same scheme as Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019, Reference Yu, Hendrickson and Yue2020). It is a pressure-projecting second-order staggered-grid finite-volume scheme, where air and water are evolved using the conservative volume-of-fluid method (Weymouth & Yue Reference Weymouth and Yue2010). Our interest is ${{We}}=U^2 L/(\sigma /\rho _w)\gg 1$ where surface tension is negligible, and Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) validate that not modelling surface tension accurately captures this limit. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) also study grid convergence to show that a grid $384^2 \times 256$ is sufficient for DNS of this flow, giving a grid size $\varDelta \approx 0.027$. We include an additional grid convergence study in the Appendix to confirm that this grid size sufficiently resolves entrainment and degassing. Following Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019), only bubbles of radius larger than $a_{res}=1.5 \varDelta$ are considered resolved.

For each ${{Fr}}^2$ studied (see table 1), we repeat the simulation with different realizations of the random initial perturbation to obtain ensemble statistics. We have at least six simulations at each ${{Fr}}^2$, and $\operatorname {O}(10^3)$ resolved entrainment and degassing events during the measured air-generating period for all but the smallest ${{Fr}}^2$ (close to ${{Fr}}^2_{cr}$). Because we normalize $U$ and $L$ to unity, ${g}$ is varied to obtain different ${{Fr}}^2$. A more intuitive interpretation of table 1 would be to consider a shear flow with fixed ${g}$ and length scale $L$, where larger ${{Fr}}^2$ mean larger $U^2$. Recall that reported units are $L$ for length and $L/U$ for time.

Table 1. Summary of FST simulations performed and values obtained during $t\in [40,70]$, including $95\,\%$ confidence intervals. Here, ${N}_{sim}$ is the number of ensemble simulations, ${N}_{I}$ is the number of entrainment events, and ${N}_{D}$ is the number of degassing events. Note that all otherwise dimensional quantities are normalized by the shear velocity $U$ and length $L$.

For comparing this shear flow FST to general FST, it is useful to use turbulence characteristics. We use a water domain spatial average

(4.2)\begin{equation} {\langle \cdot \rangle_\delta}{{\equiv} \frac{\iiint \cdot\ f \,\mathrm{d}\kern0.06em \boldsymbol{x}}{\iiint f \,\mathrm{d}\kern0.06em \boldsymbol{x}} }\quad\text{for } z>-\delta {,} \end{equation}

where $f$ is the volume of fluid ($f=1$ in water), and $\delta =0.2$ captures the near-surface HIT region (Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019). Figure 2 shows the turbulent dissipation rate $\varepsilon$ and turbulent kinetic energy $k$. As expected, there is little change in the strength of near-surface turbulence with ${{Fr}}$, and a quasi-steady period exists when $k$ remains approximately constant (Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019). We use $u_{rms}\equiv \sqrt {2 {\langle k \rangle _\delta }/3}$ and ${\langle \varepsilon \rangle _\delta }$ averaged over the quasi-steady period to calculate quantities from § 3.

Figure 2. Ensemble average (a) turbulent dissipation rate and (b) turbulent kinetic energy in the near-surface region, for solid lines: blue, ${{Fr}}^2=5$; light green, ${{Fr}}^2=8$; dark green, ${{Fr}}^2=10$; brown, ${{Fr}}^2=15$; red, ${{Fr}}^2=20$. Dashed lines indicate $t\in [40,70]$ over which we average to obtain the values in table 1.

While the values that we present are normalized by the shear flow $U$ and $L$, we note that the measured characteristic scales of the turbulence, $U_T\sim u_{rms}$ and $L_T\sim u_{rms}^3/\varepsilon$, allow comparison to general FST. For example, we calculate a turbulent Froude number squared, ${{Fr}}^2_{T} = \varepsilon / u_{rms} \,{g}$. Because ${\langle k \rangle _\delta }$ and ${\langle \varepsilon \rangle _\delta }$ do not change significantly with ${{Fr}}$ (see table 1), normalizing by the turbulence characteristic scales rather than the shear characteristic scales would represent only a change in a constant rather than a change in any scaling with ${g}$.

4.2. Measuring degassing and entrainment with ELA

The evolution of bubbles can be described using a volume tracking matrix (Gaylo et al. Reference Gaylo, Hendrickson and Yue2022), illustrated in figure 3. Consider bubbles at time $t^n$ labelled $j=1,\ldots, M^n$, where $M^{n}$ is the number of bubbles (including one for the air above the free surface). At time $t^{n+1}=t^n+T_s$, where $T_s$ is the snapshot interval, we have another set of bubbles labelled $i=1,\ldots, M^{n+1}$. The evolution of bubbles over $T_s$ is described by a volume tracking matrix, where each matrix element $q_{ij}$ gives the volume of air from a bubble $j$ at time $t^{n}$ that contributes to a bubble $i$ at time $t^{n+1}$. The volume tracking matrix provides a robust description of bubble evolution, and entrainment and degassing statistics can be extracted easily. We use ELA to track the movement of air volume on the grid level to (explicitly) calculate the volume tracking matrix. This method is proven to be volume-conservative, meaning that all air is tracked accurately (Gaylo et al. Reference Gaylo, Hendrickson and Yue2022; for a recent application, see Gaylo et al. Reference Gaylo, Hendrickson and Yue2023).

Figure 3. (a) Illustration of entrainment (red) and degassing (green). (b) The associated volume tracking matrix. At $t^{n}$, we have $j=1$ and three bubbles $j=2,3,4$. At $t^{n+1}$, we have $i=1$ and two bubbles $i=2,3$.

For ELA, we must choose the snapshot interval $T_s$. We want $T_s$ small because we miss the entrainment and degassing events of a bubble if they are separated by less than $T_s$. However, very small $T_s$ will generate spurious events because it is numerically unclear if bubbles (or a bubble and the free surface) are separate if their interfaces are within $\lesssim \varDelta$ (Chan et al. Reference Chan, Dodd, Johnson and Moin2021). To decrease spurious events, we use interface reconstruction to avoid connecting closely separated interfaces (Hendrickson, Weymouth & Yue Reference Hendrickson, Weymouth and Yue2020). In general, $T_s\, \varOmega (a_{res}) > 0.1$ avoids spurious events (Chan et al. Reference Chan, Dodd, Johnson and Moin2021; Gaylo et al. Reference Gaylo, Hendrickson and Yue2022). For our $\varepsilon$, this gives $T_s>0.4$. We perform ELA with $T_s=0.4$, $0.8$ and $1.6$. A posteriori, we see slight decreases in $Q_I$ and $Q_D$ as $T_s$ increases, as expected, due to missing temporally close entrainment and degassing events. We find that the shapes of $I(a)$ and $D(a)$ are independent of $T_s$, suggesting that spurious events are negligible. Thus we select the smallest $T_s=0.4$.