Theoretical elucidation of effect of drag force and translation of bubble on weakly nonlinear pressure waves in bubbly flows

Many forces act on gas bubbles, such as drag, lift, gravity, and virtual mass force. The drag force acting on translational bubbles, in particular, is one of the most important factors in terms of the dynamics of single- or multi-bubble flows. Many theoretical studies and subsequent numerical and experimental works on the drag force and bubble dynamics have been carried out.^1–6 Furthermore, the translation and drag forces are quite important factors when assuming a high-velocity water flow accompanied by cavitation in hydraulic machinery.

In cavitating water flows, pressure variation is always generated by bubble oscillations and evolves into pressure waves (non-void waves⁷). As in nonlinear acoustics or nonlinear wave theory for pure (non-bubbly) water flows, the pressure wave evolves into a shock wave owing to the competition between the nonlinear effect and dissipation effect of the waves. Contrarily, for bubbly flows, the volumetric oscillations of bubbles induce a dispersion effect of waves. Owing to the competition between the nonlinear and dispersion effects, the pressure wave evolves into a stable wave, the so-called (acoustic) soliton. The shock wave and the soliton are exact solutions of the Burgers and Korteweg–de Vries (KdV) equations,⁸ respectively. Hence, it is important to explore the relative ratios of the nonlinear, dissipation, and dispersion effects, because the pressure wave in bubbly flows may evolve into a shock wave or a soliton,⁹ which have quite different properties. In the framework of weakly nonlinear (i.e., finite but small amplitude¹⁰) waves in bubbly liquids, the Korteweg–de Vries–Burgers (KdVB) equation¹¹ is one of the most famous nonlinear wave equations. As the KdVB equation comprises a linear combination of the nonlinear, dissipation, and dispersion terms [see (55)], the relative magnitude of these effects determines whether the pressure waves evolve into the shock wave or the soliton. Although estimating the magnitude of these three effects is essential to predict the evolution of the pressure wave, it cannot be obtained directly from experiments or numerical analysis of the governing equations. Therefore, the theoretical derivation of weakly nonlinear wave equations, such as the KdVB equation, is an effective method of estimating the relative strength of the above-mentioned effects.

Since the pioneering work of van Wijngaarden,¹² many studies have derived the KdVB equation from the gas–liquid mixture model¹³ to describe pressure waves in bubbly liquids. Furthermore, Kuznetsov et al. showed that a waveform obtained by the KdVB equation agrees with the experimentally observed propagation of weakly nonlinear pressure waves in a bubbly flow.¹⁴ However, all existing theoretical studies have ignored the drag force acting on the bubbles. Furthermore, only some studies^15,16 incorporated the translation of bubbles and the initial velocity of bubbly liquids. This may stem from the preconception that the effect of the non-oscillating components (i.e., translation and drag forces) on the oscillating components (i.e., bubble oscillations and pressure waves) is negligible. Because momentum transport across the bubble–liquid interface should be formulated to incorporate the drag force, such complex basic equations (e.g., two-fluid model equations) are required to resolve the weakly nonlinear (or linear) wave problem. Recently, based on the two-fluid model equations,¹⁷ our group theoretically investigated the linear^17,18 and weakly nonlinear^19,20 waves and derived the KdVB equation¹⁹ under the assumption that the translation and drag forces are negligible. Later, the initial nonuniform flow velocities were also incorporated.²¹ On the other hand, our group numerically solved the original KdVB equation¹⁹ and studied the temporal evolution of waveforms taking into account the nonlinear, dissipation, and dispersion effect; however, the translation and drag forces were not considered.²² 

In this study, we employ the two-fluid model¹⁷ that accounts for the interfacial momentum transport and introduce the translation and drag forces. Unlike our previous work in which the initial flow was assumed to be at rest,¹⁹ herein, at the initial state, the gas and liquid phases have velocities in this study as our preceding study.²¹ The remainder of this paper is organized as follows: In Sec. II, we introduce the basic equations of the two-fluid model and bubble dynamics, including the drag force and translation, respectively. In Sec. III, we derive two types of KdVB equations and point out that the translation of bubbles increases the nonlinearity, while the drag force affects both the nonlinearity and dissipation. In Sec. IV, the evolution of waveforms is calculated using numerical analysis, and the effects of the translation and drag force are discussed quantitatively. Then, the causes of dissipation are decomposed and compared using long-period calculation. Finally, Sec. V concludes this paper.

This paper theoretically investigates the weakly nonlinear (i.e., finite but small amplitude) propagation of one-dimensional (plane) pressure progressive waves in a flowing compressible water. The water containing distributed small spherical gas bubbles is shown in Fig. 1. Initially, the gas and liquid phases flow with a constant velocity. We newly introduce drag force and translation to the bubble dynamics. The bubbles do not coalesce, break up, extinct, or appear. For simplicity, the gas viscosity, Reynolds stress, thermal conduction^23,24 of the gas and liquid, and phase change and mass transport across the bubble–liquid interface are ignored. We assume a laminar flow and do not use turbulent models.

To introduce the drag force in the interfacial momentum transport, we first utilize the conservation laws of mass and momentum for the gas and liquid phases based on a two-fluid model,¹⁷ 

where t^* is the time, x^* is the space coordinate, α is the void fraction (0 < α < 1), ρ^* is the density, u^* is the velocity, p^* is the pressure, and P^* is the averaged liquid pressure on the bubble–liquid interface. The subscripts G and L denote the volume averaged variables in the gas and liquid phases, respectively, and the superscript * denotes a dimensional quantity. The following model for the virtual mass force¹⁸ is introduced as the interfacial momentum transport F^*:

where β₁, β₂, and β₃ are constants and may be set as 1/2 for the spherical bubble. The Lagrange derivatives D_G/Dt^* and D_L/Dt^* are defined by

Furthermore, we introduce a model for the drag force term D^* for spherical bubbles,

where R^* is the radius of a representative bubble and C_D is the drag coefficient for a single spherical bubble (see Sec. II D). The gas viscosity and Reynolds stress induced by volume averaging are ignored. However, for simplicity, we use a Stokes-type drag force. We intend to validate these assumptions in a forthcoming paper.

The bubble dynamics is described as

where $cL0*$ is the speed of sound in pure water. We introduced the translation of bubbles^15,16 in the third term of the right-hand side of (8) into the original Keller equation.²⁵ Hence, the translation and volumetric oscillations are described by (8) as a linear combination.

To complete the set of (1)–(4) and (8), the polytropic equation of state for gas, the Tait equation of state for liquid, the mass conservation law of gas inside the bubbles, and the balance of normal stresses across the bubble–liquid interface are introduced,

where γ is the polytropic exponent, n is a material constant (e.g., n =  7.15 for water), σ^* is the surface tension, and μ^* is the liquid viscosity (see Sec. II D). It should be noted that the effect of liquid viscosity is considered only at the bubble–liquid interface. The physical quantities at the initial state are signified by the subscript 0, and they are constants.

Based on the method of multiple scales,¹⁰ four multiple scales¹⁰ are introduced as extended independent variables by assuming the finite but small nondimensional wave amplitude ε (≪1),

where the nondimensional independent variables are defined by t = t^*/T^* and x = x^*/L^*; T^* is the typical period and L^* is the typical wavelength. Here, the subscripts 0 and 1 correspond to the near and far fields,¹⁰ respectively.

The dependent variables are nondimensionalized and are expanded in power series of ε,

where U^* (= L^*/T^*) is the typical propagation speed and the initial nondimensional pressures of O(1) are defined by p_G0 = $pG0*/(ρL0*U*2)$ and p_L0 = $pL0*/(ρL0*U*2)$⁠. Furthermore, the ratio of initial densities of gas and liquid phases is negligibly small.¹⁹ 

Similar to our previous study,¹⁹ the set of sizes of the three non-dimensional ratios appropriate to the low-frequency long wave is also determined using ε,

where V, Δ, and Ω are constants of O(1), and the eigenfrequency of a single bubble $ωB*$ is given by

Equations (20)–(22) represent that the speed of sound in bubbly flows is much smaller than that in pure water, the initial bubble radius is much shorter than the typical wavelength, and the incident frequency of waves is much lower than the eigenfrequency of a single bubble, respectively.

The nondimensional liquid viscosity μ is defined by

where cases A and B correspond to low liquid viscosity (e.g., $R0*≳500 μm$⁠) and high liquid viscosity (e.g., $R0*≲10 μm$⁠), respectively.

In this paper, the drag coefficient is assumed to be $O(ε)$ and is defined as

where Λ is the constant of O(1). In case A, C_D depends on the Reynolds number Re (C_D = 16/Re).¹ 

This section focuses on the derivation of the two types of KdVB equations and theoretically discusses their physical meaning.

Substituting (13)–(25) into (1)–(4) and (8), the set of linear equations can be derived from the leading-order of approximation with the aid of (9)–(12),

Combining (26)–(30) results in the linear wave equation for the first-order variation of the bubble radius, R₁,

where v_p is the phase velocity given by

The definition of the linear Lagrange derivative D/Dt₀ is

It should be noted that the initial velocities of both phases are assumed to be the same (u_G0 = u_L0 ≡ u₀) for simplicity; however, the perturbations of each velocity are not the same (u_G1 ≠ u_L1). Setting v_p = 1 gives the explicit form of U^* as

By focusing on the right-running wave (i.e., introducing a moving coordinate $x0−vpt0$⁠), α₁, u_G1, u_L1, and p_L1 can be expressed in terms of R₁,

with

The effects of the translation and drag force are not evident in the near field. Therefore, in the absence of an initial flow velocity (i.e., u_G0 = u_L0 = u₀ = 0), D/Dt₀ becomes ∂/∂t₀, and the present result coincides with the results of our previous work.¹⁹ 

As in the case of O(ε), the following set of inhomogeneous equations of O(ε²) is derived:

The inhomogeneous terms K_3A, K_4A, and K_5A are defined as

where K_i (i = 1, 2, 3, 4, 5) are the original inhomogeneous terms (see Appendix 1 in Ref. 19) and $K5′$ is the inhomogeneous term that ignores the effect of liquid viscosity from $K5$⁠.

Then, (43)–(47) are combined into a single inhomogeneous equation

with

From the solvability condition for (51), K = 0 is required.¹⁹ From (13), the original independent variables x and t are restored,

Finally, we obtain the KdVB-A equation

through the variable transform

where Π₀ is the original advection coefficient [see Eq. (52) in Ref. 19] and the other constant coefficients are

Here, Π₁ is the nonlinear coefficient (see the explicit form in Appendix 1 of Ref. 19) and Π₃ is the dispersion coefficient. These coefficients are the same as the original ones.¹⁹ Furthermore, Π_1A is the nonlinear coefficient, Π_2A is the dissipation coefficient owing to acoustic radiation, and Π_4A is the dissipation coefficient owing to the drag force. For the preceding case without translation,²⁶ the nonlinear coefficient Π_1A reduces to Π₁, and the other coefficients remain unchanged.

In case B, the following set of inhomogeneous equations of O(ε²) can be derived by introducing the new scaling Λ. The mass conservation equations are the same as (43) and (44). However, $K3A$ in (45), $K4A$ in (46), and $K5A$ in (47) become $K3B$⁠, $K4B$⁠, and $K5B$⁠, respectively:

Finally, we obtain the KdVB-B equation as in the case of KdVB-A equation,

Here, Π₀, Π₂, and Π₃ are the same as the original coefficients,¹⁹ Π_1A is the nonlinear coefficient that is the same as that of the KdVB-A equation, and Π_4B is the nonlinear dissipation coefficient owing to the drag force, which is obtained as follows:

The second term on the right-hand side of (57) describes the effect of bubble translation, and the translation contributes to the nonlinear effect of waves. Figure 2 shows the dependence of the absolute value of the nonlinear coefficient on the initial void fraction α₀. It can be seen that the incorporation of bubble translation results in an increase in the absolute value of the nonlinear coefficient Π_1A and thus the increase of the nonlinear effect.

Then, we focus on Π_4A and Π_4B, the coefficients related to the drag force. Figures 3 and 4 show the dependence of Π_4A and Π_4B on the initial void fraction α₀, respectively. As can be seen, both Π_4A and Π_4B increase with increasing α₀ and with decreasing $R0*$⁠.

The two types of KdVB equations [(55) and (64)] may be difficult to solve analytically, thus detracting to some extent from the utility of the KdVB formalism. The most important effect is that the drag force and translation affect the dissipation and nonlinearity, respectively, as discussed below.

First, for drag force, by assuming the solution of R₁ to the KdVB equation (55) as a sinusoidal wave, we establish that the drag force (i.e., Π_4A) affects the wave dissipation effect. From an intuitive point of view, we also consider the physical reason for the impact of the drag force on dissipation. We consider that the drag force transports the momentum across the bubble–liquid interface against the flow. It is well known that the drag force contributes to the dissipation of the flow, not the pressure wave. In this study, it is newly found that the drag force also contributes to the dissipation of waves from the physico-mathematical point of view. Since this is only a theoretical prediction and is not justified experimentally owing to the complete absence of corresponding experiments and direct numerical simulations, our theory should be verified experimentally.

Second, for translation, the results of a previous numerical study²⁷ indicated that bubble translation (bubble slip) slightly affects the waveforms. Herein, we demonstrated that bubble translation affects nonlinearity from the physico-mathematical viewpoint.

In this section, we discuss the effect of the drag force and bubble translation based on the KdVB-A equation (55).

First, (55) is numerically solved via the split-step (spectrum) Fourier–Galerkin method under a periodic boundary condition; the detailed scheme is described in our previous study.²² Figure 5 shows the temporal evolution of the waveform for the case of α₀ = 0.05 and $R0*=500 μm$⁠; the blue, red, and black curves represent the solutions to the present KdVB-A equation (55), the KdVB equation with only translation, and the original KdVB equation,¹⁹ respectively. Note that the liquid viscosity in the original KdVB equation is dropped from the dissipation coefficient for consistency in the comparison.

Consequently, the tendency of the temporal evolution of the waveform remains unchanged qualitatively between the present and previous cases. The amplitude of the present case is, however, slightly smaller than that of the previous case. As the drag force contributes to the wave dissipation, the consideration of the drag force decreases the wave amplitude (i.e., wave attenuation). On the other hand, the nonlinear effect increases owing to bubble translation.

This subsection discusses the most important points in this study.

The fifth term on the left-hand side of (55), i.e., Π_4AR₁, as a new term, contributes to wave dissipation.²⁶ The liquid viscosity does not affect the dissipation coefficient, as shown in (58) for case A, because the nondimensional size of liquid viscosity [see (24)] is changed from our previous study.¹⁹ Hence, Π_2A is due to only the acoustic radiation originating from the oscillating bubbles in a compressible liquid. The dissipation effect due to the liquid viscosity then appears as an approximation of O(ε³), that is, at a region very far from the sound source in this situation.

Therefore, there are two types of dissipation effects, Π_2A and Π_4A, with different formulas. As it is impossible to know the temporal evolution of the two types of dissipation by only theoretical analysis, we use numerical analysis. Here, we define the magnitude of dissipation as

where the index $i$ represents the computational mesh of $ξ$ and numerically obtained values of R₁ in Sec. IV A are used. Equation (66) represents the spatial integration of the dissipation terms. Figures 6 and 7 show the temporal evolutions of Γ_A and Γ_D for various α₀ and $R0*$⁠, respectively. First, while Γ_A decreases with an oscillation (Fig. 6), Γ_D decreases almost monotonically (Fig. 7). Furthermore, Γ_D increases with increasing α₀ and decreasing $R0*$⁠; however, Γ_A is almost independent of $R0*$⁠. Finally, a frequency of Γ_A decreases with increasing α₀.

We have theoretically and numerically examined the weakly nonlinear propagation of pressure waves in bubbly flows, especially focusing on the effects of the translation and drag force acting on the bubbles. The main results are summarized as follows:

• From multiple scales analysis, the two types of KdVB equations describing plane progressive pressure waves with low frequency in water flows that contain uniformly distributed translational bubbles were derived. The effect of drag force acting on the bubbles appeared for both the KdVB equation with a correction linear term in (55) and the KdVB equation with a correction nonlinear term in (64).

• The translation of bubbles enhanced the wave nonlinearity and increased the absolute value of the nonlinear coefficient.

• The drag force produced new terms in the KdVB equations and contributed to the dissipation of pressure waves.²⁶ This term enhanced the dissipation for large bubbles (e.g., $R0*≳500 μm$⁠) but contributed to both the nonlinearity and dissipation for small bubbles (e.g., $R0*≲10 μm$⁠).

• The temporal evolution of the dissipation effects due to the acoustic radiation and drag force exhibited distinct trends; the acoustic radiation caused a decrease with an oscillation, whereas the drag force resulted in a nearly monotonic decrease.

• The dissipation effect due to the drag force strongly depended on both the initial void fraction and the initial bubble radius.

This work was partially carried out with the aid of JSPS KAKENHI (Grant No. 18K03942) and the Casio Science Promotion Foundation. We would like to thank the referees for their valuable comments and Editage (www.editage.com) for English language editing.

The data that support the findings of this study are available from the corresponding author upon reasonable request.