Evolution of spherical cavitation bubbles: Parametric and closed-form solutions

It is well established that the size evolution of unstable, spherical cavitation bubbles is governed by the Rayleigh-Plesset (RP) equation^1–3 

In (1)ρ_w is the density of the water, R(t) is the radius of the bubble, p and P_∞ are, respectively, the pressures inside the bubble and at large distance, σ is the surface tension of the bubble, and μ_w is the dynamic viscosity of water. In the simpler form with only the pressure difference in the right hand side, Equation (1) was first derived by Rayleigh in 1917,¹ but it was only in 1949 that Plesset developed the full form of the equation and applied it to the problem of traveling cavitation bubbles.³ In the second half of the last century, a steady progress has been achieved with driving forces from engineering, medical, sonoluminescence, microfluidics, and even pharmaceutical applications of the cavitation phenomena.⁴ In addition, the interest in the analytical and numerical solutions of cavitation dynamics remained considerable as these solutions can lead to a better control and understanding of the bubble collapse processes. In an effort to discern the peculiar features of the usual three-dimensional collapse, Prosperetti² and more recently Klotz⁵ worked out generalizations to N-dimensional bubble dynamics.

In this paper, we first review the analytic solutions of the RP equation in terms of hypergeometric functions when the surface tension is neglected, and in terms of Weierstrass elliptic functions when the surface tension is taken into account.^6,7 In the latter case, we employ an Abel equation approach which is a novel mathematical method of looking to the nonlinear evolution of cavitation bubbles. On the other hand, when the viscous term is introduced, we show that an Abel equation with a non constant invariant occurs. Since it is not yet known how to find analytical solutions to this equation (if any), we resort to numerical integration of the RP equation from t = 0 to t = t_c, where t_c is the time of collapse of the bubble.

We first consider the idealized case whereby the viscosity of the water is neglected, since μ_w ≪ 1, and we further discard the effect of surface tension σ. For this case, there are also recent analytical approximations in Ref. 8, which are further discussed in Ref. 9, but here we are concerned with the standard results.

Let us consider a vacuous p = 0 bubble of radius R which is surrounded by an infinite uniform incompressible fluid, such as water, that is at rest at infinity. We remark that “infinity” in the present context refers to distance far enough away from the initial position of the bubble, and we further assume that the pressure at infinity is constant, P_∞ = const. Neglecting the body forces acting on the bubble, we have from Equation (1)

Since $R 2 R ̇$ is an integrating factor of (1) consequently we obtain by one quadrature

Using the initial conditions R(0) = R₀ and $R ̇ ( 0 ) =0$ we find the integration constant to be

and hence, we obtain

Note that one can find a simple novel particular solution for R(t) by substituting (3) into (2) to obtain the Emden-Fowler equation

with $A=− 3 C 2$⁠, n = 0, m = − 4, and particular solution

However, this solution is obtained under the assumption of a nonzero integrating factor⁸ and therefore it fails to satisfy the second initial condition. Other solutions that do not satisfy the initial conditions have been found previously by Amore and Fernández.⁹ As noticed in Ref. 10, Equation (3) can be also viewed as a conservation law for the dynamics of the radius of the bubble, since its kinetic energy can be expressed as

To proceed with the integration of Equation (3) we will use the set of transformations as given by Kudryashov,⁶ namely, R = S^ϵ, dt = R^δdτ, where ϵ, δ are constants that depend on the dimension of the bubble, and S, τ are the new dependent and independent variables, respectively. Applying the transformations upon (3), we obtain the new dynamics in S and τ

To find S, one sets $ε= 1 N$⁠, and δ = N + 1, where N = 3 is the dimension of the bubble, which will in turn reduce (7) to the simpler equation

By integrating the above with $S ( 0 ) = R 0 3$⁠, we obtain the rational solution

where for convenience we set $B= 9 C 4 R 0 3 = 3 2 P ∞ ρ w R 0 6$⁠. Once we determine S, we can find the parametric solutions for the bubble radius R(τ) and evolution time of the bubble t(τ),^6,

To achieve the time of collapse one needs to allow τ → ∞. This leads to $lim τ → ∞ t ( τ ) =0.000 908 681$ if we use the S.I. units, [ρ_w] = 1000 kg/m³, [R₀] = 10⁻² m, and [P_∞] = 101 325 Pa. We point out that the same asymptotic value for the collapse time has been also obtained, albeit using a different approach, by Obreschkow et al.⁸ and it is also obtained from Eq. (14) below by direct integration. Once we solve for τ as a function of R from the first equation of (10), and substituting it into the second equation of (10) we obtain the closed-form solution

which is plotted as R(t) in Fig. 1.

Next, we will find the time for the total collapse t_c of the bubble by integration of Equation (3), which yields

This is just the N = 3 case of the general N-dimensional formula given in Ref. 5. If we let R = R₀ sin^2/3θ, where θ ∈ [0, π/2] then integral (13) transforms to

which by comparing with the integral relation for beta function, namely,

leads to

This solution can be also obtained from (4) by using the parametric solution given in Section 2.3.1-2 Equation 2 of Polyanin’s book.¹¹ If we insert in (14) the same S.I. numbers as used in the hypergeometric asymptotics, we obtain again t_c = 0.000 908 681 s.

When the surface tension term added, Equation (1) can be written in the form

where we defined $K 1 = P ∞ − P ρ w$ and $K 2 = 2 σ ρ w$⁠.

Proceeding as in Ref. 12, we first show that solutions to a general second order ODE of type

may be obtained via the solutions to Abel’s equation (17) of the first kind (and vice-versa)

using the substitution

which turns (16) into the Abel equation of the second kind in canonical form

Moreover, via the inverse transformation

of the dependent variable, Equation (19) becomes (17). The invariant of Abel’s equation, see Ref. 13, can be written as

and when is a constant is an indication that Abel’s equation is integrable.

By identification of Equation (15) with (16), we see that $f 1 ( R ) = 3 2 R$⁠, f₂(R) = f₀(R) = 0, and $f 3 ( R ) = K 1 R + K 2 R 2$⁠, and hence the Kamke invariant is Φ(R) = 0, therefore Abel’s equation (17) becomes the Bernoulli equation

which, by one quadrature has the solution

By using Equations (20) and (18), we obtain

and via the same initial conditions we obtain the integration constant

which gives

Notice that when σ = 0 → K₂ = 0, and $p=0→ K 1 = P ∞ ρ w$⁠, then the above becomes Equation (3). The new energy with surface tension is

Thus, at the time of collapse we have $4 π 3 R 0 2 ( P ∞ R 0 + 3 σ ) =0.424 43J$ for a surface tension [σ] = 10⁻³ N/m.

To integrate Equation (25), first let us write it in a more convenient way, as

with coefficients defined as $a 3 = R 0 2 K 2 + 2 K 1 3 R 0$⁠, a₁ = − K₂, and $a 0 =− 2 K 1 3$⁠, and we will use the same set of transformations, namely, R = S^ϵ, dt = R^δdτ which give in turn

where now we set ϵ = − 3/N = − 1, and $δ= N + 1 2 =2$⁠. Thus, we obtain the Weierstrass elliptic equation

which in standard form is

via the linear substitution¹⁴ 

The germs of the Weierstrass function g₂, g₃ are given by

where g₂ = 3.377 51 ⋅ 10⁻¹¹ m⁸/s⁴, g₃ = 1.926 45 ⋅ 10⁻⁸ m¹²/s⁶, and the constant in the front of Weierstrass function from (33) takes the value of 59 215.2 s²/m⁵. Once S is known we can find the parametric solutions for the radius and time of the bubble with surface tension as

Related plots are presented in Fig. 2, while in Fig. 3 we compare the exact solutions of RP equation with and without surface tension.

When we include the viscosity, Equation (1) reads

where $f 2 ( R ) = K 3 R 2$⁠, $K 3 = 4 μ w ρ w$ is a constant, and f₁(R), f₃(R) being the same functions as before. Thus, Abel’s equation (17) becomes

Now, we also find the Kamke invariant according to Equation (21) and we obtain

In terms of the physical variables of the system, the invariant is

and because it is not a constant we will try to reduce (36) using the Appell invariant instead.

First, we will eliminate the linear term via the transformation $y ( R ) = R 3 2 z ( R )$ to obtain the reduced Abel equation

where $h 2 ( R ) = K 3 R$⁠, and h₃(R) = (K₁R + K₂) R.

According to the book of Kamke,¹³ for equations of type (39) for which there is no constant invariant one should change the variables according to

which lead to the canonical form

where

is the Appell invariant, and the constants b_i are

In this work, we have considered the RP equation for the size evolution of a bubble in water. In the first part of the paper, we have surveyed the standard results in the absence of surface tension but in a different way from those usually pursued in the literature. We have obtained the closed form hypergeometric solutions of Kudryashov and Sinelshchikov although in a different but equivalent form. From an Emden-Fowler form of the RP equation, we have also obtained a particular solution which, however, does not satisfy the second initial condition. In the presence of surface tension and viscosity, we have employed a new approach based on Abel’s equation. When only the surface tension is included, we have obtained the known parametric rational Weierstrass solutions, whereas when viscosity is added, the corresponding Abel equation does not have a constant invariant, which explains the nonintegrability in this case. A numerical integration obtained from this nonintegrable Abel route is presented graphically.

This research was partially supported by internal funding from Embry-Riddle Aeronautical University. We thank the reviewers for their appropriate remarks which improved the quality of the paper.