2.1. Governing equations

The flow field comprising water and air bubbles is treated as a one-fluid incompressible mixture. The volume-of-fluid (VOF) method is employed to capture the bubble interfaces. The governing equations in Cartesian coordinates for the liquid-phase volume fraction $\alpha _l$ , the mass and momentum conservation read

(2.1) \begin{equation} \frac {{\partial \alpha _l }}{{\partial t}} + {u_i}\frac {{\partial \alpha _l }}{{\partial {x_i}}} = 0, \end{equation}

(2.2) \begin{equation} \frac {{\partial {u_i}}}{{\partial {x_i}}} = 0, \end{equation}

(2.3) \begin{equation} \frac {{\partial ({\rho _m}{u_i}})}{{\partial t}} + \frac {{\partial \left ( {{\rho _m}{u_i}{u_{\!j}}} \right )}}{{\partial {x_{\!j}}}} = - \frac {{\partial p}}{{\partial {x_i}}} + \frac {\partial \tau _{ij}^{\textit{eff}}}{\partial x_{\!j}} + {F_i^\sigma } + {F_i^{\textit {rep}}}, \end{equation}

where $p$ is the pressure field, $u_i$ is the velocity vector and $(u,v,w)$ denote the components of the velocity vector in the $x$ , $y$ and $z$ directions, respectively. Gravity is neglected in the present work to focus on the interaction between bubbles and the turbulent boundary layer. The effective stress tensor, including both the molecular and subgrid-scale (SGS) contributions, is modelled as

(2.4) \begin{equation} \tau _{ij}^{\textit{eff}} = \tau _{ij}+\tau _{ij}^{\textit {SGS}} = \left ( {\mu _m} + {\mu _t}\right ) \left ( \frac {\partial u_i}{\partial x_{\!j}}+ \frac {\partial u_{\!j}}{\partial x_i} \right )\!. \end{equation}

The mixture density and viscosity are calculated using a volume-fraction-weighted approach: $\rho _m=\alpha _l\rho _l+\alpha _g\rho _g$ and $\mu _m=\alpha _l\mu _l+\alpha _g\mu _g$ , with the subscripts $l$ and $g$ denoting the liquid and gas phases, respectively. The two-phase volume fraction satisfies $\alpha _l+\alpha _g=1$ . Considering the necessity to resolve both the bubbles and boundary layer in order to capture precise pressure fluctuations, we employ the wall-adapting local eddy-viscosity (WALE) LES (Nicoud & Ducros Reference Nicoud and Ducros1999) to relax the requirement for extremely fine meshes near the wall. Thus, the turbulent viscosity $\mu _t=\rho _m\nu _t$ is included in the momentum equation. A thorough validation is performed by comparing the single-phase results with DNS data for velocity statistics (Lee & Moser Reference Lee and Moser2015) and pressure fluctuations spectra (Yang & Yang Reference Yang and Yang2022), demonstrating the high fidelity of the present LES (see Appendix A). The influence of the WALE model on bubbly flows is assessed in Appendix B in terms of the distributions of turbulent viscosity and SGS stresses. The time- and space-averaged turbulent kinematic viscosity $\langle \nu _t \rangle$ remains below $0.08$ times the water kinematic viscosity, while the SGS stresses are at least three orders of magnitude smaller than the resolved Reynolds stresses throughout the channel. Furthermore, the second-order velocity moments and the wall-pressure spectra are in close agreement with the corresponding under-resolved DNS (uDNS) results obtained without the explicit SGS model. These results demonstrate that the contribution of the WALE model remains negligible under the present grid resolution and does not influence the principal conclusions of the study. Nevertheless, the LES framework is retained for future investigations at higher Reynolds numbers.

The surface tension is modelled using the continuum surface force model (Brackbill, Kothe & Zemach Reference Brackbill, Kothe and Zemach1992) as $F_i^\sigma =\sigma \kappa \partial \alpha _l/\partial x_i$ , where $\sigma$ is the constant surface-tension coefficient, $\kappa =- \partial n_i/ \partial x_i$ is the curvature of the interface and $n_i=\boldsymbol{\nabla }\alpha / \left | \boldsymbol{\nabla }\alpha \right |$ is the unit normal vector of the interfaces. The VOF method is prone to spurious numerical coalescence between bubbles (Zhang, Ni & Magnaudet Reference Zhang, Ni and Magnaudet2021; Innocenti et al. Reference Innocenti, Jaccod, Popinet and Chibbaro2021), requiring an extremely high grid resolution to prevent this unphysical phenomenon. An alternative way is to apply a short-range repulsive force on the interface through calibrating the artificial coefficient (Zhang et al. Reference Zhang, Peng, Shao and Deng2022). Considering the trade-offs, we adopt the latter approach with the coefficient $K=1\times 10^{-8}\,\textrm {J}$ to prevent unphysical bubble coalescence while avoiding excessive influence on the overall bubble motion. The repulsive force $F_i^{\textit{rep}}$ is incorporated into the source term of the momentum equation (2.3).

The simulations are performed using the open-source code OpenFOAM, employing the two-phase solver interFoam. The governing equations are discretised using the finite-volume method. Time advancement is performed with a first-order implicit Euler scheme. The divergence terms are evaluated using the Gauss theorem, which converts the cell-volume integrals into summations of face fluxes. The fluxes and the field quantities required at cell faces are obtained using second-order linear interpolation. For the phase-fraction equation, the advective flux is discretised using the second-order bounded van Leer scheme. An artificial compression term is introduced in the phase-fraction equation to counteract numerical diffusion of the interface; the corresponding compression flux is evaluated using second-order linear interpolation. The momentum convection term is discretised using a bounded flux-limited linear scheme, which retains second-order accuracy in smooth regions and locally reverts towards upwind differencing near steep gradients to improve boundedness and stability. The viscous stress divergence is also evaluated from face-flux summations, with the relevant quantities at cell faces obtained using second-order linear interpolation.

2.2. Computational set-up

The channel consists of two no-slip walls perpendicular to the direction $y$ and periodic boundaries in the streamwise $x$ and spanwise $z$ directions. The computational domain is defined as $L_x \times L_y \times L_z = 4h \times 2h \times 2h$ , where $h$ is the channel half-height. This domain size is slightly larger than the ‘minimum turbulent channel‘ ( $\pi h \times 2h \times \pi h/2$ ) in Lu et al. (Reference Lu, Fernández and Tryggvason2005). Extending the channel configuration can facilitate the discovery of large-scale flow structures (Santarelli & Fröehlich Reference Santarelli and Fröehlich2015; Cifani, Kuerten & Geurts Reference Cifani, Kuerten and Geurts2018; Nemati et al. Reference Nemati, Breugem, Kwakkel and Boersma2021). However, the large computational domain appears formidable in the present simulation. In Appendix C we calculate the two-point correlation function in space, near the wall and at the channel centre, respectively, to validate the computational domain size. The computational domain is discretised with a structured mesh. Bubbles with diameters of $D_b=2$ and 3 mm are considered in this study. From the grid-convergence analysis (see Appendix D), it is found that 30 grid cells per bubble diameter are sufficient to achieve convergence in the distribution of the gas volume fraction, velocity profiles and second-order velocity moments. The channel half-height is set to $h=0.01$ m, resulting in the ratio $D_b/h = 0.2$ and 0.3. For $D_b=2$ mm, the computational domain is meshed with $N_x \times N_y \times N_z = 600 \times 400 \times 300$ , corresponding to a wall-unit resolution of $\varDelta x^+ = \varDelta z^+ = 3.67$ in the streamwise and spanwise directions. The grid employs a uniform expansion rate in the wall-normal $y$ direction to enhance near-wall resolution, resulting in $\varDelta y^+ |_{y=0}= \varDelta y^+ |_{y=2h}=1.85$ and $\varDelta y^+|_{y=h}/ \varDelta y^+ |_{y=0}= 2.1$ . For $D_b=3$ mm, the numbers of the mesh are $N_x \times N_y \times N_z = 400 \times 300 \times 200$ with resolution $\varDelta x^+ = \varDelta z^+ = 5.5$ , $\varDelta y^+ |_{y=0}= \varDelta y^+ |_{y=2h}=1.47$ and $\varDelta y^+|_{y=h}/ \varDelta y^+ |_{y=0}= 5$ . These details of the mesh configurations and parameters of cases are listed in table 2. The number 64 or 128 in the case name denotes the bubble count and ‘D2’ or ‘D3’ denotes the bubble diameter. In this paper, the superscript ’+’ denotes variables normalised using wall units, i.e. viscous length $\delta _\nu =\nu _l /u_\tau$ and wall-friction velocity $u_\tau$ as the characteristic length and velocity scales, respectively, where $\nu _l=\mu _l/\rho _l$ is the kinematic viscosity of the liquid and $\mu _l$ is the dynamic viscosity of the liquid. The density and viscosity for water and air under standard atmospheric pressure and at 20 $^\circ$ C are employed, i.e. $\rho _l=998.18\,{\textrm {kg}\,{\textrm {m}}^{-3}}$ , $\nu _l = 1.00\times 10^{-6}\,{\textrm {m}^2\textrm {s}^{-1}}$ , $\rho _g=1.1881\, {\textrm {kg}\,{\textrm {m}}^{-3}}$ and $\nu _g = 1.83\times 10^{-5}\,{\textrm {m}^2\textrm {s}^{-1}}$ . Naturally, the interfacial tension between the gas and liquid phases adopts the standard value at ambient temperature, i.e. $\sigma = 0.072\,{\textrm {N m}^{-1}}$ .

Table 2.

Simulation parameters for all cases. The channel half-height is $h = 0.01$ m. The mesh for case B128D2 is identical to that for B64D2, and the mesh for B128D3 is identical to that for B64D3.

During the simulations, the Courant–Friedrichs–Lewy number is maintained at 0.5. The interfacial Courant number, which is defined based on the volumetric flux magnitude in the vicinity of the interface and provides a measure of the local interface advection speed relative to the cell size, is also restricted to 0.5 for all bubbly cases. In addition, surface-tension resolution is assessed using the cell-surface Weber number $We_\varDelta = \rho _l u_{\textit{rms}}^2 \varDelta / 4\pi \sigma$ (Popinet Reference Popinet2018). In the present work, $We_\varDelta \lt 0.001$ , ensuring surface-tension forces are well resolved.

Since the instantaneous turbulent flow fields, especially in bubbly cases, are not perfectly symmetric up and down, the wall shear stress $\tau _w$ is computed based on the streamwise velocity profile and averaged with the values at the bottom and top walls, i.e.

(2.5) \begin{equation} {\tau _w} = \frac {{{{\left . {{\tau _w}} \right |}_{y = 0}} + {{\left . {{\tau _w}} \right |}_{y = 2h}}}}{2} = \frac {1}{2}\left ( {{{\left . {{\mu _l}\frac {{{\mathrm{d}}\left \langle u \right \rangle }}{{{\mathrm{d}}y}}} \right |}_{y = 0}} - {{\left . {{\mu _l}\frac {{{\mathrm{d}}\left \langle u \right \rangle }}{{{\mathrm{d}}y}}} \right |}_{y = 2h}}} \right )\!, \end{equation}

and the wall-friction velocity $u_\tau$ is calculated based on the above wall shear stress as

(2.6) \begin{equation} {u_\tau } = \sqrt {\frac {{{\tau _w}}}{{{\rho _l}}}}. \end{equation}

The friction Reynolds number is calculated by $R{e_\tau } = {{u_\tau }h}/{{\nu _l}}$ . The pair of angular brackets in (2.5) denotes averaging over time and the streamwise–spanwise plane, i.e.

(2.7) \begin{equation} \left \langle \phi \right \rangle (y) = \frac {1}{{\varDelta {T}{L_x}{L_z}}}\int _t^{t + \varDelta {T}} {\int _0^{{L_x}} {\int _0^{{L_z}} {\phi \left ( {x,y,z,t} \right ){\mathrm{d}}z{\mathrm{d}}x{\mathrm{d}}t} } }, \end{equation}

where $\varDelta T$ is the sampling period. The averaging in (2.7) results in a function that varies only with the $y$ coordinate.

To facilitate the subsequent discussion, we introduce the decomposition of the field into the mean and fluctuation components from the perspective of Reynolds decomposition. A field $\phi (x,y,z)$ can thus be expressed as

(2.8) \begin{equation} \phi (x,y,z,t) = \overline {\phi }(x,y,z) + \phi '(x,y,z,t), \end{equation}

where the overbar denotes time averaging applied to the mixture field (i.e. without phase discrimination) and the prime denotes the corresponding fluctuation about this mixture-averaged mean. It should be noted that this approach differs from the phase-averaged decomposition employed in some studies listed in table 1 (e.g. Bolotnov et al. (Reference Bolotnov, Jansen, Drew, Oberai, Lahey and Podowski2011), Lakehal et al. (Reference Lakehal, Métrailler and Reboux2017), Cifani et al. (Reference Cifani, Kuerten and Geurts2018), du Cluzeau et al. (Reference du Cluzeau, Bois and Toutant2019), Cifani et al. (Reference Cifani, Kuerten and Geurts2020), Trautner et al. (Reference Trautner, Klein, Bräuer and Hasslberger2021)); consequently, the resulting fluctuations $\phi ^\prime$ contain contributions from both turbulent motions and interfacial dynamics.

Channel flow must be driven by a streamwise pressure gradient to balance the shear stress. Some studies employed a fixed pressure gradient (Lu & Tryggvason Reference Lu and Tryggvason2006; Liu et al. Reference Liu, Wang, Yang, Nemati and Chu2023), while others used temporally variable pressure gradients to maintain fixed parameters, such as the friction Reynolds number (du Cluzeau et al. Reference du Cluzeau, Bois and Toutant2019, Reference du Cluzeau, Bois, Toutant and Martinez2020), bulk velocity (Cifani et al. Reference Cifani, Kuerten and Geurts2020) or mass flux (Hasslberger et al. Reference Hasslberger, Cifani, Chakraborty and Klein2020; Klein et al. Reference Klein, Trummler and Radtke2022). In this paper, a fixed pressure gradient is applied for the single-phase channel case, while temporally variable pressure gradients are imposed for bubbly cases to ensure a consistent bulk velocity $U_b$ . Thus, the bulk Reynolds numbers $\textit{Re}_{_B} = 2h U_b /\nu _l$ are uniform for all cases.

A single-phase turbulent channel flow filled with water is first simulated under the friction Reynolds number $\textit{Re}_\tau =550$ . The velocity field is initialised using a Blasius profile with prescribed perturbations following De Villiers (Reference De Villiers2006) to accelerate the transition to a fully developed turbulent flow. Integrating $ \left \langle u \right \rangle$ in the $y$ direction yields the bulk mean velocity $U_b$ and the corresponding bulk Reynolds number is $\textit{Re}_{_B} = 20\,700$ , which will be applied to the bubbly flows. Consequently, the friction Reynolds number will vary within a narrow range among different bubbly cases. After a transient period, the turbulence reaches a statistically stationary state, and the data are sampled over a duration of $30L_x/U_b$ for further time averaging and post-processing.

The instantaneous velocity field at the end of the single-phase channel flow simulation is used as the initial velocity field for the bubbly cases. The influence of the initial bubble positions has been examined. After release, the bubbles undergo significant deformation and subsequently migrate toward the channel centreline, largely independent of their initial release locations. To shorten the transient period and reduce the simulation time, spherical bubbles are initially introduced into the flow field near the channel core in all subsequent bubbly cases. After initialisation, each simulation is advanced for an additional $7.77\,L_x/U_b$ to ensure that the system reaches a statistically stationary state. Subsequently, data are collected over an averaging duration of $38.86\,L_x/U_b$ . The statistical convergence with this averaging duration is examined in Appendix E, where the temporal stabilisation of first- and second-order quantities as well as the symmetry between the upper and lower halves of the channel are demonstrated.

4.1. Mean pressure and pressure fluctuation intensity

The pressure fluctuation is defined as the instantaneous pressure minus the mean pressure. We therefore first examine the distribution of the mean pressure, $\langle p\rangle$ . Taking the spatio-temporal average of (2.3), neglecting the repulsive force and considering the wall-normal component, the averaged momentum equation can be simplified as

(4.1) \begin{equation} 0 = - \frac {{\mathrm{d}\langle p\rangle }}{{\mathrm{d}y}} - \frac {{\mathrm{d}\langle \rho _m v'v'\rangle }}{{\mathrm{d}y}} + \frac {{\mathrm{d}\langle \tau ^{\textit{eff}}_{yy}\rangle }}{{\mathrm{d}y}} + \langle F^{\sigma }_{y}\rangle . \end{equation}

Here, the mean viscous contribution in $\langle \tau _{yy}\rangle$ is of the order of $10^{-3}\rho _l u_\tau ^2$ , while the SGS stress $\langle \tau ^{\textit{SGS}}_{yy}\rangle$ is of the order of $10^{-4}\rho _l u_\tau ^2$ , as shown in Appendix B. The contribution associated with the effective wall-normal normal stress is several orders of magnitude smaller than $\langle \rho _m v'v'\rangle /(\rho _l u_\tau ^2)$ , thus, is negligible in both the single-phase and bubbly flows.

For the single-phase flow, the surface-tension term is absent. The wall-normal gradient of mean pressure is then balanced by the gradient of the turbulent momentum. Integrating (4.1) gives

(4.2) \begin{equation} \langle p\rangle = \textrm {const.} - \rho _l \langle v'v'\rangle . \end{equation}

Thus, the mean pressure profile follows the distribution of $-\langle v'v'\rangle$ , as shown in figure 6(a), where $\langle p\rangle$ has been shifted by a constant such that $\langle p\rangle =0$ at the wall.

For the bubbly flows, $\langle \rho _m v'v'\rangle$ remains close to its single-phase counterpart, as shown in figure 5. Therefore, near the wall, where the gas volume fraction is small and the mean surface-tension force $\langle F^{\sigma }_{y}\rangle$ is weak, the shape of $\langle p\rangle$ is similar to that in the single-phase flow and exhibits a local minimum. In regions with a larger gas volume fraction, however, $\langle F^{\sigma }_{y}\rangle$ becomes important, leading to a mean pressure profile that resembles the distribution of $\langle \alpha _g\rangle$ . Physically, this behaviour reflects the pressure jump across the gas–liquid interface induced by surface tension: in regions where bubbles are more frequently present, the averaged pressure is correspondingly increased.

The profiles of the root-mean-square (r.m.s.) pressure fluctuations $p^\prime _{\textit{rms}} = {\langle \sqrt { \overline {p^{\prime 2}}}\rangle }$ , shown in figure 6(b) in solid lines, demonstrate that the pressure fluctuations in the bubbly flows are substantially larger than that in the single-phase flow except in the very near-wall region. At the wall, $p^\prime _{\textit{rms}}(0)/\tau _w$ equals 2.25 (SP), 2.09 (B64D2), 2.08 (B128D2), 2.06 (B64D3) and 2.19 (B128D3), which are similar across all cases and are in reasonable agreement with the semi-empirical prediction of wall-pressure fluctuations 2.73 for single-phase flow (Farabee & Casarella Reference Farabee and Casarella1991), which is based on experimental and numerical data. A detailed analysis of the wall-pressure fluctuations is presented in § 5 based on wavenumber–frequency spectra and the associated convection velocity. In particular, the discussion of the 1-D spectra in § 5.4 will explain why $p^\prime _{\textit{rms}}$ at the wall exhibits much smaller differences between single-phase and bubbly flows than those observed at other locations.

4.2. Decomposition of the instantaneous pressure fluctuation field

Referring to the Poisson equation of pressure fluctuations for variable-density flows derived by Chassaing et al. (Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002), the pressure fluctuation in the present multiphase incompressible flows is decomposed into the components

(4.3) \begin{equation} p' = p'_r + p'_s + p'_c + p'_v + p'_\sigma , \end{equation}

including the rapid pressure fluctuation $p'_r$ , slow pressure fluctuation $p'_s$ , variable-density-related pressure fluctuation $p'_c$ , viscosity-related pressure fluctuation $p'_v$ and the surface-tension-related pressure fluctuation $p'_\sigma$ . Each component satisfies the Poisson equation associated with the corresponding source term:

(4.4) \begin{align} \boldsymbol{\nabla} ^2 p'_r = S_r = -2\frac {\partial ^2 \left (\overline {\rho _m}\,\bar {u}_{\!j}u'_i \right )} {\partial x_i \partial x_{\!j}}, \end{align}

(4.5) \begin{align} \boldsymbol{\nabla} ^2 p'_s = S_s = -\frac {\partial ^2 \left (\overline {\rho _m}u'_i u'_{\!j} - \overline {\rho _m}\,\overline {u'_i u'_{\!j}}\right )} {\partial x_i \partial x_{\!j}}, \end{align}

(4.6) \begin{equation} \boldsymbol{\nabla} ^2 p'_c = S_c = \frac {\partial ^2 \rho '_m}{\partial t^2} -2\frac {\partial ^2 \left [\left (\rho '_m u'_i - \overline {\rho '_m u'_i}\right )\bar {u}_{\!j}\right ]} {\partial x_i \partial x_{\!j}} -\frac {\partial ^2 \left (\rho '_m u'_i u'_{\!j} - \overline {\rho '_m u'_i u'_{\!j}}\right )} {\partial x_i \partial x_{\!j}} -\frac {\partial ^2 \left (\rho '_m \bar {u}_i \bar {u}_{\!j}\right )} {\partial x_i \partial x_{\!j}}, \end{equation}

(4.7) \begin{align} \boldsymbol{\nabla} ^2 p'_v = S_v = \frac {\partial ^2 \left (\tau ^{\textit{eff}}_{ij}-\overline {\tau ^{\textit{eff}}_{ij}}\right )} {\partial x_i \partial x_{\!j}}, \end{align}

(4.8) \begin{align} \boldsymbol{\nabla} ^2 p'_\sigma = S_\sigma = \frac {\partial }{\partial x_i} \left (\sigma \kappa \frac {\partial \alpha }{\partial x_i}\right ) - \overline { \frac {\partial }{\partial x_i} \left (\sigma \kappa \frac {\partial \alpha }{\partial x_i}\right ) }. \end{align}

The right-hand-side source term in the above Poisson equations are extracted from the numerical data. At the walls, the viscous component satisfies the boundary condition $\partial p'_v / \partial y = \partial \tau _{yi} / \partial x_i$ , whereas homogeneous Neumann boundary conditions are imposed for the remaining Poisson equations.

The first two equations, (4.4) and (4.5), take the same form as the classical Poisson equation of pressure fluctuations for single-phase incompressible flows (Kim Reference Kim1989), solving for the rapid pressure $p'_r$ with the linear source term $S_r$ and the slow pressure $p'_s$ with the nonlinear source term $S_s$ , respectively. The difference from the single-phase case lies in the use of mean density $\overline {\rho _m}$ in the present formulation. The solutions to these two Poisson equations are first performed for the single-phase case. The instantaneous contours of the pressure fluctuations associated with these two components are shown in figures 7(c) and 7(d). The slow component exhibits a stronger intensity than the rapid one, consistent with the well-established dominance of the slow pressure contribution in turbulent channel flows. To further verify the correctness of the solution, the superposition of these two components produces the pressure fluctuation shown in figure 7(b). The reconstructed field agrees almost perfectly with the pressure fluctuation obtained directly from the simulation (figure 7 a). This agreement confirms the reliability of the evaluation of the two incompressible source terms.

Figure 7.

Instantaneous pressure fluctuation fields for the single-phase case: (a) instantaneous pressure fluctuations $p^\prime$ ; (b) reconstructed field $p^\prime _r + p^\prime _s$ ; (c) the rapid component $p^\prime _r$ ; (d) the slow component $p^\prime _s$ . The wall-parallel plane is located at $y/h = 0.1$ and the z-normal plane is located at $z/h = 0.1$ .

The viscous source term $S_v$ in (4.7) is computed using the same viscous stress tensor appearing in the momentum equation, thereby ensuring numerical consistency.

The surface-tension-related source $S_\sigma$ in (4.8) is derived from the formulation of the momentum equation (2.3). In the present work, the term $\boldsymbol{\nabla }\boldsymbol{\cdot }(\sigma \kappa \boldsymbol{\nabla }\alpha )$ is evaluated using a face-based finite-volume discretisation. The product of curvature $\kappa$ and surface-tension coefficient $\sigma$ is interpolated to the cell faces, and the face-normal gradient $\boldsymbol{\nabla} ^\perp _f \alpha$ is computed using the standard face-normal gradient operator. These quantities are assembled into a surface flux $(\sigma \kappa )_f (\boldsymbol{\nabla} ^\perp _f \alpha ) S_f$ , where the subscript $f$ denotes the face of the finite-volume cell and $S_f$ is the face area magnitude. The divergence is then obtained through a Gauss integration over the cell faces and normalised by the cell volume. This discretisation is identical to the finite-volume operators employed in the discrete solution of the governing equations (2.1)–(2.3), ensuring full consistency with the underlying numerical formulation.

The terms involving the density fluctuations are grouped into the variable-density source $S_c$ . The form is consistent with the compressible terms in Chassaing et al. (Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002), while here the effects arise from density variations associated with the motion of bubbles. Because the mixture density exhibits sharp variation across the gas–liquid interface, evaluating the individual density-fluctuation-related terms is extremely sensitive to discretisation used in post-processing due to the presence of a second derivative, and may lead to significant numerical errors. Therefore, the combined contribution from density variation $S_c$ is estimated through the residual of the Poisson equation as $S_c = \boldsymbol{\nabla} ^2 p' - S_r - S_s - S_v - S_\sigma$ .

Figure 8.

Instantaneous pressure fluctuation fields for case B128D3: (a) instantaneous pressure fluctuations $p^\prime$ ; (b) reconstructed field $p^\prime _r + p^\prime _s + p^\prime _c + p^\prime _v + p^\prime _\sigma$ ; (c) rapid-linear component $p^\prime _r$ ; (d) slow-nonlinear component $p^\prime _s$ ; (e) variable-density-related component $p^\prime _c$ ; (f) viscosity-related component $p^\prime _v$ ; (g) surface-tension-related component $p^\prime _\sigma$ . The wall-parallel plane is located at $y/h=0.5$ , and the $z$ -normal plane is located at $z/h=0.1$ .

Figure 8. Long description

Panel A: A heat map representing instantaneous pressure fluctuation. The color scale ranges from -20 to 20, with darker colors indicating lower values and lighter colors indicating higher values. The heat map shows a distribution of pressure fluctuations with some regions of high intensity. Panel B: A heat map showing the reconstructed field of pressure fluctuations. The color scale ranges from -20 to 20, similar to Panel A, with a similar distribution of high and low intensity regions. Panel C: A heat map depicting the rapid-linear component of pressure fluctuations. The color scale ranges from -10 to 10, with a more uniform distribution and fewer high-intensity regions compared to Panels A and B. Panel D: A heat map illustrating the slow-nonlinear component of pressure fluctuations. The color scale ranges from -10 to 10, showing a more dispersed pattern of fluctuations. Panel E: A heat map showing the variable-density-related component of pressure fluctuations. The color scale ranges from -20 to 20, with distinct regions of high and low intensity. Panel F: A heat map depicting the viscosity-related component of pressure fluctuations. The color scale ranges from -1 to 1, indicating a more subtle variation in pressure fluctuations. Panel G: A heat map illustrating the surface-tension-related component of pressure fluctuations. The color scale ranges from -20 to 20, with prominent regions of high intensity.

Figure 8 present the instantaneous pressure fluctuation decomposition for case B128D3 as an example. The reconstructed pressure fluctuation obtained from the Poisson equations, $p'_r + p'_s + p'_c + p'_v + p'_\sigma$ (figure 8 b), is compared with the pressure fluctuation directly obtained from the numerical simulation (figure 8 a). The two fields exhibit excellent agreement, both in their spatial structures and in quantitative pointwise comparisons, thereby demonstrating the accuracy of the Poisson solutions in (4.4)–(4.8) and the resulting pressure fluctuation components. In general, the pressure fluctuations are primarily associated with the presence of bubbles, with the largest magnitudes occurring in the region surrounding the bubbles.

The rapid component $p'_r$ originates from the interaction between velocity fluctuations and the mean shear. Strong velocity fluctuations occur across the bubble interfaces due to the discontinuous flow characteristics inside and outside the bubbles. As a result, the rapid component is significantly stronger around the bubbles than in the pure liquid phase. In comparison, the slow component $p'_s$ provides the dominant contribution to the pressure fluctuations within the liquid regions, although the largest magnitudes still occur near the bubbles.

The variable-density component $p'_c$ is closely associated with density fluctuations induced by bubble motion. A characteristic pattern is observed in figure 8(e) around the bubbles: the pressure fluctuation is negative on the upstream side of the bubble and positive on the downstream side, consistent with the direction of bubble motion. Even in regions without bubbles, $\rho '_m$ remains non-zero because the mean density is less than the liquid density.

The viscous contribution $p'_v$ , as shown in figure 8(f), remains small throughout the channel, with an amplitude more than one order of magnitude smaller than those of the other components. This behaviour is expected because it is associated with the Stokes pressure, which vanishes in incompressible flows and remains negligible in compressible flows (Gerolymos, Sénéchal & Vallet Reference Gerolymos, Sénéchal and Vallet2013).

The surface-tension-related component $p'_\sigma$ clearly reveals the pressure difference between the interior and exterior of the bubbles, shown in figure 8(g), which is physically consistent with the Laplace pressure jump across the gas–liquid interface. This surface-tension contribution accounts for the dominant part of the pressure fluctuations in the bubbly flow; a similar prominence has been found in the momentum-equation budgets (du Cluzeau et al. Reference du Cluzeau, Bois and Toutant2019). As a result, the r.m.s. pressure fluctuations $p'_{\textit{rms}}$ in the bubbly flows becomes significantly larger over most of the channel compared with the single-phase case (figure 6). Both a higher gas volume fraction and a larger surface-tension force increase the capillary pressure jump across the bubble interface, which elevates the level of $p'_{\textit{rms}}$ .

4.3. Comparison of the source magnitudes of pressure fluctuations

The instantaneous pressure fluctuation decomposition discussed above reveals the spatial organisation of the different pressure fluctuation components. To further connect these pressure fields with their corresponding Poisson sources and to assess the relative importance of the different source mechanisms, we next examine the source terms in the pressure fluctuation Poisson equation.

We first consider their instantaneous spatial distributions. Figure 9 shows the instantaneous source terms for case B128D3 at the same time instant as that used in figure 8. The colour ranges are clipped to improve the visibility of the spatial structures. The total source term, $\boldsymbol{\nabla} ^2 p'$ , is mainly concentrated in the vicinity of the gas–liquid interfaces. Additional localised structures can also be observed in the small bubble wake regions and near the wall. The strong interfacial distribution results from the combined contributions of the different source terms, whereas the localised wake structures mainly arise from $S_r$ and $S_c$ . The near-wall distribution of $\boldsymbol{\nabla} ^2 p'$ is mainly originated from the slow source term $S_s$ . The viscous source term $S_v$ has a much smaller magnitude, consistent with the weak contribution of $p'_v$ shown above. In particular, $S_\sigma$ is sharply localised at the gas–liquid interfaces and forms ring-like structures exactly along the bubble interfaces. These interfacial source structures, when inverted through the Poisson equation, give rise to the spatial distribution of $p'_\sigma$ , which represents the surface-tension-induced pressure jump across the bubble interface.

Figure 9.

Instantaneous source terms for case B128D3 at the same time instant as the snapshot shown in figure 8. (a) Laplacian of the pressure fluctuation, $\boldsymbol{\nabla} ^2 p^\prime$ ; (b) rapid source term, $S_r$ ; (c) slow source term, $S_s$ ; (d) variable-density source term, $S_c$ ; (e) viscous source term, $S_v$ ; and (f) surface-tension-related source term, $S_\sigma$ . The wall-parallel plane is located at $y/h=0.5$ , and the $z$ -normal plane is located at $z/h=0.1$ .

Figure 10.

Profiles of the mean square of the source terms in the Poisson equations (4.4)–(4.8) for case B128D2 (solid lines) and B128D3 (dashed lines).

Figure 10. Long description

Panel: A line graph showing the profiles of the mean square of the source terms in the Poisson equations for cases B 128 D 2 and B 128 D 3. The horizontal axis represents y plus, ranging from 10 to the power of 0 to 10 to the power of 2. The vertical axis represents the mean-square of source terms divided by tau w squared delta n v to the power of 4, ranging from 10 to the power of -8 to 10 to the power of 0. The graph includes multiple lines representing different source terms: B 128 D 2 ((S r plus S s) squared) in red solid line, B 128 D 2 (S c squared) in blue solid line, B 128 D 2 (S v squared) in orange solid line, B 128 D 3 ((S r plus S s) squared) in red dashed line, and B 128 D 2 (S sigma squared) in yellow dashed line. The graph is divided into three regions: Viscous sublayer, Buffer layer, and Log layer. Each line shows different trends across these regions.

Figure 11.

Profiles of the mean square of the incompressible source terms in the Poisson equations (4.4) and (4.5) for case B128D2, B128D3 and SP.

We then examine the intensity of the corresponding source terms in a statistical sense, taking cases B128D2 and B128D3 as examples. The spatio-temporally averaged mean-square values of the source terms as a function of the wall-normal coordinate $y^+$ are shown in figure 10. The following discussion is organised according to the classical division of the near-wall region into the viscous sublayer, buffer layer and logarithmic layer.

Because only three grid points fall within the viscous sublayer, the detailed variation of the source terms in this region cannot be fully resolved, and the corresponding trends should therefore be interpreted with caution. Nevertheless, the available results suggest that the viscous contribution is comparable to the other terms very close to the wall. Moving away from the viscous sublayer, however, its magnitude rapidly decreases and becomes several orders of magnitude smaller than the other contributions. This behaviour is consistent with the fact that viscous effects are confined to the smallest scales and are primarily important in the immediate vicinity of the wall. The surface-tension source $S_\sigma$ becomes the dominant contribution among all source terms. The surface-tension term is determined by the bubble distribution and the interface curvature, which vary with bubble diameter and deformation across the channel. Bubbles in the near-wall region exhibit evident deformation, as shown in figure 2(c), resulting in large values of $S_\sigma$ . Although the mean gas volume fraction $\langle \alpha _g \rangle$ is low in this region, instantaneous flow-field visualisations in figure 1 reveal sporadic bubble passages, which may generate high local values in $S_\sigma$ . Consequently, the layer averaged $\langle S^2_\sigma \rangle$ in this near-wall region exceeds that of other source terms because of these intermittent yet intense contributions.

In the buffer layer, all source terms increase in magnitude compared with the near-wall region. A clear difference between cases B128D3 and B128D2 can already be observed, with larger values in B128D3. This trend is broadly consistent with the higher gas volume fraction and indicates a stronger bubble-induced modulation of the pressure fluctuation field.

In the logarithmic layer, the incompressible contribution (including both the rapid and slow components), the variable-density terms and the surface-tension-related terms are of comparable magnitude and together dominate the pressure fluctuations. This balance highlights that, in bubble-rich regions, the pressure fluctuation is governed by a combination of mechanisms, including shear-induced turbulent interactions, density variations associated with bubble motion and interfacial surface-tension effects. In contrast, the viscous term remains negligible throughout this region.

Overall, these results demonstrate that, unlike in single-phase flows where pressure fluctuations are governed by the incompressible rapid and slow terms, bubbly flows exhibit a multi-mechanism balance in which variable-density and surface-tension effects play equally important roles in the generation of pressure fluctuations.

Figure 11 shows the mean-square profiles of the incompressible source terms for cases B128D2, B128D3 and the single-phase flow (SP). In the bubbly flows, the incompressible sources are smaller than those in the single-phase case within the viscous sublayer, but they exhibit a markedly different trend farther away from the wall. In the single-phase case, the incompressible source terms attain their maximum within the buffer layer and subsequently decrease towards the channel centre. In contrast, for the bubbly cases, the source terms continue to increase beyond the buffer layer and reach their largest values in the logarithmic and core regions. In the channel core, the incompressible source terms become significantly larger than those in the single-phase flow. This behaviour is consistent with the instantaneous observations discussed in figures 7 and 8, where $p'_r$ and $p'_s$ remain strong in the core region of the bubbly flow, while in single-phase flow they both decay in the channel core.

A further distinction arises in the relative importance of the linear and nonlinear contributions. In bubbly cases, the linear term becomes larger in the logarithmic region ( $y^+\gt 55$ for case B128D3 and $y^+\gt 40$ for case B128D2), contrary to the common understanding that the nonlinear term dominates throughout the channel in single-phase flows (Kim Reference Kim1989; Chang III et al. Reference Chang, Piomelli and Blake1999). Referring to its expression, the larger linear term may arise from the gradient of the mean density in combination with the mean velocity. This finding emphasises that the bubble-induced contribution to pressure fluctuations originates primarily from variable-density effects, involving not only density fluctuations but also the spatial gradient of the mean density.

These results indicate that the pressure fluctuations induced by bubbles originate primarily from two mechanisms: (i) a density-related effect, including the spatial variation of mean density and the interfacial jump of density, and (ii) a surface-tension-related effect, inducing a pressure difference across bubble interfaces. These combined effects constitute a fundamental distinction from single-phase flows, in which pressure fluctuations are predominantly governed by incompressible turbulence dynamics.

5.1. The method for calculating the wavenumber–frequency spectrum

The wall-pressure fluctuations are analysed in terms of the wavenumber–frequency spectrum, which is calculated using the method described by Choi & Moin (Reference Choi and Moin1990) and then by Yang & Yang (Reference Yang and Yang2022). In this subsection we briefly review the calculation method.

The sequence of the 3-D wall-pressure fluctuations, denoted by $p'(x, z, t)$ , is discretised into $M$ segments in time, each with a 50 % overlap with its adjacent segment. Each segment contains a number of successive snapshots with $N_t = 256$ for the single-phase case and $N_t = 500$ for bubbly cases. Subsequently, a Fourier transform is performed within each segment in both the streamwise and spanwise directions, as well as in time, to extract the Fourier components. Given that the data are discretised on a grid, the discrete Fourier transform is employed to calculate the Fourier modes $\widehat {p^\prime }( {k_x},{k_z},\omega )$ as

(5.1) \begin{eqnarray} \widehat {p^\prime }\left ( {{k_x},{k_z},\omega } \right ) & = & \frac {1}{{{N_x}{N_z}{N_t}{\sqrt {C_w}}}}\sum \limits _{l = 0}^{{N_x} - 1} {\sum \limits _{m = 0}^{{N_z} - 1} {\sum \limits _{n = 0}^{{N_t} - 1} w } } \left ( {\frac {{nT}}{{{N_t}}}} \right )p^\prime \left ( {\frac {{l{L_x}}}{{{N_x}}},{\mkern 1mu} \frac {{m{L_z}}}{{{N_z}}},{\mkern 1mu} \frac {{nT}}{{{N_t}}}} \right ) \nonumber \\ && \times \exp \left ( { - {\textrm {i}}\left ( {\frac {{l{L_x}}}{{{N_x}}}{k_x} + \frac {{m{L_z}}}{{{N_z}}}{k_z} + \frac {{nT}}{{{N_t}}}\omega } \right )} \right )\!, \end{eqnarray}

(5.2) \begin{equation} C_w=\frac {1}{{{N_t}}}\sum \limits _{n = 0}^{{N_t} - 1} \; {w^2}\left ( {{\mkern 1mu} \frac {{nT}}{{{N_t}}}} \right )\!, \end{equation}

where $k_x$ and $k_z$ are the streamwise and spanwise wavenumbers, $\omega$ is the frequency, $C_w$ is the window constant to keep the r.m.s. of $p^\prime$ unchanged after applying the standard Hanning window function $w(t)$ , $T = {N_t} \times {T_s}$ is the time duration of each segment, with $T_s$ being the sampling interval. The 3-D wavenumber–frequency spectrum of wall-pressure fluctuations is defined as

(5.3) \begin{equation} {\varPhi _{p^\prime p^\prime }}({k_x},{k_z},\omega ) = \frac {{\overline {\widehat {p^\prime }({k_x},{k_z},\omega ){{\widehat {p^\prime }}^*}({k_x},{k_z},\omega )} }}{{\varDelta {k_x}\varDelta {k_z}\varDelta \omega }}, \end{equation}

where the superscript $*$ denotes the complex conjugate, $\varDelta k_x=2\pi /L_x$ and $\varDelta k_z=2\pi /L_z$ are the resolution of the streamwise wavenumber and spanwise wavenumber, respectively, and $\varDelta \omega =2\pi /T$ is the frequency resolution. Here, the overbar symbol represents the mean value over the entire set of $M$ segments. The reduced-dimension spectra are calculated by integrating the 3-D wavenumber–frequency spectrum, for example, applying

(5.4) \begin{equation} {\varPhi _{p^\prime p^\prime }}({k_x},\omega ) = \int {{\varPhi _{p^\prime p^\prime }}({k_x},{k_z},\omega ){\mathrm{d}}{k_z}} \end{equation}

to obtain the 2-D spectrum and

(5.5) \begin{equation} {\varPhi _{p^\prime p^\prime }}\left ( {{k_x}} \right ) = \int _{ - \infty }^\infty {\int _{ - \infty }^\infty {{\varPhi _{p^\prime p^\prime }}} } \left ( {{k_x},{k_z},\omega } \right ){\mkern 1mu} {\mathrm{d}}{k_z}{\mkern 1mu} {\mathrm{d}}\omega \end{equation}

to obtain the 1-D spectrum. The validation of the wavenumber–frequency spectrum for single-phase flows is provided in Appendix A.

5.2. Two-dimensional spectra

An overview of the characteristics of wavenumber–frequency spectra for pressure fluctuations through three representative cases, B64D3, B128D3 and SP, is illustrated in figure 12. Among all spectra, particular attention is given to the streamwise wavenumber–frequency spectrum $\varPhi _{p^\prime p^\prime }(k_x, \omega )$ , as it captures the distribution characteristics of wall-pressure fluctuations across different frequency components along the flow direction. The ridge in the spectrum denotes the convection process of energy-containing wall-pressure fluctuations. Compared with the spectrum in case SP (figure 12 c), which exclusively exhibits a single ridge, it is evident that the spectra in bubbly cases retain the fundamental characteristics of turbulent flow in the liquid phase. Additionally, a new convective ridge emerges, as can be more clearly observed from the contour lines in the $(k_x,\omega )$ plane shown in figures 12(a) and 12(b). This feature, hereafter referred to as the bubble convection ridge, will later be shown to be associated with bubble-induced dynamics, distinct from the turbulence convection ridge in single-phase flows. The solid black line along the convection ridge represents the most energetic pressure fluctuation convection with its expression $\varPhi _{p^\prime p^\prime }(k_x^*, \omega ^*)$ , where $\omega ^* / k_x^* = U_{bc}$ with $U_{bc}$ denoting the bulk convection velocity. On the two vertical planes ( $k_x-$ and $\omega -$ planes), variations of spectral magnitude with respect to frequency $\omega$ at a fixed streamwise wavenumber $k_x^*$ and wavenumber $k_x$ at a fixed frequency $\omega ^*$ are illustrated. Projections of the convection ridge onto the two vertical coordinate planes indicate variations in the approximate maximum pressure fluctuation intensities as functions of wavenumber and frequency. For all cases, the maximum intensity decreases with increasing wavenumber and frequency.

Figure 12.

Three-dimensional profiles of wavenumber–frequency spectrum $\varPhi _{p^\prime p^\prime }( k_x,\omega )$ and the contours for case (a) B64D3, (b) B128D3, and (c) SP. The solid black line is the ridge line with its horizontal coordinates $U_{bc}=\omega /k_x$ . The projections of the ridge line onto the $k_x$ plane (dot–dashed black line) and the $\omega$ plane (dashed black line) are plotted. The spectrum for a fixed streamwise wavenumber $\varPhi _{p^\prime p^\prime }(k_x^*,\omega )$ is plotted as a dot–dashed line and for a fixed frequency, $\varPhi _{p^\prime p^\prime }(k_x,\omega ^*)$ is plotted as a dashed line.

Figure 12. Long description

Panel C: A 3D plot depicting the wavenumber-frequency spectrum profiles and contours for case SP. The horizontal axis represents the streamwise wavenumber kxh, the vertical axis represents the frequency omega h divided by Ub, and the third axis represents the logarithm of the spectrum. The solid black line indicates the ridge line, with its projections onto the kxh plane shown as a dot-dashed black line and onto the omega h divided by Ub plane shown as a dashed black line. The spectrum for a fixed streamwise wavenumber is plotted as a dot-dashed line, and for a fixed frequency, it is plotted as a dashed line. The plot shows a peak around kxh of 100 and omega h divided by Ub of 0, with the spectrum values decreasing as you move away from this peak.

For clarity, figure 13(a–d) presents the contours of the 2-D wavenumber–frequency spectrum $\varPhi _{p^\prime p^\prime } ( k_x,\omega )$ for bubbly cases. The general features observed in cases B64D2 and B128D2 are similar to those of case SP (figure 22), primarily exhibiting the turbulence convection ridge. In cases B64D3 and B128D3, a distinct secondary convection ridge, that is, the bubble convection ridge, emerges alongside the primary turbulence-related ridge. This feature is attributed to bubble-induced pressure fluctuations based on the following evidence. First, as the gas volume fraction increases, the greater concentration of bubbles near the wall leads to a more pronounced ridge with higher spectral magnitudes and broader bandwidth, consistent with intuitive physical expectations. Second, the bubble convection ridge exhibits a markedly steeper slope and extends toward higher wavenumbers and frequencies, clearly distinguishing it from the turbulence convection ridge.

Figure 13.

The 2-D spectrum in the streamwise wavenumber and frequency normalised as $\textrm {log}_{10}[\varPhi _{p^\prime p^\prime }(k_x,\omega )u_\tau /(\tau _w^2h^2)]$ for cases (a) B64D2, (b) B128D2, (c) B64D3, and (d) B128D3. The 2-D spectrum $\varPhi _{p^\prime H}(k_x,\omega )$ of case (e) B64D3 and (f) B128D3.

Figure 14.

The wavenumber–frequency spectrum of pressure fluctuations at fixed wavenumbers $\varPhi _{p^\prime p^\prime }(k_x^*,\omega )$ with $k_x^* h=27,35,50$ and 100.

Figure 14. Long description

Four line graphs depict the wavenumber-frequency spectrum of pressure fluctuations at fixed wavenumbers. Panel A: The line graph shows the spectrum for kxh equal to 27. The x-axis represents the non-dimensional frequency, ωh/U_b, ranging from -120 to 20. The y-axis represents the non-dimensional spectrum, Φp’p’(kx*, ω)uτ/(τw^2h^2), ranging from -12 to -4. Different colored lines represent different datasets: B64D2, B128D2, B64D3, B128D3, and SP. Panel B: The line graph shows the spectrum for kxh equal to 35. The axes and datasets are the same as in Panel A. Panel C: The line graph shows the spectrum for kxh equal to 50. The axes and datasets are the same as in Panel A. Panel D: The line graph shows the spectrum for kxh equal to 100. The axes and datasets are the same as in Panel A. Annotations indicate bubble-induced and turbulence-induced regions.

To further substantiate the bubble-induced nature of this secondary ridge, we employ a cross-correlation analysis to explicitly link the observed convective ridge in $\varPhi _{p^\prime p^\prime } ( k_x,\omega )$ for cases B64D3 and B128D3 to actual bubble motion. In order to calculate the cross-correlation function between the wall-pressure fluctuation and the distance of bubble from the wall, we define a function $H(x,z,t)$ to measure the height of the bubble surface such that a larger value corresponds to a bubble closer to the wall, taking the $y = 0$ wall as an example,

(5.6) \begin{equation} H(x,z,t)=2h-y_b(x,z,t), \end{equation}

where $y_b(x,z,t)$ denotes the $y$ coordinate of the lowest bubble interface. As bubbles move, the function $H$ and $y_b$ vary with spatial coordinates $(x, z)$ and time $t$ . If no bubble exists at a specific $(x, z)$ coordinate, $y_b$ is set to $2h$ , resulting in $H=0$ .

The 2-D cross-correlation function between the wall-pressure fluctuation $p^\prime$ and the bubble height function $H$ in the streamwise distance $\xi$ and time interval $\tau$ is defined as

(5.7) \begin{equation} {R_{p^\prime H}}\left ( {\xi ,\tau } \right ) = \frac {1}{{{N_z}}}\sum \limits _z {\left \langle {p^\prime \left ( {x,y=0,z,t} \right )H\left ( {x + \xi ,z,t + \tau } \right )} \right \rangle } , \end{equation}

where the angular brackets denote the ensemble average over all $x$ and $t$ . The normalised 2-D cross-correlation function $R_{p^\prime H}(\xi ,\tau )/R_{p^\prime H}(0,0)$ is converted into the spectrum $\varPhi _{p^\prime H} ( k_x,\omega )$ through a Fourier transform. Figures 13(e) and 13(f) present the pseudocolour plots of the spectrum $\varPhi _{p^\prime H} (k_x,\omega )$ for cases B64D3 and B128D3. For case B64D3, the spectrum $\varPhi _{p^\prime H} ( k_x,\omega )$ exhibits a narrow ridge whose shape and slope closely resemble those of the bubble convection ridge in $\varPhi _{p^\prime p^\prime } ( k_x,\omega )$ in figure 13(c), while the broader turbulence convection ridge is excluded through the cross-correlation procedure. For case B128D3, the ridge in $\varPhi _{p^\prime H} ( k_x,\omega )$ correspondingly becomes broader and exhibits higher intensity like the changes in $\varPhi _{p^\prime p^\prime } ( k_x,\omega )$ in figure 13(d). This supports the association of the secondary ridge in $\varPhi _{p^\prime p^\prime } ( k_x,\omega )$ with bubble-induced pressure fluctuations.

The contour in figure 13 presents the morphology of the bubble convection ridge and turbulence convection ridge, while figure 14 provides quantitative cross-sectional views of $\varPhi _{p^\prime p^\prime }(k_x, \omega )$ at fixed streamwise wavenumbers $k_x^*$ . Overall, the magnitude of the bubble convection ridge exceeds that of the turbulence ridge roughly at a certain wavenumber. The bubble convection ridge becomes more pronounced in the superconvective region (higher wavenumber or frequency region). In contrast, the ridge magnitudes in cases with a low gas volume fraction, i.e. B64D2 and B128D2, are lower than those in the SP case, with minor differences in the low wavenumber or frequency region but significant deviation in the superconvective region. The influence of bubbles is evident in two main aspects. First, a higher volume fraction leads to more bubbles towards the wall, resulting in stronger bubble-induced wall-pressure fluctuation, which is specifically reflected in the higher spectral magnitude and broader bandwidth observed in case B128D3 compared with case B64D3. Second, as shown in figures 14(a) and 14(c), the bubble ridge in B128D3 emerges at a lower wavenumber, surpassing the turbulence ridge at $k_x h \approx 27$ , whereas in B64D3, it is $k_x h \approx 50$ . These wavenumbers correspond to scales of $k_x D_b = \textit {O}(8{\sim} 15)$ , i.e. wavelengths of roughly 0.8 $D_b$ –0.4 $D_b$ . This indicates that bubble motion enhances pressure fluctuations at scales comparable to the bubble size, further supporting the interpretation that this spectral ridge originates from bubble-induced motions.

For the low gas volume fraction case, B64D2, the magnitude of the wavenumber–frequency spectrum of pressure fluctuation is found to be suppressed. This behaviour can be partly explained by the distribution of turbulent coherent structures over the channel wall. The $Q$ -criterion is utilised to identify vortices, with results presented in figure 15 for the cases SP, B64D2 and B128D3. Under an identical $Q$ value, the present numerical results show that introducing bubbles to the channel disrupts and reduces considerably the coherence of the longitudinal vortex structures near the wall. In the core region of the channel, bubble-induced turbulence dominates over the shear-induced turbulence. Here, the vortices in both the wake region and near the bubble surface are smaller in size but denser and more intense than those in the single-phase case. In addition, as shown by the visualisation of vortices generated on the bubble surfaces, the vortex layers are relatively thin. The interfacial vorticity is confined in a boundary layer of thickness $\delta /D_b \approx \textit{Re}_b^{-1/2}$ (Riboux, Legendre & Risso Reference Riboux, Legendre and Risso2013). Accordingly, the present mesh resolution of $D_b /\varDelta \geqslant 30$ captures the interfacial vorticity with reasonable accuracy for $\textit{Re}_b \sim 50$ , as discussed in Appendix D.

The anisotropic characteristics of flow structures can be qualitatively revealed according to the 2-D wavenumber spectrum $\varPhi _{p^\prime p^\prime }(k_x,k_z)$ , as shown in figure 16. For case SP, the isopleths are elongated in the spanwise direction, indicating that the flow structures contributing to wall-pressure fluctuations are more extended spanwise than streamwise – a feature also reported by Choi & Moin (Reference Choi and Moin1990). In the bubbly cases, the isopleth shapes in low gas volume fraction cases (B64D2 and B128D2) remain similar to those in the single-phase flow. In contrast, cases B64D3 and B128D3 exhibit noticeable changes in the isopleth patterns: the low-magnitude contours bulge toward the $k_x$ direction at low $k_z$ wavenumbers, indicating that bubbles amplify the spectral components associated with large spanwise scales and small streamwise scales.

Figure 15.

Isosurfaces of the $Q$ -criterion ( $Q = 7 \times 10^4$ ), coloured in teal, are shown for (a) SP, (b) B64D2, (c) B128D3, with a zoomed-out view of two adjacent bubbles interacting with wall-generated vortical structures. The red plane indicates the plane at $y^+=30$ . Translucent spheres represent the instantaneous bubble positions.

Figure 16.

The 2-D wavenumber spectrum ${\textrm {log}_{10}}[\varPhi _{p^\prime p^\prime }(k_x,k_z)/(\tau _w^2h^2)]$ for cases: (a) B64D2, (b) B128D2, (c) B64D3, (d) B128D3, and (e) SP.

5.3. Convection velocity of the 2-D spectrum

Convection velocity is a key parameter that can be extracted from the convection ridge in the wavenumber–frequency spectrum. It represents the propagation speed of pressure fluctuations correlated with coherent motions. In practical applications, it can also provide a mapping between spatial and temporal quantities. According to the early experimental data (Willmarth & Wooldridge Reference Willmarth and Wooldridge1962), the convection velocity of the maximal energy-contained eddies ranges from 0.56 to 0.83 times the free-stream velocity. Since Taylor’s hypothesis is not applicable in the near-wall region, the convection velocity, which should depend on the scale of flow structures, is not constant. The distinct characteristics of the bubble-induced convection ridge, as compared with the turbulence-induced ridge, indicate that multiple types of coherent motions coexist. In this section we investigate the scale-dependent convection velocity of wall-pressure fluctuation.

The streamwise-wavenumber-dependent convection velocity can be calculated by

(5.8) \begin{equation} {U_c}\!\left ( {{k_x}} \right ) = -\frac {{\int {\omega \varPhi ({k_x},\omega ){\mkern 1mu} } {\mathrm{d}}\omega }}{{{k_x}\int {\varPhi ({k_x},\omega ){\mkern 1mu} } {\mathrm{d}}\omega }}, \end{equation}

which is equivalent to searching the position of the centre of gravity in the frequency spectrum at a given wavenumber $k_x^*$ (del Álamo & Jiménez Reference del Álamo and Jiménez2009). As shown in figure 17(a), the velocities in case SP decrease with increasing wavenumber, which agrees with previous studies (e.g. Choi & Moin (Reference Choi and Moin1990)) and align with the empirical curve in Panton & Linebarger (Reference Panton and Linebarger1974). This variation is attributed to the higher convection velocity of larger-scale flow structures originating in the outer layer, while smaller-scale structures near the wall exhibit a lower velocity. For case B64D2, the shape of the 2-D spectrum and the resulting convection velocity ${U_c} (k_x )$ are similar to those in the single-phase flow, indicating that dilute bubbles have limited influence on the convection velocity of wall-pressure fluctuations. In contrast, for cases B128D2, B64D3 and B128D3, the presence of the bubble convection ridge shifts the centre of gravity of $\varPhi _{p^\prime p^\prime }(k_x^*, \omega )$ toward higher frequencies, leading to a noticeable increase in ${U_c} (k_x )$ beyond a certain wavenumber. The higher magnitude of ${U_c} (k_x )$ in case B128D3 compared with B64D3 corresponds to the stronger bubble convection ridge observed in the high gas volume fraction case. The higher convection velocity combines the effects of the bubble motions in the channel, with smaller-scale flow structures moving faster than the turbulent structure near the wall.

Alternatively, the convection velocity that depends on the streamwise wavenumber can also be determined by searching for a maximum along $\omega$ for fixed $k_x^*$ (del Álamo & Jiménez Reference del Álamo and Jiménez2009) as

(5.9) \begin{equation} {U_c}{\left ( {{k_x}} \right )_{\textit{peak}}} = - \frac {{{\omega _c}\left ( {{k_x}} \right )}}{{{k_x}}},\quad {\left . {\frac {{\partial \varPhi ({k_x},\omega )}}{{\partial \omega }}} \right |_{\omega = {\omega _c}}} = 0 . \end{equation}

This definition extracts the convection velocity closely related to the bubble motion since the peaks are located within the bubble convection ridge. The results are illustrated in figure 17(b), which is not as smooth as figure 17(a). Searching for the peak numerically, rather than integration, may lead to data scattering. In cases B128D2, B64D3 and B128D3, a distinct and nearly constant convection velocity is observed. For case B64D3, the convection velocity within this range is $0.917U_b$ , larger than that of case B128D3, $0.862U_b$ . These values reflect that the transport of near-wall bubbles contributes more to wall-pressure fluctuations. Therefore, case B128D3, which contains more bubbles closer to the wall, exhibits a smaller peak convection velocity. It is noteworthy that the constant convection velocity in $U_c(k_x)_{\textit{peak}}$ begins at $k_x h \gt 50$ for case B64D3 and $k_x h \gt 27$ for case B128D3. These values can be regarded as the thresholds beyond which the intensity of the bubble convection ridge overwhelms that of the turbulence convection ridge, consistent with the results in figure 14. In comparison, the transition points in the ${U_c}(k_x)$ profiles (figure 17 a) mark the onset of the bubble convection ridge, occurring at $k_x h \approx 35$ for case B64D3 and $k_x h \approx 20$ for case B128D3, respectively. The gaps between the transition points defined by different convection velocities represent the wavenumber range over which the bubble convection ridge develops.

Figure 17.

Comparisons on the different convection velocities versus wavenumbers: (a) $U_c(k_x)$ and (b) $U_c(k_x)_{\textit{peak}}$ .

Finally, the bulk convection velocity is calculated by

(5.10) \begin{equation} {U}_{bc} = - \frac { {\iint {{k_x}} } \omega \varPhi ({k_x},\omega ){\mkern 1mu} {\mathrm{d}}\omega {\mkern 1mu} {\mathrm{d}}{k_x}}{ {\iint {k_x^2} } \varPhi ({k_x},\omega ){\mkern 1mu} {\mathrm{d}}\omega {\mkern 1mu} {\mathrm{d}}{k_x}} \end{equation}

to characterise the overall propagation speed, as listed in table 3. Bubbly flows enhance the propagation speed of wall-pressure fluctuation and higher volume fraction cases tend to exhibit larger values of $U_{bc}$ overall. The effective sources of the pressure field can be localised by relating the convection velocity to the mean velocity (Schewe Reference Schewe1983). As indicated by the markers in figure 6, the larger convection velocities observed in the bubbly cases suggest that the dominant sources are located farther from the wall. Referring to the magnitudes of the pressure fluctuations across the channel, the motions of the bubbles induce the dominant pressure fluctuation. These sources, compared with the single-phase case, are shifted away from the wall, leading to larger convection velocities associated with the bubbly flows.

Table 3.

The bulk convection velocities normalised by the bulk velocity $U_{bc} / U_b$ for all cases.

5.4. One-dimensional spectrum

The 1-D spectra, including the streamwise-wavenumber spectrum $\varPhi _{p^\prime p^\prime }(k_x)$ , spanwise-wavenumber spectrum $\varPhi _{p^\prime p^\prime }(k_z)$ and frequency spectrum $\varPhi _{p^\prime p^\prime }(\omega )$ of wall-pressure fluctuations are derived by integrating the spectrum function (5.3) over the remaining two variables. Figure 18 illustrates the spectra of the bubbly and single-phase cases. The spectral magnitudes are normalised using their respective wall units for each individual case. The results obtained for case SP are consistent with the DNS results reported by Yang & Yang (Reference Yang and Yang2022). Overall, the differences between the bubbly flow and single-phase flow cases are relatively minor. The spectra at low streamwise-wavenumber range ( $k_x h \lt 10$ ) exhibit a non-monotonic behaviour, rising first before declining, for both bubbly and single-phase cases, while the spanwise-wavenumber spectra remain consistent with expected trends. The frequency spectra show oscillations in the low-frequency range ( $\omega h / U_b \lt 10$ ) for the bubbly cases, reflecting the characteristic time scales of bubble motion. The scaling laws observed for the streamwise-wavenumber spectrum $\varPhi _{p^\prime p^\prime }(k_x)$ agree with those presented by Choi & Moin (Reference Choi and Moin1990), specifically the slopes of −1 and −5, which were obtained from the DNS database. The slopes of −0.7, −7/3 and −5 in the $\varPhi _{p^\prime p^\prime }(\omega )$ spectrum are in agreement with the spectral characteristics outlined in Hwang, Bonness & Hambric (Reference Hwang, Bonness and Hambric2009). It should be noted, however, that in case B128D3 the spectral decay becomes slightly shallower and does not strictly follow the −5 scaling, due to the influence of the bubble-induced convective ridge.

Figure 18.

One-dimensional spectra of wall-pressure fluctuations with respect to (a) streamwise wavenumber $k_x$ , (b) spanwise wavenumber $k_z$ , and (c) frequency $\omega$ . The vertical dashed lines denote $k_xD_b$ and $k_z D_b$ , respectively, with $D_b=2\ \text{and}\ 3$ mm.

Figure 19.

Wall-normal distribution of the Green’s function $G(y=0, y^\prime , k)$ at different $kh$ values.

Although the newly identified bubble convection ridge in the 2-D spectra of bubbly flows clearly distinguishes from that of single-phase flow (see § 5.2), such distinct differences are insignificant in the 1-D spectra, as shown in figure 18. The integration from the 2-D spectrum to 1-D spectrum may conceal certain physical characteristics.

The diminished discrepancy in 1-D spectra can be interpreted from another perspective. The wall-pressure fluctuation can be expressed as the integral of the Fourier-transformed source terms over the entire channel height, weighted by the Green’s function of the Poisson equation (Kim Reference Kim1989; Anantharamu & Mahesh Reference Anantharamu and Mahesh2020). The Green’s function $G(y=0, y^\prime , k)$ quantifies the contribution of a source at height $y^\prime$ with wavenumber $k = \sqrt {k_x^2 + k_z^2}$ , imposing an exponential decay of the form ${\rm e}^{-k y^\prime }$ (Kraichnan Reference Kraichnan1956). As shown in figure 19, the influence on the wall-pressure fluctuations decreases as the source moves away from the wall, and depends on the dimensionless wavenumber $kh$ . Although the source terms in bubbly flows are stronger than those in the single-phase flow within the channel core, their contributions to the wall-pressure fluctuations are substantially attenuated with the increasing distance from the wall owing to the associated exponential decay rate. This explains why the pressure fluctuations exhibit significant discrepancies in the channel core, yet the differences appear less pronounced in the 2-D spectra, remain insignificant in the 1-D spectra and become minimal in the wall value of $p^\prime _{\textit{rms}}$ . Nevertheless, the 2-D spectra are still able to capture the essential characteristics of wall-pressure fluctuations that are induced by bubbles.