Chemical Oscillations of Air-Seeded Bubbles in Water Driven by Ultrasound

Joachim Holzfuss
Institut für Angewandte Physik, TU Darmstadt, Schloßgartenstr. 7, 64289 Darmstadt, Germany

Phys. Rev. E 78, 025303(R) (2008)


Chemical oscillations are shown to be responsible for very low frequency modulations of a bubble oscillating nonlinearly in a high intensity ultrasound field. In the parameter space of incomplete dissociation near the onset of sonoluminescence a small bubble is shown to grow on a long time scale by the intake of dissolved air. Bubble collapses get hotter and more dense, noninert gases are dissociated and removed and a small growing argon bubble is left behind continuing the circle.

A single bubble levitated in water is driven ultrasonically to exhibit extremely nonlinear oscillations. An enormous energy concentration leading to ps-light emission [1] [2] is observed, termed single bubble sonoluminescence (SBSL). Stable SBSL of air-seeded bubbles in water has been linked to a requirement of the presence of a noble gas (argon) [3]. At higher driving pressures, when temperatures exceeding 104 deg Kelvin and pressures of 105 bars exist in the bubble only the chemically inert argon, water vapor and some reaction products remain inside [45]. At slightly lower pressures, an equilibrium exists between rectified diffusion of air into the bubble and partial dissociation of noninert air constituents (curve B in [4]). In this region, interesting dynamical behavior has been found after catastrophic events [6].

Parameter regions have been found experimentally [78], where the clock-like collapse of an air seeded bubble in water looses its regularity. The variability of the collapse time, defined as the interval between some fixed phase of the driving and the moment of light emission has been determined to be mostly of nano second magnitude [72]. But sometimes the distribution spans a couple of micro seconds (Fig. 7 in [8]). The collapse times may not be distributed randomly as there is a frequency associated with this process much lower and quasiperiodic with respect to the driving. The reason of this is unknown.

To solve this phenomenon, a numerical model is developed that includes as much details as is needed to include most relevant physical aspects while maintaining the possibility to calculate several ten-thousand cycles (seconds) of oscillations. As the measured frequency of the modulations is 5 orders of magnitude lower than the main bubble oscillation frequency, computational time restrictions imply simulation of the most important features only. The radial dynamics of a bubble in a compressible liquid is calculated with the Gilmore model [6].

  (     ˙)     (     ˙ )   (    ˙ )    (     ˙)
RR¨ 1- R- + 3R˙2 1 - R-- =H  1 + R-+ R-H˙ 1- R-
       C    2       3C          C   C       C

    ∫ p(R,˙R) -1     p + B   {  ρ}n
H =        ρ  dp , p0 +-B-=  ρ0

            ∘ --||       (           ) (n- 1)∕2n
C = c|   =    dp||   = c0  p(R,R˙)+-B-
      r=R      dρ|r=R         p0 +B

p(R,R˙) =   pg(R )- 2σ-- 4ηR˙        (4)
            (      R) ( R3   3)γ
  pg(R ) =    p0 +-2σ   R-0 --b      (5)
                 R0    R3 - b3
R is the bubble radius, C, ρ, and p(R,) are the speed of sound in the liquid, its density, and the pressure at the bubble wall, respectively. The Tait equation is taken as the equation of state for water using n=7.025, B=3046 bars [9] as parameters. pg is the pressure in the bubble. H is the enthalpy difference of the liquid at pressure p and p(R,) at the bubble wall. p is the pressure at infinity taken as p = p0 + pA cos(2πft), p0 is the ambient pressure increased by hydrostatic pressure above the bubble. The driving frequency f is 14.62 kHz. σ is the surface tension and η is the viscosity of the liquid. The fluid parameters are from tabulated values [10]. b is a van der Waals hard-core radius and γ a polytropic exponent. Its value is set between 1 (=isothermal) and the adiabatic exponent of the gas according to an instantaneous Peclet number [411] Pe = R02|(t)|R(t)κ, reflecting thermal conduction at the involved time scales. κ is the thermal diffusivity of the gas. To avoid a change to isothermal behavior at the bubble wall turning point during maximum compression, is replaced by its maximum during positive accelerations around collapse. It is smoothly approaching the real velocity during the rest of the cycle. The polytropic exponent γ and the hard-core b are taken as an average over the values reflecting the instantaneous content of species, and are thus updated continuously. κ depends on the varying density ρg of the gas κ = ρgρ(g(RR0)) 2 10-5 m2s-1 and is scaled with the ratio of ambient gas density to actual density.

The temperature TB is taken to be uniform within the bubble. It is calculated via compression of a van der Waals gas by TB = T0[(R30 - b3)∕(R3 - b3)]γ-1 with the ambient liquid temperature T0 and the effective polytropic exponent γ.

The number of moles of gases in the bubble ni, i = N2,O2,Ar is changed by diffusion of air through the bubble wall. Because of the slow diffusional time scale an adiabatic approximation [131214] can be employed and the change per period T is

                     (     ⟨        ⟩ )
Δndiiff   4πDiC0iRmax-  ∞     ni
  T    =    Mip0      pi -   n0pg(R ) 4

n0 is the sum of moles of all molecules in the bubble, Mi is the molar mass, Di and Ci are the diffusion constants and concentration fields of the gas species. The concentration of species at the bubble wall is assumed to connect to the partial pressures pi inside according to Henry’s law: Ci|r=R = Ci0pi(R)∕p0. The same law holds true for r = . ⟨f(t)⟩i = 0Tf(t)Ri(t)dt(∫T  i   )
  0 R (t)dt-1 are weighted time averages.

Evaporation and condensation of water molecules at the bubble wall [15316] are included in the model for the bubble dynamics, as experimental results [3] stress the importance of a decrease of the polytropic exponent induced by water vapor at bubble collapse. A Hertz/Knudsen model for the change of moles nH2O of water vapor in the bubble is H2O = H2Oevap -H2Ocond and

˙n      =  4πR2--α¯c(Ts)-(ρsat   - Γ ρ   ) (7)
 H2O      MH2O    4     g,H2O     g,H2O
with the constant evaporation coefficient α (also called accommodation coefficient or sticking probability), and the average velocity of molecules of a Maxwell-Boltzmann distribution ¯c (Ts) = [(8RgasTs)∕(πMH2O )]12. ρg,H2O is the density of water vapor of molar weight MH2O in the bubble, ρg,H2Osat the saturated vapor density [10] and Rgas = 8.3143 J mol-1K-1 is the gas constant. The bubble surface temperature is taken as Ts = T0. The simple model (7) takes the temperature distributions in the bubble and liquid as fixed. Γ is a correction factor for nonequilibrium conditions induced by mass motion of vapor and bubble wall movement [171618]. Calculations show, that Γ varies by as much as 20% around unity during the collapses. However, in the observed parameter range almost no notable difference exists in the amount of water vapor at collapse time when compared to a fixed value. Therefore constant values of Γ = 1, α = 0.4 [6] are taken in the following calculations.


Figure 1: Histogram of collapse times. From data of Fig. 2a

Chemical dissociation occurs for noninert gases [416153196], reaction products are immediately diffused into the liquid. The dissociation per period is calculated as a second order reaction by a modified Arrhenius law:

  diss       ⟨            Ei  ⟩
Δni---= - ni [n0]AiT βie- RgaAsTB-
  T                 B          0

EAi are activation energies, [n0] = n043πR3 is the molar concentration of all molecules, i = N2,O2. Ai and βi are Arrhenius constants [20] . The values of the Arrhenius constant Ai has been shown to depend on the gas mixture (3rd body), especially on the argon content and is changed accordingly.

A drop of the bubble temperature due to water vapor condensation and endothermic chemical reactions is neglected. This contribution is small compared to the heating over a temperature range of several thousand degrees Kelvin. The same holds true for reaction enthalpies, as only small amounts are dissociated during a single collapse. Spatial translations of the bubble in the sound field are implemented [6].

Experiments [8] have determined parameter regions, where a quasiperiodic modulation of bubble characteristics occurs: Low frequency oscillations are seen below the pressure range for stable argon bubbles. Calculations show that in this region incomplete dissociation of noninert gases exists. The rates and the change of rates of diffusion and dissociation may oscillate anti-symmetrically and never reach a static equilibrium. Fig. 1 shows a normalized distribution of calculated collapse times, defined as the interval between the moment the driving pressure crosses zero going to rarefaction and the moment of minimum radius. The broad distribution spanning a couple of micro seconds is substantially the same as in experiments [7], Fig. 7 in [8].


Figure 2: Quasiperiodic modulation on long time scales of an oscillating bubble. From top to bottom as a function of time is shown: a) collapse time, b) maximal radius (straight line) and equilibrium radius (dashed), c) levitation position (height) in the resonator cell.

Fig. 2 shows changes of some characteristics of a bubble on a long time scale. The numerical calculations use parameter settings as in published experimental results [8]. The driving amplitude is 1.194 bars, the ambient pressure 1.033515 bars augmented by a hydrostatic pressure of a water column of 4.36 cm. The water is containing air degassed to 20% ambient concentration. The wave number of the standing wave mode in the resonator equals the published value k = --2π--
1.7179 λ, λ = c0∕f.


Figure 3: Radius-time curves of bubble at minimum and maximum of the oscillation shown in Fig. 2. The small bubble contains 94% argon, the big bubble 79% nitrogen.

In consistency with the experiments, a slow modulation with a period of 3 s is visible in the collapse time and in the maximum radius taken during a single driving period. The collapse time varies by 3μs, the maximum radius has a variability of about 20 μm. The levitation position [6] in the resonator varies by approximately 15 μm in Fig. 2c. This accounts for a very small change in driving pressure of 0.2 Pa during the modulation period. This change cannot be made responsible for the modulation, in agreement with arguments in [8]. The radius-time graphs in Fig. 3 show that the bubble oscillates inertially well above the blake threshold.


Figure 4: Quasiperiodic modulation on long time scales of an oscillating bubble. From top to bottom as a function of time is shown: a) molecules of H2O (straight) and O2 (dotted), b) molecules of N2 (straight) and Ar (dashed), c) dissociation (straight) and diffusion rate (dashed) of N2 and O2 (dotted), d) diffusion rate of Ar, e) temperature (straight) and density (dashed) in the bubble at collapse.

To get further insight into the reasons for these modulations, the gas content of the bubble is analyzed (Fig. 4). Here, modulations of different species and long term changes of the maximum temperature are observed: The variation of gas content species in the bubble is shown in Fig. 4a and b. Non-noble gas contents vary synchronous with the maximum radius. However, the number of molecules of argon changes in an anti correlated manner during parts of the modulation oscillation (Fig. 4b, dashed line). The same is true for the diffusion and dissociation rates of nitrogen in Fig. 4c. Their values cross several times and the changes of rates have different signs: while the diffusion rate decreases, the dissociation rate increases by a large factor. The oxygen diffusion and dissociation rates are equal (dotted line in Fig. 4c) and do not dynamically influence the oscillations. Fig. 4d shows the diffusion rate of argon as it oscillates around zero. The temperature and the density in the bubble (Fig. 4e) show vast changes during the modulation oscillations.

When analyzing Fig. 4 in detail the following explanation of the long term modulations is plausible: When the bubble is small (at t 0 s) (see also Fig . 3) nitrogen is diffusing into and argon out of the bubble because only the nitrogen concentration is large enough to establish a net influx. As the temperature decreases due to the smaller average polytropic exponent within the bubble, the dissociation rate of nitrogen does not keep up with the diffusion rate. Due to the increasingly larger ambient and maximal bubble radius (t 1 s) the argon diffusion rate increases to a net influx. This results in a temperature and nitrogen dissociation rate increase. Shortly before t 3 s the nitrogen dissociation rate overtakes its diffusion rate while the argon rate still increases as the bubble oscillations are fairly large now. At t 3.2 s the temperature increases to 9000 K which, together with a high density of 0.7 g/cm3 results in a sudden nitrogen loss of 2 orders of magnitude. The bubble now contains mostly argon. Its oscillations get small, the temperature, the density and the nitrogen dissociation rate drop sharply. The cycle repeats itself with an increasing influx of nitrogen. The modulation oscillation occurs as a limit cycle in a bistable system, where both equilibrium states are unstable.


Figure 5: Slow quasiperiodic oscillations of the maximum radius as a function of time for decreasing driving amplitude (from top to bottom: 119800 Pa, 119600 Pa, 119500 Pa, 119400 Pa and 119200Pa (dotted)) The driving frequency is 14.620 kHz, ambient pressure 103351 Pa augmented by the hydrostatic pressure at water depth 4.36 cm.

The dependence of the slow oscillations on the driving pressure is seen in the time dependence of the modulation of the maximum radius of the bubble taken during a single driving period (Fig. 5). With decreasing driving pressure the frequency of the modulation decreases and its amplitude increases. The same behavior has been reported in the experiments in [8]. In Fig. 5a the bubble settles to a non modulated oscillating volume containing an N2∕Ar mixture (ratio 6040) with some O2 and H2O. With decreasing driving amplitude modulation oscillations set in. The dotted line in Fig. 5d shows the lower pressure limit of the modulations. The long term transient of the maximum radius shows that the bubble is getting so large that it enters the region of parametric instability. Here, micro bubble pinchoff can occur [6]. Calculations verify, that the reported modulations are robust with respect to changes in the amount of allowed water vapor and to a more complete model of thermal damping [3]. The observed effects may well play a role during path instabilities of bubbles in sulfuric acid [21].

Oscillations occur in a large number of chemical systems in different configurations [22]. The system described here is a forced system which is open by continuous inflow of reactants, non-isothermal and non-isobaric, whereby the Arrhenius rates are changed nonlinearly. It has analogies to the nonlinear dynamics in a continuously stirred tank reactor (CSTR), where Hopf-bifurcations, saddle-node bifurcations and multistabilities have been found [23]. This system is unique in that it shows a slow physico-chemical oscillation initiated by nonlinear dynamics on a fast time scale of a bubble.

The author thanks C. R. Thomas and R. G. Holt for making available their data prior to publication.


[1]   D. F. Gaitan, L. A. Crum, C. C. Church, and R. A. Roy, J. Acoust. Soc. Am. 91, 3166 (1992).

[2]   B. P. Barber et al., Phys. Rep. 281, 65 (1997), B. Gompf et al., Phys. Rev. Lett. 79, 1405 (1997).

[3]   M. P. Brenner, S. Hilgenfeldt, and D. Lohse, Rev. Mod. Phys. 74, 425 (2002).

[4]   D. Lohse and S. Hilgenfeldt, J. Chem. Phys. 107, 6986 (1997).

[5]   Y. T. Didenko and K. S. Suslick, Nature 418, 394 (2002).

[6]   J. Holzfuss, Phys. Rev. E 71, 026304 (2005).

[7]   R. G. Holt, D. F. Gaitan, A. A. Atchley, and J. Holzfuss, Phys. Rev. Lett. 72, 1376 (1994).

[8]   C. R. Thomas, R. A. Roy, and R. G. Holt, Phys. Rev. E 70, 066301 (2004).

[9]   J. Holzfuss, M. Rüggeberg, and A. Billo, Phys. Rev. Lett. 81, 5434 (1998).

[10]   At 20 deg C the following parameters have been taken in the calculations: surface tension σ = 0.0725 N m-1, viscosity η = 0.001 N s m-2, saturated vapor density ρg,H2Osat = 0.0173 kg m-3, diffusion constant DAr = 2.0 × 10-9 m2s-1, DO2 = 1.9 × 10-9 m2s-1, DN2 = 1.9 × 10-9 m2s-1, solubility Ar = 0.061 kg m-3, O2 = 0.0444 kg m-3, N2 = 0.020 kg m-3, water density ρ = 998.23 kg m-3, sound velocity c0 = 1483 ms-1.

[11]   A. Prosperetti, J. Acoust. Soc. Am. 61, 17 (1977).

[12]   D. Lohse et al., Phys. Rev. Lett. 78, 1359 (1997).

[13]   R. Löfstedt, B. P. Barber, and S. J. Putterman, Phys. Fluids A 5, 2911 (1993), M. M. Fyrillas and A. J. Szeri, J. Fluid Mech. 277, 381 (1994).

[14]   R. Löfstedt, K. R. Weninger, S. Putterman, and B. P. Barber, Phys. Rev. E 51, 4400 (1995), S. Hilgenfeldt, D. Lohse, and M. P. Brenner, Phys. Fluids 8, 2808 (1996).

[15]   V. Kamath, A. Prosperetti, and F. N. Egolfopoulos, J. Acoust. Soc. Am. 94, 248 (1993).

[16]   B. D. Storey and A. J. Szeri, Proc. R. Soc. London A 456, 1685 (2000), K.Yasui, T. Tuziuti, T. Kozuka, A. Towata, Y. Iida, J. Chem. Phys. 127, 154502 (2007).

[17]   R. W. Schrage, A theoretical study of interphase mass transfer (Columbia University Press, New York, 1953).

[18]   J. Holzfuss and M. T. Levinsen, Phys. Rev. E 77, 046304 (2008).

[19]   R. Toegel and D. Lohse, J. Chem. Phys. 118, 1863 (2003), L. S. Bernstein, M. R. Zakin, E. B. Flint, and K. S. Suslick, J. Phys. Chem. 100, 6612 (1996).

[20]   Arrhenius constants: (Ai[m3mole-1s-1]i,EAi∕Rgas[K]) for i = N2 in N2 (2.3 1023,-3.5,113200), N2 in Ar (0.4 2.3 1023,-3.5,113200), O2 in N2 (3.4 1012,-1,59380), O2 in Ar (1.6 1012,-1,59380) [24].

[21]   D. J. Flannigan and K. S. Suslick, Phys Rev Lett 95, 044301 (2005), R. Toegel, S. Luther, and D. Lohse, Phys. Rev. Lett. 96, 114301 (2006).

[22]    P. Gray and S. K. Scott, Chemical Oscillations and Instabilities. Non-linear Chemical Kinetics., Vol. 21 of International Series of Monographs on Chemistry (Clarendon Press, Oxford, 1990).

[23]   J. Guckenheimer, Physica D 20, 1 (1986), I. G. Kevrekidis, R. Aris, and L. D. Schmidt, Physica D 23, 391 (1986).

[24]   D. J. Kewley and H. G. Hornung, Chem. Phys. Lett. 25, 531 (1974), D. L. Baulch, D. D. Drysdale, and D. G. Horne, Evaluated kinetic data for high temperature reactions (Butterworths, London, 1972), Vol. 2, L. Jerig, K. Thielen, and P. Roth, AIAA Jornal 29, 1136 (1991).