Numerical modeling of the dynamics of bubble oscillations subjected to fast variations in the ambient pressure with a coupled level set and volume of fluid method

A numerical method for modeling and understanding the dynamics of bubble oscillations subjected to fast variations in the ambient pressure is proposed under low Mach number conditions. In the present work, the method uses a single-fluid continuum formalism of weakly compressible axisymmetric Navier-Stokes equations for the numerical simulation of liquid-gas flows with surface tension and adopts the interface capturing approach based on a coupled level set and volume of fluid (CLSVOF) method for describing the moving and deformed interfaces. To demonstrate the efficacy of the proposed method, first, the numerical results of the radial oscillations of a spherical gas bubble are tested with the numerical solutions of Rayleigh-Plesset equation. Then, the numerical method is applied to reproduce the growth and subsequent collapse of a bubble in an infinite liquid medium observed in experiments. Finally, the numerical simulation of the interaction of two oscillating bubbles at small separation distance is evaluated in response to a moderate step change in the ambient pressure. It is shown that two deformable bubbles undergo coupled radial and oscillatory translational motions which eventually results in the bubbles' attraction and coalescence caused by the secondary Bjerknes forces. The numerical predictions show very good accuracy with the experimental and numerical results reported in the literature.

Figure 1

(a) Schematic diagram of the problem geometry and axisymmetric computational domain on a $(r,z)$ plane used to model the behavior of a single gas bubble oscillations due to changes in the ambient pressure $p_{\infty }$. (b) Computational two-phase cell with piecewise-linear interface segment and distribution of volume of liquid fraction $\alpha$ around the interface. The width and height of the computational domain are denoted as $R_{\infty }$ and $L_{\infty }$ (see in Sec. 2d), respectively.

Figure 2

Comparison of numerically computed radial oscillations of a spherical gas bubble with the numerical solution of RP equation [1] using the parameters (a) $\varepsilon ={10}^{-1}$ and $R_0=1\mathrm{mm}$, (b) $\varepsilon ={10}^0$ and $R_0=50\mu \mathrm{m}$, (c) $\varepsilon ={10}^{-1}$ and $R_0=1\mathrm{mm}$, and $f_f=5$ kHz, and (d) $\varepsilon ={10}^0$ and $R_0=50\mu \mathrm{m}$ and $f_f=80$ kHz. The other input parameters taken as follows: $p_{\text{atm}}={10}^5$ Pa, ${\rho }_l=1000$ kg ${\mathrm{m}}^{-3},{\rho }_{g,\text{ref}}=1$ kg ${\mathrm{m}}^{-3},{\mu }_l={10}^{-3}$ kg ${\mathrm{m}}^{-1}{\mathrm{s}}^{-1},{\mu }_g={10}^{-5}$ kg ${\mathrm{m}}^{-1}{\mathrm{s}}^{-1},\sigma =0.072$ kg ${\mathrm{s}}^{-2}$, and $\gamma =1.4$, with $800\times 800$ grid cells.

Figure 3

Effect of computational grid cells on the temporal variation of the spherical bubble radius using the same parameters as in Fig. 2.

Figure 4

Effect of grids on the radius profile comparison with RP solution during the first oscillation with the same parameters as in Fig. 2. Arrow indicates the zoom view of the variations of the minimum radius for different grids. The minimum radius of the bubble is $R_{\min ,\mathrm{RP}}=9.554\times {10}^{-1}\mathrm{mm}$ at $t=1.42\times {10}^{-4}\mathrm{s}$.

Figure 5

The time evolution of relative gas-mass error under the same conditions of Fig. 2. Black dash-dotted line: $\text{CFL}=1.0$; blue dash line: $\text{CFL}=0.5$; red solid line: $\text{CFL}=0.25$.

Figure 6

(a) Pressure distribution with velocity field during first contraction phase. For the first cycle, the time period for contraction phase is 0 s $\le t\le 1.42\times {10}^{-4}$ s; (b) pressure distribution with velocity field during first expansion phase of a spherical gas bubble. For the first cycle, the time period for expansion phase is $1.42\times {10}^{-4}$ s $<t\le 2.84\times {10}^{-4}$ s. The parameters used here are the same as in Fig. 2.

Figure 7

(a) Imposed oscillatory pressure field at the far boundary and (b) Comparison of the history of the equivalent radius of the bubble between experimental measurements [9], Rayleigh-Plesset equation for oscillating spherical bubble system, and present numerical results. Here, $\mathrm{A}$ and $\mathrm{B}$ inside the figure refer to the end of first expansion phase with maximum bubble radius and the end of first collapse phase with minimum bubble radius, respectively.

Figure 8

Snapshots of the time evolution of the bubble shapes during (a) first expansion phase at times $t$ ($t^{*}$) (1) 0.0000 ms (0), (2) 0.1165 ms (6.47), (3) 0.1731 ms (9.62); (b) first collapse phase at times $t$ ($t^{*}$) (1) 0.1731 ms (9.62), (2) 0.2004 ms (11.13), (3) 0.2066 ms (11.48), (4) 0.2105 ms (11.69), (5) 0.2115 ms (11.75); (c) magnified snapshots of the evolution of the bubble shapes at the final stages of the collapse phase at times $t$ ($t^{*}$) (1) 0.2097 ms (11.65), (2) 0.2105 ms (11.69), (3) 0.2111 ms (11.73), (4) 0.2113 ms (11.74), (5) 0.2115 ms (11.75). The bubble reaches the first maximum radius at the end of the expansion phase and the first minimum radius at the end of the collapse phase at $t=0.1731\mathrm{ms}$ or $t^{*}=9.62$ and $t=0.2115\mathrm{ms}$ or $t^{*}=11.75$, respectively. Note that the horizontal black marks of the left side near vertical axis denote the initial position of the bubble. These black marks indicate the displacement of the bubble center from its initial position. The input parameters are same as used in Fig. 7.

Figure 9

Snapshots of the evolution of the bubble shapes during the rebound phase (second expansion phase). Here, the dimensionless time $t^{*}$ can be written as (a) 0.11.75, (b) 11.77, (c) 11.78, (d) 11.84, (e) 12.12, and (f) 12.22.

Figure 10

Bubble interface with pressure contours and velocity vectors during (a)–(c) first expansion phase, (d)–(f) first collapse phase, and (g)–(i) early stages of rebound phase under the same conditions as in Fig. 7.

Figure 11

The time evolution of (a) the position of of the extremes of the bubble located on the symmetric axis of the computational geometry, (b) relative velocity of the top and bottom extremes of the bubble, (c) the center position of the bubble, and (d) the bubble aspect ratio, $\beta =\frac{d_z}{d_r}$. The input parameters are the same as in Fig. 7.

Figure 12

Schematic view of the problem geometry used to consider the interaction of two coaxial bubbles in response to fast variations in the ambient pressure $p_{\infty }$ in an axisymmetric computational domain on a $(r,z)$ with independent $\theta$.

Figure 13

(a) Comparison of the present computed history of the equivalent radius of the lowermost bubble $R_2^{*}$ with the numerical result of Chatzidai et al. [23]. The time evolution of (b) $Z_{\text{center},1}^{*},Z_{\text{center},2}^{*}$ and $Z_c^{*}$, (c) $D^{*}$, and (d) $U_1^{*}$ and $U_2^{*}$. The pertinent input parameters are $p_{\text{atm}}^{*}=100,\varepsilon =1,\gamma =1.4,R_{10}^{*}=R_{20}^{*}=1.0,D^{*}\left(0\right)=2.8,\text{Oh}=0.1,\lambda =3.652\times {10}^{-4}$, and ${\eta }_{\text{init}}={10}^{-3}$. Here, in (b) and (d) the solid line is for the lowermost bubble and the dashed line for the uppermost bubble.

Figure 14

Bubble shapes with pressure contours and velocity vectors during the interactions of two equal sized bubbles with the same parameters as in Fig. 13. Here, (a) and (d) expansion phase, and (b) and (c) contraction phase.

Figure 15

Time instant of the coalescence of the two coaxial bubbles and an oscillating compound bubble formation with the same parameters as in Fig. 13.

Figure 16

(a) The computed history of the equivalent radius of two oscillating bubbles. The time evolution of (b) $Z_{\mathrm{center},1}^{*},Z_{\mathrm{center},2}^{*}$ and $Z_c^{*}$, (c) $D^{*}$, (d) $U_1^{*}$ and $U_2^{*}$, and (e) bubble shapes at times $t^{*}$ (1) 0.02, (2) 0.44, and (3) 0.68, under a step change in surrounding pressure $(f_f^{*}=0$). The input parameters are $p_{\mathrm{atm}}^{*}=100,\varepsilon =1,\gamma =1.4,R_{10}^{*}=1,R_{20}^{*}=0.7,D^{*}\left(0\right)=2.8,\mathrm{Oh}=0.1,\lambda =3.652\times {10}^{-4}$ and ${\eta }_{\mathrm{init}}={10}^{-3}$.

Figure 17

(a) The computed history of the equivalent radius of two oscillating bubbles. The time evolution of (b) $Z_{\text{center},1}^{*},Z_{\text{center},2}^{*}$, (c) $D^{*}$, (d) $U_1^{*}$ and $U_2^{*}$, and (e) bubble shapes at times $t^{*}$ (1) 0.22, (2) 0.85, (3) 1.64, and (4) 2.12 under an oscillatory pressure field $(f_f^{*}=4.07$). The pertinent input parameters are $p_{\text{atm}}^{*}=100,\varepsilon =0.3,\gamma =1.4,R_{10}^{*}=1,R_{20}^{*}=0.7,D^{*}\left(0\right)=4.0,\text{Oh}=0.1,\lambda =3.652\times {10}^{-4}$, and ${\eta }_{\text{init}}={10}^{-3}$.

Figure 18

(a) The computed history of the equivalent radius of two oscillating bubbles. The time evolution of (b) $Z_{\text{center},1}^{*},Z_{\text{center},2}^{*}$ and $Z_c^{*}$, (c) $D^{*}$, (d) $U_1^{*}$ and $U_2^{*}$, and (e) bubble shapes at times $t^{*}$ (1) 0.27, (2) 0.82, and (3) 2.02 under an oscillatory pressure field $(f_f^{*}=4.07$). The input parameters are $p_{\text{atm}}^{*}=100,\varepsilon =0.3,\gamma =1.4,R_{10}^{*}=1,R_{20}^{*}=0.7,D^{*}\left(0\right)=2.8\text{Oh}=0.1,\lambda =3.652\times {10}^{-4}$, and ${\eta }_{\text{init}}={10}^{-3}$.