© The Author(s), published by EDP Sciences, 2022

Nomenclature

θ _α : Azimuth of bubble α

θ _β : Azimuth of bubble β

S _α : Radius of bubble α

S _β : Radius of bubble β

P₂(cosθ): Second-order Legendre function

a₀(t): Spherical radius component of bubble α

b₀(t): Spherical radius component of bubble β

a₂(t): Coefficient of nonspherical radius component of bubble α

b₂(t): Coefficient of nonspherical radius component of bubble β

ϵ : Small parameter

r _α : Distance from a point in the liquid to the center of bubble α

r _β : Distance from a point in the liquid to the center of bubble β

L : Distance between bubble centers

f _s : Electric stress

ε ₀ : Vacuum permittivity

ε _rl : Relative permittivity

E : Electric field strength in the liquid region

p ₀ : Static pressure of the liquid

p _v : Vapor pressure in the bubble

p _a : Acoustic pressure amplitude

σ : Surface tension at the liquid–gas interface

η : Viscosity of the liquid

γ : Adiabatic index

R ₀ : Initial bubble radius

R : Bubble radius

Ṙ : Interface velocity

ω : Angular frequency of the acoustic wave

1 Introduction

As an efficient, environmentally-friendly and economical technology, the combination of an electric field with ultrasound has been widely used in laboratories and industrial fields in recent years. A study by Jung et al. [1] revealed that this combined treatment can be used to disintegrate waste activated sludge. After 60 min treatment time, the average particle size of the sludge decreased significantly more, compared to ultrasonic treatment alone. In addition, Yang et al. [2] used the ultrasound-synergized electric field extraction (UEE) to extract flavonoids from hemerocallis citrina baroni (a citron daylily). They found that the yield of flavonoids was 10% higher using UEE than conventional ultrasound extraction, and it was 18% higher than water bath extraction. A work by Yang et al. [3] reveals the mechanism how electric stress strengthens ultrasonic cavitation. However, many details of interactions between nonspherical bubbles are not known. In fact, when cavitation occurs in the liquid, it can generate a high number of cavitation bubbles, instead of a single one. The interactions between these adjacent bubbles affect the oscillation of the individual bubbles significantly, which should be considered in the mathematic model. To describe the effect of an adjacent bubble on the cavitation of an investigated bubble, we need to develop a dynamic equation-set including the interaction between double bubbles.

Many researchers helped advance the theoretical model for two interacting bubbles. Li et al. [4], Pityuk et al. [5], and Wang et al. [6] provided three different mathematical models for two adjacent spherical bubbles in an ultrasonic field. In the study of CCl₄ sono-pyrolysis, Merouani et al. [7–9] used modified Keller-Miksis equations to analyze the thermal effect and the interaction between bubbles. Shima and Fujiwara [10] investigated the behavior of two spherical bubbles with different sizes near a solid wall and derived the motion equations for two bubbles. However, these models are unsuitable for two interacting nonspherical bubbles in the presence of an electric force. The variation of nonspherical radius components of bubbles cannot be obtained by simply modifying the models for spherical bubbles. Liang et al. [11] developed a dynamic model for two native nonspherical bubbles in an incompressible liquid. However, their derivation includes some important inaccuracies, which leads to incorrect mathematic equations. Furthermore, in the calculation of the total velocity potential near the interface of a bubble, they did not consider the perturbing velocity potential induced by the interaction between bubbles, which is an additional term of the total velocity potential. In addition, none of the above-mentioned studies considered the liquid’s compressibility. If the liquid’s compressibility is neglected in the numerical calculation, the maximum bubble-radius can easily be larger, and the minimum bubble-radius can be smaller, compared to the actual situation.

The present study aims to develop a dynamic model for nonspherical double-bubbles in a compressible liquid under the coupling effect of ultrasound and an electrostatic field. We used perturbation analysis, the sphere theorem, and the basic fluid-mechanics equations to derive the required model. We mainly analyze the effect of the interaction force on the investigated bubble using numerical solutions to the theoretical model. The variable parameters are: the distance between bubble centers, and the initial radius of the adjacent bubble. Furthermore, in the derivation, we offer an explanation why the effect of the liquid’s compressibility cannot transfer from the area near the interface of one bubble to the one of another bubble, when the two bubbles are far apart from each other. Finally, we applied approximate equations to investigate temporal evolutions of the gas pressure and the temperature inside a bubble. These results can provide better insights into the oscillation mechanism of interacting nonspherical bubbles in a compressible liquid.

2 Model for double nonspherical bubbles

2.1 Interface equations for double bubbles

To describe a double-bubble system, we introduce a nonspherical function S(θ, t) with rotational symmetry, which is independent of the rotation angle ϕ. In spherical coordinates, the interface equations for double bubbles can be expressed as

$$S_{\alpha }=S_{\alpha }\left({\theta }_{\alpha },t\right)=a_0\left(t\right)+ϵa_2\left(t\right)P_2\left(\cos {\theta }_{\alpha }\right),$$(1a)

$$S_{\beta }=S_{\beta }\left({\theta }_{\beta },t\right)=b_0\left(t\right)+ϵb_2\left(t\right)P_2\left(\cos {\theta }_{\beta }\right),$$(1b)

where

$$P_2\left(\cos {\theta }_{\alpha }\right)=\frac{1}{2}\left(3{\cos }^2{\theta }_{\alpha } -1\right),$$(2a)

$$P_2\left(\cos {\theta }_{\beta }\right)=\frac{1}{2}\left(3{\cos }^2{\theta }_{\beta } -1\right).$$(2b)

The small parameter ϵ indicates that, for both bubbles, the nonspherical components are much smaller than the spherical components.

Figure 1 shows the sections of two nonspherical bubbles, which are elliptic. The distance L is set to hundreds of microns in this study, namely L ≫ S_α + S_β. The reason is that if two bubbles are too close, the coalescence of these two bubbles probably occurs. This is out of the scope of our discussion.

Figure 1 Sections of two nonspherical bubbles.

2.2 Total pressure-difference on the whole liquid region outside a bubble

The direction of the electric field is set to parallel to the line of bubble centers. For an insulating bubble, the expression for electric stress f_s [12, 13] is as follows:

$$f_{\mathrm{s}}={\epsilon }_0\frac{{\left({\epsilon }_{\mathrm{r}\mathrm{l}}-1\right)}^2}{3}{E}^2 P_2\left(\cos \theta \right).$$(3)

The electric stress f_s exerts on the bubble wall, stretching the bubble in the electric-field direction and compressing the bubble in the direction vertical to the electric field.

If a bubble undergoes adiabatic oscillation, the total pressure-difference on the whole liquid region outside this bubble can be written as

$$\Delta P=\left(p_0+\frac{2\sigma }{R_0}-p_{\mathrm{v}}-ϵf_{\mathrm{s}}\right){\left(\frac{R_0}{R}\right)}^{3\gamma }+p_{\mathrm{v}}-p_0-\frac{2\sigma }{R}-\frac{4\eta \dot{R}}{R}+ϵf_{\mathrm{s}}-p_a\cos \omega t=\Delta P_1+ϵf_{\mathrm{e}}{\mathrm{P}}_2\left(\cos \theta \right),$$(4)

where

$$\Delta P_1=\left(p_0+\frac{2\sigma }{R_0}-p_{\mathrm{v}}\right){\left(\frac{R_0}{R}\right)}^{3\gamma }+p_{\mathrm{v}}-p_0-\frac{2\sigma }{R}-\frac{4\eta \dot{R}}{R}-p_{\mathrm{a}}\cos \omega t,$$(5)

$$f_{\mathrm{e}}=\left[1-{\left(\frac{R_0}{R}\right)}^{3\gamma }\right]{\epsilon }_0\frac{{\left({\epsilon }_{\mathrm{r}\mathrm{l}}-1\right)}^2}{3}{E}^2.$$(6)

For a nonspherical bubble, R = a₀(t) or b₀(t). To simplify the calculation, the bubble is considered to be approximately spherical for the force analysis (Eqs. (6)–(8)). The rationality of this simplification is based on that the nonspherical term is much smaller than the spherical term for an elliptical bubble.

2.3 Derivation of the theoretical model

Let coordinates α be the spherical coordinates with the origin at the center of bubble α, and coordinates β be the spherical coordinates with the origin at the center of bubble β.

In coordinates α, the velocity potential of bubble α can be written as

$${\mathrm{\Phi }}_{\alpha }\left(r_{\alpha },{\theta }_{\alpha },t\right)={\mathrm{\Phi }}_{\alpha 0}\left(r_{\alpha },t\right)+ϵ{\mathrm{\Phi }}_{\alpha 2}\left(r_{\alpha },{\theta }_{\alpha },t\right),$$(7)

where

$${\mathrm{\Phi }}_{\alpha 0}\left(r_{\alpha },t\right)\approx \frac{A_0\left(t\right)}{k{r}_{\alpha }}-\frac{{\dot{A}}_0\left(t\right)}{\omega },$$(8)

$${\mathrm{\Phi }}_{\alpha 2}\left(r_{\alpha },{\theta }_{\alpha },t\right)\approx P_2\left(\cos {\theta }_{\alpha }\right)\left[\frac{3A_2\left(t\right)}{k^3{r}_{\alpha }^3}+\frac{2A_2\left(t\right)}{k{r}_{\alpha }}+\frac{{\dot{A}}_2\left(t\right)}{\omega }\right].$$(9)

Here, Φ_α0 (r_α, t) and Φ_α2(r_α, θ_α, t) are the spherical component and the nonspherical component of the velocity potential of bubble α. k is the wave number, k = ω/c, and c is the speed of sound in the liquid. A₀(t) and A₂(t) are two functions of time. The derivations of equations (10) and (11) (see Ref. [3]) are based on the assumption that ultrasonic wave propagates within an infinitely large liquid region. Hence, no standing wave occur in the liquid. This means that translational motions of bubbles cannot occur. Thus, the distance L can be considered constant.

In the dynamic model for bubble α, the interaction terms reflect the effects of bubble β on bubble α. The effect of the liquid’s compressibility only occurs in the vicinity of the interface of a bubble, which cannot transfer to the liquid region far away from the interface of a bubble. Importantly, this effect cannot transfer to the vicinity of the interface of another bubble, when the two bubbles are far apart from each other. The reason for this is as follows: In the radial direction, as the distance from a bubble center to a point in the liquid increases, the fluid velocity decreases sharply. The effect of the liquid’s compressibility is positively correlated with the fluid velocity. For a liquid region far away from the bubble interface, the fluid velocity is very small, the effect of the liquid’s compressibility can therefore be neglected. Hence, in the derivation of the dynamic model for bubble α, the liquid’s compressibility needs not to be considered when calculating the velocity potential of bubble β.

In coordinates β, the velocity potential of bubble β is

$${\mathrm{\Phi }}_{\beta }\left(r_{\beta },{\theta }_{\beta },t\right)={\mathrm{\Phi }}_{\beta 0}\left(r_{\beta },t\right)+ϵ{\mathrm{\Phi }}_{\beta 2}\left(r_{\beta },{\theta }_{\beta },t\right),$$(10)

where

$${\mathrm{\Phi }}_{\beta 0}\left(r_{\beta },t\right)=\frac{B_0\left(t\right)}{r_{\beta }},$$(11)

$${\mathrm{\Phi }}_{\beta 2}\left(r_{\beta },{\theta }_{\beta },t\right)=\frac{B_2\left(t\right)}{r_{\beta }^3} P_2\left(\cos {\theta }_{\beta }\right).$$(12)

Here, Φ_β0(r_β, t) and Φ_β2(r_β, θ_β, t) denote the spherical component and the nonspherical component of the velocity potential of bubble β. B₀(t) and B₂(t) are two functions of time. Equations (11) and (12) can be obtained from a study by Chen et al. [14].

To calculate the total velocity potential for a point in the liquid, we need to convert the expression of the velocity potential of bubble β from coordinates β to coordinates α. When r_α is near the interface of bubble α, namely r_α ≪ L, we have the following relations [15–19]:

$$\frac{1}{r_{\beta }}=\frac{1}{\sqrt{L^2-2L{r}_{\alpha }\cos {\theta }_{\alpha }+r_{\alpha }^2}}=\frac{1/L}{\sqrt{1-2\left(r_{\alpha }/L\right)\cos {\theta }_{\alpha }+{\left(r_{\alpha }/L\right)}^2}}=\sum _{i=0}^{\mathrm{\infty }}\mathrm{}\frac{r_{\alpha }^i}{L^{i+1}}{P}_i\left(\cos {\theta }_{\alpha }\right),$$(13)

$$\frac{P_2\left(\cos {\theta }_{\beta }\right)}{r_{\beta }^3}=\frac{1}{L^3}\sum _{i=0}^{\mathrm{\infty }}\mathrm{}{C}_{2+i}^2{\left(\frac{r_{\alpha }}{L}\right)}^i{P}_i\left(\cos {\theta }_{\alpha }\right),$$(14)

where

$$C_{2+i}^2=\frac{(2+i)!}{2! i!}=\frac{(i+2)(i+1)}{2},$$

and P_i(cosθ_α) is the ith-order Legendre function.

Substituting these relations into equations (11) and (12), we obtain the velocity potentials of bubble β in coordinates α:

$${\mathrm{\Phi }}_{\beta 0}\left(r_{\alpha },t\right)=\frac{B_0\left(t\right)}{L}\sum _{i=0}^{\mathrm{\infty }}\mathrm{}{\left(\frac{r_{\alpha }}{L}\right)}^i{P}_i\left(\cos {\theta }_{\alpha }\right),$$(15)

$${\mathrm{\Phi }}_{\beta 2}\left(r_{\alpha },{\theta }_{\alpha },t\right)=\frac{B_2\left(t\right)}{L^3}\sum _{i=0}^{\mathrm{\infty }}\mathrm{}{C}_{2+i}^2{\left(\frac{r_{\alpha }}{L}\right)}^i{P}_i\left(\cos {\theta }_{\alpha }\right).$$(16)

In addition, according to the sphere theorem by Weiss [20], when bubble α exists in an arbitrary irrotational flow, it can cause a disturbance to this flow. Thus, an extra perturbing velocity potential can occur in the vicinity of the interface of bubble α. When the flow is undisturbed by bubble α, the unperturbed velocity potential of the flow should be Φ_β. Hence, under the disturbance of bubble α, we denote the perturbing velocity potential as Φ_β→α(r_α, θ_α, t), which is also the additional velocity potential induced by bubble β in the vicinity of the interface of bubble α.

The perturbing velocity potential can be calculated by the following equality [20]:

$${\mathrm{\Phi }}_{\beta \to \alpha }\left(r_{\alpha },{\theta }_{\alpha },t\right)=\frac{S_{\alpha }}{r_{\alpha }}{\mathrm{\Phi }}_{\beta }\left(\frac{S_{\alpha }^2}{r_{\alpha }},{\theta }_{\alpha },t\right)-\frac{2}{S_{\alpha }{r}_{\alpha }}{\int }_0^{S_{\alpha }}\mathrm{}\lambda {\mathrm{\Phi }}_{\beta }\left(\frac{{\lambda }^2}{r_{\alpha }},{\theta }_{\alpha },t\right)\mathrm{d\lambda }.$$

Substituting equations (15) and (16) into the above equality, we obtain

$${\mathrm{\Phi }}_{\beta \to \alpha }(r_{\alpha },{\theta }_{\alpha },t)=\frac{B_0\left(t\right)}{L}\sum _{i=0}^{\mathrm{\infty }}\mathrm{}\frac{i}{i+1}\frac{S_{\alpha }^{2i+1}}{r_{\alpha }^{i+1}{L}^i}{P}_i\left(\cos {\theta }_{\alpha }\right)+ϵ\frac{B_2\left(t\right)}{L^3}\sum _{j=0}^{\mathrm{\infty }}\mathrm{}{C}_{2+j}^2 \frac{j}{j+1}\frac{S_{\alpha }^{2j+1}}{r_{\alpha }^{j+1}{L}^j}{P}_j\left(\cos {\theta }_{\alpha }\right).$$(17)

Similarly, Φ_α→β(r_α, θ_α, t) represents the perturbing velocity potential induced by bubble α in the vicinity of the interface of bubble β. However, we assume that L ≫ S_α + S_β. Thus, the effect of Φ_α→β(r_α, θ_α, t) on the vicinity of the interface of bubble α is extremely small. By calculation, we can formulate Φ_α→β(r_α, θ_α, t) = O(1/L¹⁰) in coordinates α. Hence, the term Φ_α→β(r_α, θ_α, t) can be neglected for the calculation of the total velocity potential.

In addition, the total velocity potential Φ(r_α, θ_α, t) should be consistent with the radius function of bubble α (see Eqs. (1a)) in form. Thus, in the nonspherical term of Φ(r_α, θ_α, t), only the terms about the second-order Legendre function P₂(cosθ_α) can be retained. The expression for the total velocity potential is:

$$\mathrm{\Phi }(r_{\alpha },{\theta }_{\alpha },t)={\mathrm{\Phi }}_{\alpha }(r_{\alpha },{\theta }_{\alpha },t)+{\mathrm{\Phi }}_{\beta }(r_{\alpha },{\theta }_{\alpha },t)+{\mathrm{\Phi }}_{\beta \to \alpha }(r_{\alpha },{\theta }_{\alpha },t)\approx \frac{A_0\left(t\right)}{k{r}_{\alpha }}-\frac{{\dot{A}}_0\left(t\right)}{\omega }+\frac{B_0\left(t\right)}{L}+ϵP_2\left(\cos {\theta }_{\alpha }\right)\{\frac{3A_2\left(t\right)}{k^3{r}_{\alpha }^3}+\frac{2A_2\left(t\right)}{k{r}_{\alpha }}+\frac{{\dot{A}}_2\left(t\right)}{\omega }+\frac{B_2\left(t\right)}{L^5}\left[6r_{\alpha }^2+\frac{4S_{\alpha }^5}{r_{\alpha }^3}\right]\}.$$(18)

At the interface of bubble α, the total velocity potential satisfies the following two equations:

$${\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }{r}_{\alpha }}|}_{r_{\alpha }=S_{\alpha }}=\frac{\mathrm{d}{S}_{\alpha }}{\mathrm{d}t},$$(19)

$${\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }t}|}_{r_{\alpha }=S_{\alpha }}+{\frac{1}{2}(\nabla \mathrm{\Phi }{)}^2|}_{r_{\alpha }=S_{\alpha }}=-\frac{\Delta P}{\rho },$$(20)

where the expression of ΔP see equation (4). Equation (20) is obtained when the gravity of the liquid is negligible and the liquid flow is irrotational. Substituting equations (1a) and (18) into equations (19) and (20), using the Taylor series expansion of ϵ, and only retaining the first two order terms for ϵ, we obtain

$${\dot{a}}_0\left(t\right)+\frac{1}{k{a}_0^2\left(t\right)}{A}_0\left(t\right)=0,$$(21)

$${\dot{a}}_2\left(t\right)-\frac{2a_2\left(t\right)}{k{a}_0^3\left(t\right)}{A}_0\left(t\right)+\left[\frac{9}{k^3{a}_0^4\left(t\right)}+\frac{2}{k{a}_0^2\left(t\right)}\right]A_2\left(t\right)=0,$$(22)

$$\frac{{\dot{A}}_0\left(t\right)}{k{a}_0\left(t\right)}+\omega A_0\left(t\right)+\frac{1}{2}{\dot{a}}_0^2\left(t\right)+\frac{{\dot{B}}_0\left(t\right)}{L}=-\frac{\Delta P_1}{\rho },$$(23)

$$\frac{\left[3+2k^2{a}_0^2\left(t\right)\right]{\dot{A}}_2\left(t\right)}{k^3{a}_0^3\left(t\right)}-\omega A_2\left(t\right)-\frac{a_2\left(t\right){\dot{A}}_0\left(t\right)}{k{a}_0^2\left(t\right)}+{\dot{a}}_0\left(t\right){\dot{a}}_2\left(t\right)+\frac{20a_0\left(t\right){\dot{a}}_0\left(t\right)}{L^5}{B}_2\left(t\right)+\frac{10a_0^2\left(t\right)}{L^5}{\dot{B}}_2\left(t\right)=-\frac{f_{\mathrm{e}}}{\rho }.$$(24)

At the interface of bubble β, the total velocity potential satisfies the following equation [20]:

$${\frac{\mathrm{\partial }\mathrm{\Phi }}{\mathrm{\partial }{r}_{\beta }}|}_{r_{\beta }=S_{\beta }}={\frac{\mathrm{\partial }{\mathrm{\Phi }}_{\beta }}{\mathrm{\partial }{r}_{\beta }}|}_{r_{\beta }=S_{\beta }}=\frac{\mathrm{d}{S}_{\beta }}{\mathrm{d}t}.$$(25)

Substituting equations (1b) and (10)–(12) into equation (25), using the Taylor series expansion of ϵ, and neglecting the Peano remainder $O\left({ϵ}^2\right)$, we obtain

$${\dot{b}}_0\left(t\right)+\frac{1}{b_0^2\left(t\right)}{B}_0\left(t\right)=0,$$(26)

$${\dot{b}}_2\left(t\right)-\frac{2b_2\left(t\right)}{b_0^3\left(t\right)}{B}_0\left(t\right)+\frac{3}{b_0^4\left(t\right)}{B}_2\left(t\right)=0.$$(27)

Here, the liquid’s compressibility is not considered.

Let R_α0 and R_β0 be initial radii of bubble α and bubble β, namely R_α0 = a₀(0) and R_β0 = b₀(0). After eliminating A₀(t), A₂(t), B₀(t), B₂(t), and their derivatives from the equations consisted by equations (21)–(24) and equations (26)–(27), we obtain the model for bubble α under the action of bubble β:

$$a_0\left(t\right){\ddot{a}}_0\left(t\right)+\frac{3}{2}{\dot{a}}_0^2\left(t\right)+\mathrm{k\omega }{a}_0^2\left(t\right){\dot{a}}_0\left(t\right)+\frac{b_0\left(t\right)}{L}{I}_{b0}=\frac{\Delta P_1}{\rho },$$(28a)

$$\frac{N}{M}{a}_0\left(t\right){\ddot{a}}_2\left(t\right)+Z_{\beta }{\dot{a}}_2\left(t\right)+Z_{\mathrm{\varpi }}{a}_2\left(t\right)+\frac{20a_0\left(t\right){\dot{a}}_0\left(t\right)}{3L^5}{I}_{b1}+\frac{10a_0^2\left(t\right)b_0^2\left(t\right)}{3L^5}{I}_{b2}=\frac{f_{\mathrm{e}}}{\rho },$$(28b)

where

$$\{\begin{array}{ll}& k=\frac{\omega }{c}\\ & M=9+2k^2{a}_0^2\left(t\right)\\ & N=3+2k^2{a}_0^2\left(t\right)\\ & Z_{\beta }=\left[\frac{6N-M}{M}-\frac{4k^2{a}_0^2\left(t\right)N}{M^2}\right]{\dot{a}}_0\left(t\right)-\frac{k^3\omega a_0^4\left(t\right)}{M}\\ & Z_{\varpi }=\left(\frac{6N}{M}-2\right)\frac{{\dot{a}}_0^2\left(t\right)}{a_0\left(t\right)}+\left(\frac{2N}{M}-1\right){\ddot{a}}_0\left(t\right)\\ & -\frac{8k^2{a}_0\left(t\right){\dot{a}}_0^2\left(t\right)N}{M^2}-\frac{2k^3\omega a_0^3\left(t\right){\dot{a}}_0\left(t\right)}{M}\\ & I_{b0}=2{\dot{b}}_0^2\left(t\right)+b_0\left(t\right){\ddot{b}}_0\left(t\right)\\ & I_{b1}=b_0^3\left(t\right)\left[2{\dot{b}}_0\left(t\right)b_2\left(t\right)+b_0\left(t\right){\dot{b}}_2\left(t\right)\right]\\ & I_{b2}=6{\dot{b}}_0^2\left(t\right)b_2\left(t\right)+2b_0\left(t\right){\ddot{b}}_0\left(t\right)b_2\left(t\right)\\ & +6b_0\left(t\right){\dot{b}}_0\left(t\right){\dot{b}}_2\left(t\right)+b_0^2\left(t\right){\ddot{b}}_2\left(t\right)\\ & \Delta P_1=\left(p_0+\frac{2\sigma }{R_{\alpha 0}}-p_{\mathrm{v}}\right){\left[\frac{R_{\alpha 0}}{a_0\left(t\right)}\right]}^{3\gamma }+p_{\mathrm{v}}-p_0\\ & -\frac{2\sigma }{a_0\left(t\right)}-\frac{4\eta {\dot{a}}_0\left(t\right)}{a_0\left(t\right)}-p_{\mathrm{a}}\cos \omega t\\ & f_{\mathrm{e}}=\left\{1-{\left[\frac{R_{\alpha 0}}{a_0\left(t\right)}\right]}^{3\gamma }\right\}{\epsilon }_0\frac{{\left({\epsilon }_{\mathrm{r}\mathrm{l}}-1\right)}^2}{3}{E}^2.\end{array}$$

Equation (28a) is a modified Rayleigh-Plesset equation for a₀(t). The term $\mathrm{k\omega }{a}_0^2\left(t\right){\dot{a}}_0\left(t\right)$ is related to the oscillation resistance induced by the liquid’s compressibility. Equation (28b) can be regarded as a damped vibration equation for a₂(t) with external forces, without consideration of time-varying coefficients. In addition, all terms for L are related to the forces that bubble β exerts on bubble α.

Following similar steps, we can also obtain the model for bubble β under the action of bubble α, which has a similar form as the model for bubble α.

$$b_0\left(t\right){\ddot{b}}_0\left(t\right)+\frac{3}{2}{\dot{b}}_0^2\left(t\right)+\mathrm{k\omega }{b}_0^2\left(t\right){\dot{b}}_0\left(t\right)+\frac{a_0\left(t\right)}{L}{I}_{a0}=\frac{\Delta P_1^{\mathrm{\prime }}}{\rho },$$(29a)

$$\frac{N^{\mathrm{\prime }}}{M^{\mathrm{\prime }}}{b}_0\left(t\right){\ddot{b}}_2\left(t\right)+Z_{\beta }^{\mathrm{\prime }}{\dot{b}}_2\left(t\right)+Z_{\mathrm{\varpi }}^{\mathrm{\prime }}{b}_2\left(t\right)+\frac{20b_0\left(t\right){\dot{b}}_0\left(t\right)}{3L^5}{I}_{a1}+\frac{10b_0^2\left(t\right)a_0^2\left(t\right)}{3L^5}{I}_{a2}=\frac{f_{\mathrm{e}}^{\mathrm{\prime }}}{\rho },$$(29b)

where

$$\{\begin{array}{ll}& M^{\mathrm{\prime }}=9+2k^2{b}_0^2\left(t\right)\\ & N^{\mathrm{\prime }}=3+2k^2{b}_0^2\left(t\right)\\ & Z_{\beta }^{\mathrm{\prime }}=\left[\frac{6N^{\mathrm{\prime }}-M^{\mathrm{\prime }}}{M^{\mathrm{\prime }}}-\frac{4k^2{b}_0^2\left(t\right)N^{\mathrm{\prime }}}{{\left(M^{\mathrm{\prime }}\right)}^2}\right]{\dot{b}}_0\left(t\right)-\frac{k^3\omega b_0^4\left(t\right)}{M^{\mathrm{\prime }}}\\ & Z_{\varpi }^{\mathrm{\prime }}=\left(\frac{6N^{\mathrm{\prime }}}{M^{\mathrm{\prime }}}-2\right)\frac{{\dot{b}}_0^2\left(t\right)}{b_0\left(t\right)}+\left(\frac{2N^{\mathrm{\prime }}}{M^{\mathrm{\prime }}}-1\right){\ddot{b}}_0\left(t\right)\\ & -\frac{8k^2{b}_0\left(t\right){\dot{b}}_0^2\left(t\right)N^{\mathrm{\prime }}}{{\left(M^{\mathrm{\prime }}\right)}^2}-\frac{2k^3\omega b_0^3\left(t\right){\dot{b}}_0\left(t\right)}{M^{\mathrm{\prime }}}\\ & I_{a0}=2{\dot{a}}_0^2\left(t\right)+a_0\left(t\right){\ddot{a}}_0\left(t\right)\\ & I_{a1}=a_0^3\left(t\right)\left[2{\dot{a}}_0\left(t\right)a_2\left(t\right)+a_0\left(t\right){\dot{a}}_2\left(t\right)\right]\\ & I_{a2}=6{\dot{a}}_0^2\left(t\right)a_2\left(t\right)+2a_0\left(t\right){\ddot{a}}_0\left(t\right)a_2\left(t\right)\\ & +6a_0\left(t\right){\dot{a}}_0\left(t\right){\dot{a}}_2\left(t\right)+a_0^2\left(t\right){\ddot{a}}_2\left(t\right)\\ & \Delta P_1^{\mathrm{\prime }}=\left(p_0+\frac{2\sigma }{R_{\beta 0}}-p_{\mathrm{v}}\right){\left[\frac{R_{\beta 0}}{b_0\left(t\right)}\right]}^{3\gamma }+p_{\mathrm{v}}-p_0\\ & -\frac{2\sigma }{b_0\left(t\right)}-\frac{4\eta {\dot{b}}_0\left(t\right)}{b_0\left(t\right)}-p_{\mathrm{a}}\cos \omega t\\ & f_{\mathrm{e}}^{\mathrm{\prime }}=\left\{1-{\left[\frac{R_{\beta 0}}{b_0\left(t\right)}\right]}^{3\gamma }\right\}{\epsilon }_0\frac{{\left({\epsilon }_{\mathrm{r}l}-1\right)}^2}{3}{E}^2.\end{array}$$

When c →∞, namely k → 0, equations (28a)–(28b) and equations (29a)–(29b) are reduced to the equations for an incompressible liquid,

$$a_0\left(t\right){\ddot{a}}_0\left(t\right)+\frac{3}{2}{\dot{a}}_0^2\left(t\right)+\frac{b_0\left(t\right)}{L}{I}_{b0}=\frac{\Delta P_1}{\rho },$$(30a)

$$\frac{1}{3}{a}_0\left(t\right){\ddot{a}}_2\left(t\right)+{\dot{a}}_0\left(t\right){\dot{a}}_2\left(t\right)-\frac{1}{3}{\ddot{a}}_0\left(t\right)a_2\left(t\right)+\frac{20a_0\left(t\right){\dot{a}}_0\left(t\right)}{3L^5}{I}_{b1}+\frac{10a_0^2\left(t\right)b_0^2\left(t\right)}{3L^5}{I}_{b2}=\frac{f_{\mathrm{e}}}{\rho },$$(30b)

$$b_0\left(t\right){\ddot{b}}_0\left(t\right)+\frac{3}{2}{\dot{b}}_0^2\left(t\right)+\frac{a_0\left(t\right)}{L}{I}_{a0}=\frac{\Delta P_1^{\mathrm{\prime }}}{\rho },$$(31a)

$$\frac{1}{3}{b}_0\left(t\right){\ddot{b}}_2\left(t\right)+{\dot{b}}_0\left(t\right){\dot{b}}_2\left(t\right)-\frac{1}{3}{\ddot{b}}_0\left(t\right)b_2\left(t\right)+\frac{20b_0\left(t\right){\dot{b}}_0\left(t\right)}{3L^5}{I}_{a1}+\frac{10b_0^2\left(t\right)a_0^2\left(t\right)}{3L^5}{I}_{a2}=\frac{f_{\mathrm{e}}^{\mathrm{\prime }}}{\rho }.$$(31b)

When L → ∞, equations (30a) and (30b) can be further reduced to the equations for a single nonspherical bubble, see equations (32a) and (32b).

$$a_0\left(t\right){\ddot{a}}_0\left(t\right)+\frac{3}{2}{\dot{a}}_0^2\left(t\right)=\frac{\Delta P_1}{\rho },$$(32a)

$$\frac{1}{3}{a}_0\left(t\right){\ddot{a}}_2\left(t\right)+{\dot{a}}_0\left(t\right){\dot{a}}_2\left(t\right)-\frac{1}{3}{\ddot{a}}_0\left(t\right)a_2\left(t\right)=\frac{f_{\mathrm{e}}}{\rho }.$$(32b)

3 Numerical solutions and analysis

3.1 Effect of the distance between bubble centers on the radius of bubble α

By solving equations (28a)–(28b) and equations (29a)–(29b) numerically, we obtain the temporal evolutions of the spherical coefficient a₀(t) and the nonspherical coefficient a₂(t) of bubble α at different distances between bubble centers (see Figs. 2 and 3). In these two figures, the condition L → ∞ represents the oscillation of an isolated nonspherical bubble. The other parameters see Table 1.

Figure 2 Temporal evolution of a0(t) at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm).

Figure 3 (a) Temporal evolution of a2(t) at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm); (b) Partial enlarged detail of (a).

In Figure 2, the variations of a₀(t) with time, at different distances L, have a similar pattern. Within one period, a₀(t) first increases rapidly to its maximum during the rarefaction cycle, then, it decreases sharply to its minimum in a rather short time. This is followed by a series of rebounds. Furthermore, as the distance between two bubbles increases, the maximum of a₀(t) increases, and, correspondingly, the minimum of a₀(t) decreases. When L → ∞, a₀(t) reaches its highest maximum and lowest minimum. In Figure 3, a₂(t) shows exponential increases with time during the first period. Theoretically, a₂(t) shows a larger increase at a longer distance L. In other words, when the distance between bubble centers is longer, bubble α is subject to lower oscillation resistance from bubble β, and appears to be more ellipsoidal. The similar time-varying curves of a₂(t) were also reported by Liang et al. [11, 21, 22] and Chen et al. [14].

3.2 Effect of the initial radius of bubble β on the radius of bubble α

When the distance is fixed at L = 350 μm, and the value of R_β0 is variable, we can solve the equations (28a)–(28b) and equations (29a)–(29b) numerically and obtain time-varying curves of two radius components of bubble α. All parameters are the same as in Section 3.1. As depicted in Figure 4, when R_β0 ≤ 3.5 μm, the wave peaks of a₀(t) increase, while R_β0 is increasing. However, as the initial radius R_β0 increases from 3.5 μm to 6 μm, the wave peaks of a₀(t) decrease significantly. Hence, the highest maximum and the lowest minimum of a₀(t) occur when R_β0 = 3.5 μm. Moreover, the spherical component a₀(t) undergoes faster oscillation damping for smaller R_β0 values. Figure 5 reveals that, when bubble β has a higher initial radius R_β0, the coefficient a₂(t) shows a lower increase during one oscillation period. This means that the nonspherical feature of bubble α is more significant for smaller R_β0 values.

Figure 4 Temporal evolution of a0(t) for different initial radii of bubble β (T = 50 μs, Rα0 = 4.5 μm).

Figure 5 Temporal evolution of a2(t) for different initial radii of bubble β (T = 50 μs, Rα0 = 4.5 μm).

3.3 Analyses of the cavitation of bubble α using approximate equations

These rapid increases of the coefficient a₂(t), as shown in Figures 3 and 5, are unlikely to occur physically. In fact, it is impossible for a₂(t) to increase without limits. This can be interpreted as follows: A nonspherical bubble has already collapsed before its nonspherical coefficient a₂(t) reaches such a large value. The variation curves of a₂(t) are only the results of the numerical solution to dynamic equations. Thus, to further investigate the effect of bubble β during the cavitation of bubble α, we introduce approximate equations. The method is to take the mean of electric stress f_s in half period, and approximately transform a nonspherical bubble into a spherical bubble. The reason for adopting this method is as follows: the presence of electric stress reduces the surface tension at a bubble wall, which causes a bubble to expand more during the rarefaction cycle. Therefore, it is more suitable to take the mean of electric stress f_s in half period, rather than the mean of the bubble radius S(θ, t) in half period. A nonspherical bubble should be approximately treated as a spherical bubble with smaller surface tension, instead of a spherical bubble with a larger initial radius.

The mean of f_s in half period is

$$\overline{f_{\mathrm{s}}}=\frac{2}{\pi }{\int }_0^{\frac{\pi }{2}}\mathrm{}{f}_{\mathrm{s}} \mathrm{d}\theta ={\epsilon }_0\frac{{\left({\epsilon }_{\mathrm{r}\mathrm{l}}-1\right)}^2}{12}{E}^2.$$(33)

Let R_α and R_β be radii of approximate spherical bubbles of bubble α and bubble β. By inserting the term $ϵ\overline{f_{\mathrm{s}}}$ into the right sides of equations (28a) and (29a) and replacing a₀(t) and b₀(t) with R_α and R_β, we can rewrite equations (28a) and (29a) as

$$R_{\alpha }{\ddot{R}}_{\alpha }+\frac{3}{2}{\dot{R}}_{\alpha }^2+\mathrm{k\omega }{R}_{\alpha }^2{\dot{R}}_{\alpha }+\frac{R_{\beta }}{L}{I}_{\mathrm{R\beta }}=\frac{\Delta P_{\mathrm{M}}}{\rho },$$(34a)

$$R_{\beta }{\ddot{R}}_{\beta }+\frac{3}{2}{\dot{R}}_{\beta }^2+\mathrm{k\omega }{R}_{\beta }^2{\dot{R}}_{\beta }+\frac{R_{\alpha }}{L}{I}_{\mathrm{R\alpha }}=\frac{\Delta P_{\mathrm{M}}^{\mathrm{\prime }}}{\rho },$$(34b)

where

$$\{\begin{array}{ll}& I_{\mathrm{R\beta }}=2{\dot{R}}_{\beta }^2+R_{\beta }{\ddot{R}}_{\beta }\\ & I_{\mathrm{R\alpha }}=2{\dot{R}}_{\alpha }^2+R_{\alpha }{\ddot{R}}_{\alpha }\\ & \Delta P_{\mathrm{M}}=\left(p_0+\frac{2\sigma }{R_{\alpha 0}}-p_{\mathrm{v}}-ϵ\overline{f_{\mathrm{s}}}\right){\left(\frac{R_{\alpha 0}}{R_{\alpha }}\right)}^{3\gamma }+p_{\mathrm{v}}\\ & -p_0-\frac{2\sigma }{R_{\alpha }}-\frac{4\eta {\dot{R}}_{\alpha }}{R_{\alpha }}+ϵ\overline{f_{\mathrm{s}}}-p_a\cos \omega t\\ & \Delta P_{\mathrm{M}}^{\mathrm{\prime }}=\left(p_0+\frac{2\sigma }{R_{\beta 0}}-p_{\mathrm{v}}-ϵ\overline{f_{\mathrm{s}}}\right){\left(\frac{R_{\beta 0}}{R_{\beta }}\right)}^{3\gamma }+p_{\mathrm{v}}\\ & -p_0-\frac{2\sigma }{R_{\beta }}-\frac{4\eta {\dot{R}}_{\beta }}{R_{\beta }}+ϵ\overline{f_{\mathrm{s}}}-p_{\mathrm{a}}\cos \omega t.\\ & \end{array}$$

The equations (34a) and (34b) are the approximate equations that describe the radius variations of bubble α and bubble β with the interaction between two bubbles. To simplify the calculation, we assume that ϵ = 0.1 in the following analysis. As long as ϵ > 0, all conclusions below do not change with the ϵ value.

3.3.1 Effect of bubble β on the radius of bubble α

In Section 3.2, we concluded that, when R_β0 = 3.5 μm, a₀(t) reaches its highest maximum and lowest minimum. Thus, the numerical solutions of equations (34a) and (34b) are performed for the bubble pair with initial radii R_α0 = 4.5 μm and R_β0 = 3.5 μm. All other parameters are the same as in Section 3.1. Figure 6 shows the temporal evolution of R_α at different distances between two bubbles. The condition L → ∞ reflects the oscillation of an isolated bubble. With increasing distance between bubble centers, the oscillation amplitude of the radius of bubble α increases. This can be explained as follows: A longer distance between bubble centers leads to less oscillation resistance and higher expansion during the rarefaction cycle for bubble α. During the compression cycle, bubble α with a larger maximum radius undergoes a more violent collapse, and it therefore reaches a lower minimum radius. Hence, the highest maximum radius and the lowest minimum radius are obtained during the oscillation of an isolated bubble.

Figure 6 Temporal evolution of Rα at different distances between bubble centers (L) (T = 50 μs, Rα0 =4.5 μm).

3.3.2 Effect of bubble β on the gas pressure and the temperature inside bubble α

When the van der Waals hard-core radius h is taken into account, the expressions for the gas pressure P_g and the temperature T_g in the center of bubble α are given by [23, 24]

$$\{\begin{array}{c}\begin{array}{c}{P}_{\mathrm{g}}=\left(p_0+\frac{2\sigma }{R_{\alpha 0}}-p_{\mathrm{v}}-ϵ\overline{f_{\mathrm{s}}}\right){\left(\frac{R_{\alpha 0}^3-h^3}{{R_{\alpha }}^3-h^3}\right)}^{\gamma }\\ \begin{array}{ll}& \\ T_{\mathrm{g}}& =T_0{\left(\frac{R_{\alpha 0}^3-h^3}{{R_{\alpha }}^3-h^3}\right)}^{\gamma -1},\end{array}\end{array}\end{array}$$

where T₀ is the ambient temperature, T₀ = 298.15 K. For bubble α with an initial radius R_α0 = 4.5 μm, h = R_α0/8.54 ≈ 0.53 μm. As shown in Figures 7 and 8, under all conditions, the maximum gas-pressure and temperature inside bubble α occur at the end of the collapse, when bubble α reaches its minimum radius. Furthermore, during the negative pressure phase of the ultrasound wave, the maximum internal pressure and temperature of an isolated bubble are significantly higher, compared to a bubble with another adjacent bubble. In general, regardless of the distance between two bubbles and the initial radius of an adjacent bubble, the interaction between two bubbles always acts as resistance to the oscillation of a bubble, when parameters of ultrasound and an electrostatic field are constant.

Figure 7 Temporal evolution of Pg at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm, p0 = 1.013 × 105 Pa).

Figure 8 Temporal evolution of Tg at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm, T0 = 298.15 K).

4 Conclusion

In the present study, we developed a model for two interacting nonspherical bubbles in a compressible liquid under the coupling effect of ultrasound and electrostatic field. When the distance between bubble centers is shorter, bubble α is subject to higher oscillation resistance from bubble β, and appears to be more ellipsoidal. When the initial radius R_β0 equals to 3.5 μm, a₀(t) reaches its highest maximum and lowest minimum. For a larger R_β0 value, the nonspherical feature of bubble α becomes more significant after several periods. By solving approximate equations numerically, we found that, at the end of the collapse, the maximum pressure and temperature inside an isolated bubble are significantly higher than a bubble with an adjacent bubble. However, this obtained model can only depict the interaction between two adjacent bubbles. In fact, when a number of bubbles undergo ultrasonic cavitation simultaneously, the interaction between them can be more complicated. Therefore, a subject of our future research could be the oscillation of a nonspherical bubble cluster, as well as the effect of the interaction on the resonance frequency of the investigated bubble.

Acknowledgments

This research were supported by Provincial Major Project in Basic and Applied Research of Natural Science (Department of Education of Guangdong Province) (No.2017GKZDXM015), National Natural Science Foundation of China (No.11464002), Innovation Team Project of Colleges and Universities in Guangdong Province (Natural Science) (No.2020KCXTD054), Key Scientific Research Platform of Colleges and Universities (Department of Education of Guangdong Province) (No.2021GCZX017), Special Project in Key Fields of Universities in Guangdong Province (Science and Technology Services for Rural Revitalization) (No.2022ZDZX4118), and Youth Innovation Talent Project of Universities in Guangdong Province (No.2022KQNCX261).

Data Availability Statement

The research data associated with this article are included within the article.

References

All Tables

All Figures

	Figure 1 Sections of two nonspherical bubbles.
In the text

	Figure 2 Temporal evolution of a0(t) at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm).
In the text

	Figure 3 (a) Temporal evolution of a2(t) at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm); (b) Partial enlarged detail of (a).
In the text

	Figure 4 Temporal evolution of a0(t) for different initial radii of bubble β (T = 50 μs, Rα0 = 4.5 μm).
In the text

	Figure 5 Temporal evolution of a2(t) for different initial radii of bubble β (T = 50 μs, Rα0 = 4.5 μm).
In the text

	Figure 6 Temporal evolution of Rα at different distances between bubble centers (L) (T = 50 μs, Rα0 =4.5 μm).
In the text

	Figure 7 Temporal evolution of Pg at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm, p0 = 1.013 × 105 Pa).
In the text

	Figure 8 Temporal evolution of Tg at different distances between bubble centers (L) (T = 50 μs, Rα0 = 4.5 μm, T0 = 298.15 K).
In the text