Experiments often use gap width averaged velocities to study statistical properties of the flow. Following the procedure outlined in [20,21] and assuming density ϱ to be constant along wall-normal direction, and starting from the NSHS equations, we get the following equations for the gap averaged indicator function c and horizontal components of the velocity ² : ∂ t + u ⋅ ∇ c = 0 ,   ∇ ⋅ u = 0 , a n d ρ ∂ t + u ⋅ ∇ u = − ∇ ⋅ ρ U ′ U ′ ̄ + ∇ ⋅ 2 μ ̄ S + F ̄ d − ∇ P ̄ + F ̄ g + F ̄ σ . (2) Here, ( . ) ̄ ≡ ( 1 / H ) ∫ 0 H ( . ) d y denotes gap averaging, ∇ ≡ (∂ _x , ∂ _z ), u ≡ (u _x (x, z), u _z (x, z)) is the gap averaged velocity field with u x = U ̄ x and u z = U ̄ z , U′(x) = (U _x − u _x , U _z − u _z ) are the three-dimensional residual velocity fluctuations, P ̄ ( x , z ) is the gap-averaged pressure field, S = ∇ u + ∇ u T is the gap averaged strain-rate tensor, ρ ≡ ϱ ̄ is the gap averaged density field, F ̄ g = [ ρ a − ρ ] g z ̂ is the buoyancy force, and F ̄ σ is the surface tension force. The viscous dissipation contributes in two parts: (a) small-scale dissipation ∇ ⋅ 2 μ ̄ S , and (b) viscous drag due to walls F ̄ d = [ μ ( ∇ U + ∇ U T ) ⋅ y ̂ ] 0 H . Previous studies [13,17,18] have argued that the gap-averaged fields appropriately describe the dynamics of bubbly Hele-Shaw flows provided H/d ≤ 1/2. In particular, they model the dynamics of an isolated bubble using two-dimensional Navier-Stokes equations with a linear drag which has similar structure to (2). We will revisit this point later in Section 3.3.

We use a second-order finite-volume solver PARIS [22] to simulate NSHS (1). For bubble tracking and updating the indicator function PARIS employs a front-tracking method, and the time marching is performed using the first order Euler method.

We consider a cuboid of breadth L _y = H, and with equal length and height (L _x = L _z = L) [see Figure 1]. The simulation domain (L _x , L _y , L _z ) is discretized with (N _x , N _y , N _z ) equispaced collocation points. We use periodic boundary conditions in the x and z directions, and impose no-slip velocity boundary U = 0 at the walls (y = 0 and y = H) [23]. We place N _b bubbles in random positions and initialize each one as an ellipsoid of volume V = 4.73 × 10³ (mono-disperse suspension). The front tracking module of the PARIS solver [22] employs a symmetry boundary condition for the surface tension force at the walls such that the net force at the wall would be tangent to the wall. This ensures a thin layer of fluid (approximately O(0.13)δy, where δy = H/N _y ) between the bubble and the wall. Our grid resolution is comparable (or higher) to a previous numerical study on Hele-Shaw with bubbles [23]. The bubbles are allowed to relax in the absence of gravity until they achieve the equilibrium pan-cake-like configuration [23] with diameter d = 2H.

In Table 1 we summarize the parameters used in our simulations. We choose parameters such that the dimensionless numbers (Ga, Bo, H/d, and ϕ) are comparable to experiments [15,17,18]. We simulate low At = 0.08 and high At = 0.9, and verify that the spectral properties are insensitive to density contrast [6,11,12].

From (2), we obtain the following balance equation for the gap-averaged kinetic energy E ∂ t ⟨ ρ u 2 2 ⟩ ⏟ E = − 2 ⟨ μ ̄ S : S ⟩ ⏟ ϵ μ + ⟨ ρ a − ρ u z g ⟩ ⏟ ϵ inj + ⟨ F ̄ σ ⋅ u ⟩ ⏟ ϵ σ + ⟨ F ̄ d ⋅ u ⟩ ⏟ ϵ d , (3) where ϵ _μ is the gap-averaged viscous energy dissipation, ϵ _d is the dissipation due to drag, ϵ _inj is the gap-averaged energy injected due to buoyancy, ϵ _σ is the contribution due to the surface tension, and the angular brackets denote spatial averaging.

In Figure 3, we plot the time-evolution of the kinetic energy E and observe that a statistically steady state is achieved for t > 0.8τ _s . Furthermore, in Table 2 we show that the energy injected by buoyancy is primarily balanced by the dissipation due to drag in the steady state (∂ _t E ≈ 0).

The energy spectrum and co-spectra for the gap-averaged velocity field are defined as: E k ≡ ∑ k − 1 / 2 < m < k + 1 / 2 | u ̂ m | 2 , E ρ u k ≡ ∑ k − 1 / 2 < m < k + 1 / 2 R ρ u ̂ − m u ̂ m . Here, ( ⋅ ) ̂ denotes the Fourier transformed fields.

In Figure 4, we plot the energy spectra E(k) and cospectra E ^ρu (k) for our simulations H 1 − H 3 ³ . From the plots, we can identify different scaling regimes: (a) For k ≪ k _d we observe E(k) ∼ k, where k _d is the wavenumber corresponding to the bubble diameter; (b) Around k ∼ k _d , we find a short E(k) ∼ k ⁻³ scaling regime followed by a steeper decay of the spectrum E(k) ∼ k ⁻⁵.

The k ⁻³ scaling range observed in our simulaiton is consistent with earlier experiments that also observe an intermediate k ⁻³ scaling subrange for 0.2⪅k/k _d ⪅1 [15]. We verify this by overlaying the energy spectrum obtained in [15] over our data in Figure 4.

Risso [24] argues that the k ⁻³ scaling could be modelled as a signal consisting of a sum of localized random bursts. Although this explanation is consistent with Figure 2, it does not highlight the underlying mechanisms that generate the observed scaling. Lance and Bataille [2] take an alternate viewpoint and argue that the balance of energy production and viscous dissipation leads to the k ⁻³ scaling.

The scaling of the energy spectrum we observe differs from earlier studies on two-dimensional unbounded flows [10,12] at comparable Ga. They find an inverse energy cascade with E(k) ∼ k ^−5/3 for k < k _d and a E(k) ∼ k ⁻³ scaling for k > k _d due to the balance of energy injected by surface tension with viscous dissipation.

In what follows, we present an energy budget analysis to explain the observed scaling of the energy spectrum.

In this section we investigate whether two-dimensional Navier-Stokes equations with a linear drag coefficient (5) are able to model the confined bubbly flows In the following, we assume all the fields are two-dimensional and for comparison with the gap-averaged quantities, we choose the same symbols. D t c = 0 , a n d ∇ ⋅ u = 0 , ρ c D t u = ∇ ⋅ 2 μ c S − ∇ P + F g + F σ − α u . (5) Earlier studies [13,17,18,28] have used NSD Eq. 5 to investigate the dynamics of an isolated bubble in a Hele-Shaw setup with H/d ≤ 0.5 and found the results to be consistent with the experiments. In the following, we use NSD equations to study bubbly flows.

We perform direct numerical simulation of NSD equations on a square domain of area L ² and discretize it with 2048² equi-spaced points. The bubbles are initialized as circles of diameter d = 24 and all the parameters of the simulation are identical to our run H 1 and we fix the drag coefficient α = 0.04. Our choice for the value of α is motivated by Figure 5. We use a front-tracking-pseudo-spectral method to evolve (5). For details of the numerical scheme, we refer the reader to Ramadugu et al. [12]. Below we discuss the statistical properties of the flow in the steady state.

In Figure 7A, we plot the bubble configuration and the flow streamlines. Clearly the large scale flow properties resemble those observed for the NSHS simulation. The flow disturbances are localized in the vicinity of the bubbles and we also observe horizontal alignment of bubbles.

The plot in Figure 7B shows a comparison of the gap-averaged energy spectrum E(k) obtained from the NSHS equation with that obtained from NSD Eq. 5. We find that the energy spectrum are nearly identical for k < k _d , E(k) ∼ k. However, discrepancies are observed for k > k _d , in contrast to E(k) ∼ k ⁻⁵ for the NSHS simulations we find E(k) ∼ k ⁻³ for the NSD simulations. As pointed out in earlier section, the oscillation in the energy spectrum appear with a period k _d .

Using (5), and ignoring the inertial contributions, we obtain the following energy balance F k + T σ k = ν k 2 E k ⏟ D + α E k ⏟ D . (6)

In Figure 8, we plot the contribution of different terms in (6) towards energy balance. For k < k _d , similar to NSHS, we observe that energy injected by buoyancy is balanced by the linear drag. However, a different balance appears for k > k _d . In contrast to NSHS, a dominant balance is observed in the NSD equations. The energy transfer by surface tension to small scales balances viscous dissipation leading to the observed E(k) ∼ k ⁻³ scaling in the energy spectrum. Similar small-scale balance has also been reported in earlier two-dimensional unbounded bubbly flow simulations [12]. Therefore, we conclude that although the NSD model captures the large scale dynamics of the Hele-Shaw flow (NSHS), it is unable to correctly capture the small scale physics.