Abstract

The homotopy analysis method is used to obtain analytical solutions of the Rayleigh equation for the radial oscillations of a multielectron bubble in liquid helium. The small order approximations for amplitude and frequency fit well with those computed numerically. The results confirm that the homotopy analysis method is a powerful and manageable tool for finding analytical solutions of strongly nonlinear dynamical systems.

1. Introduction

Nonlinear equations are widely used for modeling complex phenomena in various fields of sciences and engineering. Nonlinear problems are in general difficult to solve analytically. In recent years, the tool for finding analytical solutions of nonlinear equations known as the homotopy analysis method (HAM) has been developed by Liao [13]. This nonperturbation technique, independent of small/large physical parameters, allows to effectively control the convergence and accuracy of the series solution to the model under consideration. HAM has been applied successfully to many nonlinear problems such as free oscillations of self-excited systems [4], heat radiation [5], finding travelling-wave solutions of the Kawahara equation [6], finding solitary wave solutions for the fifth-order KdV equation [7], exponentially decaying boundary layers [8], free vibrations of tapered beams [9], and finding solutions for Duffing equations with cubic and quintic nonlinearities [10].

HAM method is applied in this paper for finding a periodic solution of the Rayleigh equation, which describes the radial free oscillations of a multielectron bubble (MEB) in liquid Helium [11, 12]: 𝑅5𝑑2𝑅𝑑𝑡2+32𝑅4𝑑𝑅𝑑𝑡21=𝜌2𝜎𝑅3𝑒2𝑍232𝜋2𝜀𝜀0,(1.1) subject to the initial conditions 𝑅(0)=𝜉𝑅𝑐,𝑑𝑅(0)𝑑𝑡=0,(1.2) where 𝑅 is the bubble radius, 𝜌 and 𝜎 are the mass density and surface tension of the liquid Helium, respectively, 𝑒 is the electron charge, 𝑍 is the number of electrons in the bubble, 𝜀 is the dielectric constant of helium, 𝜀0 is thevacuum permittivity, 𝑅𝑐 is the equilibrium Coulomb radius of the MEB, and 𝜉 is a dimensionless parameter.

Another aim of this work is to obtain analytical solutions of the Rayleigh equation (1.1) with an additional external force driving term, a pressure step function. With this kind of forcing, Tempere et al. [13, 14] have shown that the angularly undeformed MEB undergoes cyclic collapses and reexpansions with high frequency (in the order of MHz). During the last stage of collapse, the radial acceleration is sharp and large (greater than 106 m/s2) resulting in pulses of electromagnetic radiation or sonoluminescence.

2. Foundations of the HAM

The text sequence used in [9, 10] to show the basic ideas of the HAM is closely followed here. Consider a nonlinear differential equation expressed by the following: 𝑁[]𝑟(𝑡)=0,(2.1) where 𝑁 is a nonlinear differential operator, 𝑡 denotes the independent variable; 𝑟(𝑡) is an unknown function.

Applying the transformation 𝜏=𝜔𝑡, (2.1) can be expressed by the following: 𝑁[]𝑟(𝜏),𝜔=0.(2.2) Now a homotopy in general form can be constructed as follows: 𝐻(𝜙,𝑞,,𝐻(𝜏))=(1𝑞)𝐿𝜙(𝜏,𝑞)𝑟0[],(𝜏)𝑞𝐻(𝜏)𝑁𝜙(𝜏,𝑞),𝜔(𝑞)(2.3) where 𝑞[0,1] is an embedding parameter, 𝜙 is a function of 𝜏 and 𝑞, is a nonzero auxiliary parameter, 𝐻(𝜏) is a nonzero auxiliary function, and 𝐿 denotes an auxiliary linear operator. When the parameter 𝑞 increases from 0 to 1, 𝜔(𝑞) varies from 𝜔0 to 𝜔, and the solution 𝜙(𝜏,𝑞) varies from the initial approximation 𝑟0(𝜏) to the exact solution 𝑟(𝜏). Put differently, 𝜙(𝜏,0)=𝑟0(𝜏) is the solution of 𝐻(𝜙,𝑞,,𝐻(𝜏))|𝑞=0=0, and 𝜙(𝜏,1)=𝑟(𝜏) is the solution of 𝐻(𝜙,𝑞,,𝐻(𝜏))|𝑞=1=0. Setting 𝐻(𝜙,𝑞,,𝐻(𝜏))|=0, the zero-order deformation equation is obtained as 𝜙(1𝑞)𝐿(𝜏,𝑞)𝑟0[𝜙](𝜏)=𝑞𝐻(𝜏)𝑁(𝜏,𝑞),𝜔(𝑞),(2.4) under the initial conditions: 𝜙(0,𝑞)=𝑘,𝑑𝜙(0,𝑞)𝑑𝜏=0,(2.5) where 𝑘 is a constant. The functions 𝜙(𝜏,𝑞) and 𝜔(𝑞) can be expanded as power series of 𝑞 using Taylor’s theorem as follows 𝜙(𝜏,𝑞)=𝜙(𝜏,0)+𝑚=11𝜕𝑚!𝑚𝜙(𝜏,𝑞)𝜕𝑞𝑚|||||𝑞=0𝑞𝑚=𝑟0(𝜏)+𝑞=1𝑟𝑚(𝜏)𝑞𝑚,𝜔(𝑞)=𝜔0+𝑚=11𝜕𝑚!𝑚𝜔(𝑞)𝜕𝑞𝑚|||||𝑞=0𝑞𝑚=𝜔0+𝑞=1𝜔𝑚𝑞𝑚,(2.6) where 𝑟𝑚(𝜏) and 𝜔𝑚 are called the mth-order deformation derivatives.

Differentiating the zero-order deformation equation with respect to 𝑞 and then setting 𝑞=0, yields the first-order deformation equation which gives the first-order approximation of 𝑟(𝜏) as follows: 𝐿𝑟1𝑟(𝜏)=𝐻(𝜏)𝑁0(𝜏),𝜔0,(2.7) subject to the initial conditions: 𝑟1(0)=0,𝑑𝑟1(0)𝑑𝜏=0.(2.8) Differentiating Equations (2.4) and (2.5) 𝑚 times with respect to 𝑞, then setting 𝑞=0 and finally divided them by m!, yields the so called 𝑚th-order (𝑚>1) deformation equation: 𝐿𝑟𝑚(𝜏)𝜒𝑚𝑟𝑚1(𝜏)=𝐻(𝜏)𝐷𝑚𝑟𝑚1,𝜔𝑚1,(2.9) subject to the initial conditions 𝑟𝑚(0)=0,𝑑𝑟𝑚(0)𝑑𝜏=0,(2.10) where 𝐷𝑚(𝑟𝑚1,𝜔𝑚1), 𝑟𝑚1 and 𝜔𝑚1 are defined as follows: 𝐷𝑚𝑟𝑚1,𝜔𝑚1=1𝑑(𝑚1)!𝑚1𝑁[]𝜙(𝜏,𝑞),𝜔(𝑞)𝑑𝑞𝑚1||||𝑞=0,(2.11)𝑟𝑚1=𝑟0,𝑟1,𝑟2,,𝑟𝑚1,(2.12)𝜔𝑚1=𝜔0,𝜔1,𝜔2,,𝜔𝑚1,𝜒(2.13)𝑚=0,𝑚1,1,𝑚>1.(2.14)

3. Application of the HAM

3.1. Radial Free Oscillations

A nondimensional radius 𝑅 and time 𝑡 are defined as follows: 𝑅=𝑅[𝑅],𝑡=𝑡[𝑡],(3.1) where [𝑅] and [𝑡] are characteristic scales.

In terms of the dimensionless variables (3.1), the equation of motion (1.1) becomes 𝑅5𝑑2𝑅𝑑𝑡2+32𝑅4𝑑𝑅𝑑𝑡2=𝜃𝛼𝑅3,(3.2) and the initial conditions (1.2) become 𝑅𝑑(0)=𝐴,𝑅(0)𝑑𝑡=0,(3.3) where α, 𝜃, and 𝐴 are dimensionless parameters given by the following: [𝑡]𝛼=2𝜎2𝜌[𝑅]3𝑒,𝜃=2𝑍2[𝑡]232𝜋2𝜀𝜀0𝜌[𝑅]6,𝐴=𝜉𝑅𝑐[𝑅].(3.4) Introducing the new independent variable 𝜏=𝜔𝑡, where 𝜔 is the dimensionless natural frequency, (3.2) and the initial conditions (3.3) can be rewritten as follows: 𝜔2𝑅5𝑑2𝑅𝑑𝜏2+32𝜔2𝑅4𝑑𝑅𝑑𝜏2=𝜃𝛼𝑅3,(3.5)𝑅𝑑(0)=𝜉,𝑅(0)𝑑𝜏=0.(3.6) In order to reduce the number of recursions and save computational efforts, this will be clarified later, an additional transformation is introduced 𝑅(𝜏)=3𝑟(𝜏),(3.7) then, (3.5) becomes 13𝜔2𝑟𝑑2𝑟𝑑𝜏21𝜔182𝑑𝑟𝑑𝜏2+𝛼𝑟𝜃=0,(3.8) under the initial conditions 𝑟(0)=𝜉3,𝑑𝑟(0)𝑑𝜏=0.(3.9) The free oscillation of the conservative system represented by (3.8) and (3.9) is a periodic motion, which can be expressed by the following base functions [1]: {cos(𝑚𝜏)𝑚0},(3.10) in the compact form 𝑟(𝜏)=+𝑚=0𝑐𝑚cos(𝑚𝜏).(3.11) With the purpose of satisfying the initial conditions (3.9), the initial guess of 𝑟(𝜏) for zero-order deformation equation is chosen as follows: 𝑟0𝜉(𝜏)=31cos𝜏+1.(3.12) Under the Rule of Solution Expression denoted by (3.11), the linear operator is selected as: 𝐿[]𝑟(𝜏;𝑞)=𝜔20𝜕2𝑟(𝜏;𝑞)𝜕𝑟2,+𝑟(𝜏;𝑞)(3.13) with the property 𝐿𝐶1sin𝜏+𝐶2cos𝜏=0,(3.14) where 𝐶1 and 𝐶2 are the integral constants.

From (3.8), the nonlinear operator is written as 𝑁[]=1𝑟(𝜏;𝑞),𝜔(𝑞)3𝜔(𝑞)2𝜕𝑟(𝜏;𝑞)2𝑟(𝜏;𝑞)𝜕𝜏2118𝜔(𝑞)2𝜕𝑟(𝜏;𝑞)𝜕𝜏2+𝑟(𝜏;𝑞)𝜃.(3.15) For complying with the general form of the base functions (3.10), the auxiliary function 𝐻(𝜏) is chosen as follows: 𝐻(𝜏)=1.(3.16) Substituting (3.15) into (2.11) and performing the homotopy-derivatives [3], yields: 𝐷𝑚1(𝜏)=3𝑚1𝑖=0𝑟𝑚1𝑖(𝜏)𝑖𝑗=0𝑟𝑖𝑗(𝜏)𝑗=0𝜔𝑗𝜔118𝑚1𝑖=0𝑟𝑚1𝑖(𝜏)𝑖𝑗=0𝑟𝑖𝑗(𝜏)𝑗=0𝜔𝑗𝜔+𝛼𝑟𝑚1(𝜏)𝜃1𝜒𝑚,(3.17) where a prime indicates first derivative with respect to 𝜏.

At this point, it is worth comparing (3.17) with the 𝑚th-order homotopy derivative corresponding with (3.5), which reads as follows: 𝐷𝑚(𝜏)=𝑚1𝑖=0𝑅𝑖𝑚1𝑖𝑗=0𝑅𝑗𝑖𝑗=0𝑅𝑗𝑚=0𝑅𝑚𝑚𝑝=0𝑅𝑝𝑚𝑝𝑞=0𝑅𝑞𝑝𝑞𝑟=0𝜔𝑞𝑟𝜔𝑟+32𝑚1𝑖=0𝑅𝑖𝑚1𝑖𝑗=0𝑅𝑗𝑖𝑗=0𝑅𝑗𝑚=0𝑅𝑚𝑚𝑝=0𝑅𝑝𝑚𝑝𝑞=0𝑅𝑞𝑝𝑞𝑟=0𝜔𝑞𝑟𝜔𝑟+𝛼𝑚1𝑖=0𝑅𝑖𝑚1𝑖𝑗=0𝑅𝑖𝑗𝑅𝑗𝜃1𝜒𝑚,(3.18) Clearly, the number of recursions for computing (3.18) is bigger than that using (3.17). Thus, this last equation is used to save long calculations.

According to the property (3.14) of the auxiliary linear operator 𝐿, if the term cos(𝜏) exists in 𝐷𝑚(𝜏), the so-called secular term 𝜏sin(𝜏) will appear in the final solution disobeying the rule solution expression (3.11). To avoid the secular terms, the coefficient of cos(𝜏) in 𝐷𝑚(𝜏) must be zero.

From (3.17), the first-order approximation (𝑚=1) of HAM for 𝐷1(𝜏) is as follows: 𝐷1=7𝛼𝜃𝜔3620+14𝜉363𝜔207𝜉366𝜔20+𝜉31𝛼𝛼+3𝜔2013𝜉3𝜔20+cos(𝜏)10𝜉363𝜔205𝜔36205𝜉366𝜔20cos(2𝜏).(3.19) Thus, 𝜔0 results in: 𝜔0=3𝛼.(3.20) Solving (2.7) subject to conditions (2.8) 𝑟1(𝜏) is obtained as: 𝑟1(𝜏)=sin2(𝜏/2)54(𝜃𝛼)+𝜔20𝜉312(13+5cos(𝜏))27𝜔20.(3.21) Similarly, from the coefficient of cos(𝜏) in 𝑅2(𝜏), 𝜔1 is obtained as follows: 𝜔1=1296𝜔30𝜉311944𝛼𝜃1944𝛼2+1584𝛼𝜔201224𝜉3𝛼𝜔20+288𝜉6𝛼𝜔201296𝜃𝜔20+648𝜉3𝜃𝜔20137𝜔40+315𝜉3𝜔40219𝜉6𝜔40+41𝜉9𝜔40.(3.22) Solving (2.9) under conditions (2.10) for 𝑚=2, 𝑟2(𝜏) results in: 𝑟21(𝜏)=11664𝜔40×sin2𝜏223328𝛼(𝛼𝜃)+𝜔20144Γ8+Γ21𝜔0Γ7𝜔0+864Γ6𝜔1,(3.23) where Γ1=𝜉3,Γ12=Γ129𝜉3,Γ1071623=27+13Γ1,Γ4=𝛼5+𝜉3Γ6𝜃,5=2167+103𝜉3,Γ6Γ=(13+5cos(𝜏)),7=56161379+1859𝜉3+10Γ5cos(𝜏)+195Γ1,Γcos(2𝜏)8=Γ2𝛼+6Γ3𝜃5Γ1Γ4.cos(𝜏)(3.24) The higher-order approximations of 𝑟𝑚(𝜏) and 𝜔𝑚1(𝜏), for 𝑚>2, are calculated similarly.

3.2. Driven Oscillations

The Rayleigh equation (1.1) with a driving pressure step 𝑝 is written as follows: 𝑅5𝑑2𝑅𝑑𝑡2+32𝑅4𝑑𝑅𝑑𝑡21=𝜌2𝜎𝑅3+𝑝𝑅4𝑒2𝑍232𝜋2𝜀𝜀0,(3.25) subject to the initial conditions 𝑅(0)=𝑅𝑐,𝑑𝑅(0)𝑑𝑡=0.(3.26) A nondimensional radius 𝑥 and time 𝑡 are defined as follows: 𝑅𝑥=[𝑅],𝑡=𝑡[𝑡].(3.27) In terms of the dimensionless variables (3.27), the equation of motion (3.25) becomes 𝑥5𝑑2𝑥𝑑𝑡2+32𝑥4𝑑𝑥𝑑𝑡2=𝜃𝛼𝑥3𝛽𝑥4,(3.28) and the initial conditions (3.26) become 𝑥(0)=𝑎,𝑑𝑥(0)𝑑𝑡=0,(3.29) where 𝛼, 𝛽, 𝜃, and 𝑎 are dimensionless parameters given by the following: [𝑡]𝛼=2𝜎2𝜌[𝑅]3𝑝[𝑡],𝛽=2𝜌[𝑅]2𝑒,𝜃=2𝑍2[𝑡]232𝜋2𝜀𝜀0𝜌[𝑅]6𝑅,𝑎=𝑐[𝑅].(3.30) Under the transformation 𝜏=𝜔𝑡, where 𝜔 is the dimensionless natural frequency, (3.28) and the initial conditions (3.29) can be rewritten as follows: 𝜔2𝑥5𝑑2𝑥𝑑𝜏2+32𝜔2𝑥4𝑑𝑥𝑑𝜏2+𝛽𝑥4+𝛼𝑥3𝜃=0,(3.31)𝑥(0)=𝑎,𝑑𝑥(0)𝑑𝜏=0.(3.32) Equation (3.31) can be integrated once [13], yielding: 𝜔2𝑥4𝑑𝑥𝑑𝜏2+23𝛽𝑥4+𝛼𝑥32𝛼+3𝛽+2𝜃𝑥+2𝜃=0.(3.33) The form of (3.33) is simpler than that of (3.31) implying less computations. In other words, the homotopy-derivative of the first term of (3.31) (i.e., 𝑚1𝑖=0𝑥𝑚1𝑖𝑖𝑗=0𝑥𝑖𝑗𝑗=0𝑥𝑗𝑚=0𝑥𝑚𝑚𝑝=0𝑥𝑚𝑝𝑝𝑞=0𝑥𝑝𝑞𝑞𝑟=0𝜔𝑞𝑟𝜔𝑟, which is at the order of magnitude 𝑛8) is no longer present in (3.33). Only remains the second term of (3.31) (without the numeric factor 3/2) whose homotopy-derivative (i.e., 𝑚1𝑖=0𝑥𝑚1𝑖𝑖𝑗=0𝑥𝑖𝑗𝑗=0𝑥𝑗𝑚=0𝑥𝑚𝑚𝑝=0𝑥𝑚𝑝𝑝𝑞=0𝑥𝑝𝑞𝑞𝑟=0𝜔𝑞𝑟𝜔𝑟) is at the order of magnitude 𝑛8. Even though the term (𝛼+(2/3)𝛽+2𝜃)𝑥 appears in (3.33), its homotopy-derivative (i.e., (𝛼+(2/3)𝛽+2𝜃)𝑥𝑚1) is only at the order of magnitude 𝑛1.

An additional simplification is found when (3.33) is substituted into (3.31), resulting in 𝜔2𝑥5𝑑2𝑥𝑑𝜏212𝛼𝑥3+32𝛼+𝛽+3𝜃𝑥4𝜃=0.(3.34) Now, the first term of (3.31) is recovered instead of the first term of (3.33). This is not advantageous because both terms involve the same number of recursions, but the second term in (3.33) and its corresponding recursions (i.e., (2/3)𝛽(𝑚1𝑖=0𝑥𝑚1𝑖𝑖𝑗=0𝑥𝑖𝑗𝑗=0𝑥𝑗𝑥)) are absent, thus saving computations.

From now on, the system represented by (3.34) and the initial conditions (3.32), is adopted to be solved by HAM.

The oscillation of the conservative system (3.34) and (3.32) after the application of the step pressure is a periodic motion which can be expressed by the following base functions: {cos(𝑚𝜏)𝑚0},(3.35) in the form 𝑥(𝜏)=+𝑚=0𝑐𝑚cos(𝑚𝜏).(3.36) It should be noted that the system (3.34) and (3.32) oscillates between the initial radius 𝑥(0)=𝑎 and the minimum radius 𝑥min, that is, when bubble wall velocity becomes zero (𝑑𝑥/𝑑𝜏)=0; then from (3.33): 𝑥min1=3𝛼+2𝛽(3𝛼4𝛽)(3𝛼+2𝛽)+Φ6𝛽Φ,6𝛽(3.37) where Φ is given by the following: Φ=324𝛽2𝜃(3𝛼7𝛽)(3𝛼+2𝛽)2+9𝛽8𝜃(3𝛼+2𝛽)2(7𝛽3𝛼)(𝛼2𝛽)(3𝛼+2𝛽)3+1296𝛽2𝜃21/21/3.(3.38) To satisfy the initial conditions (3.32), the initial guess of 𝑥(𝜏) for zero-order deformation equation is chosen as follows: 𝑥0(𝜏)=𝑎𝑥min2cos𝜏+𝑎+𝑥min2.(3.39) Under the Rule of Solution Expression denoted by (3.36), the linear operator for this case is identical to that chosen for free oscillations, that is, (3.13) and its corresponding property denoted by (3.14).

From (3.34), the nonlinear operator is defined by the following: 𝑁[]𝑥(𝜏;𝑞),𝜔(𝑞)=𝜔2(𝑞)𝑥5𝜕(𝜏;𝑞)2𝑥(𝜏;𝑞)𝜕𝜏212𝛼𝑥33(𝜏;𝑞)+2𝛼+𝛽+3𝜃𝑥(𝜏;𝑞)4𝜃.(3.40) The solution must comply with the general form of the base functions (3.35). Hence, the auxiliary function 𝐻(𝜏) must be chosen as follows: 𝐻(𝜏)=1.(3.41) Substituting (3.40) into (2.11) and performing the homotopy-derivatives, yields: 𝐷𝑚(𝜏)=𝑚1𝑖=0𝑥𝑖𝑚1𝑖𝑗=0𝑥𝑗𝑖𝑗=0𝑥𝑗𝑚=0𝑥𝑚𝑚𝑝=0𝑥𝑝𝑚𝑝𝑞=0𝑥𝑞𝑝𝑞𝑟=0𝜔𝑞𝑟𝜔𝑟12𝛼𝑚1𝑖=0𝑥𝑖𝑚1𝑖𝑗=0𝑥𝑖𝑗𝑥𝑗+32𝑥𝛼+𝛽+3𝜃𝑚14𝜃1𝜒𝑚.(3.42) As before, to avoid the secular terms, the coefficient of cos(𝜏) in 𝐷𝑚(𝜏) must be zero. Thus, the frequency 𝜔0 is obtained by setting to zero the coefficient of cos(𝜏) in 𝐷1(𝜏), yielding: 𝜔0=22𝑎𝑥min35𝑎216𝛼32(𝛽+3𝜃)+3𝛼𝑥min6𝑎+5𝑥min93𝑎668𝑎5𝑥min15𝑎4𝑥2min+15𝑎2𝑥4min+68𝑎𝑥5min93𝑥6min.(3.43) Then solving (2.7) under the conditions (2.8), 𝑥1(𝜏) is obtained as follows: 𝑥1(𝜏)=sin2(𝜏/2)215040𝜔201680𝐾1+𝑎6𝐾2𝜔20+5𝑎2𝐾3𝑥4min𝜔204𝑎𝐾4𝑥5min𝜔20+𝐾5𝑥6min𝜔204𝑥min×420𝐾6+𝑎5𝐾7𝜔20+5𝑎𝑥2min336𝛼𝐾8+𝑎3𝐾9𝜔20+80𝑥3min21𝛼𝐾1016𝑎3𝐾11𝜔20sin4𝜏2,(3.44) where 𝐾1=192𝑎𝛼+49𝑎3𝛼128𝑎𝛽+1024𝜃384𝑎𝜃+18𝑎3𝛼cos(𝜏)+𝑎3𝐾𝛼cos(2𝜏),2𝐾=71087+53974cos(𝜏)+9974cos(2𝜏)+1674cos(3𝜏)+199cos(4𝜏)+12cos(5𝜏),3𝐾=2255310cos(𝜏)+74cos(2𝜏)+38cos(3𝜏)103cos(4𝜏)+36cos(5𝜏),4𝐾=10013+5326cos(𝜏)1454cos(2𝜏)+516cos(3𝜏)139cos(4𝜏)+18cos(5𝜏),5𝐾=56737+25274cos(𝜏)4726cos(2𝜏)+974cos(3𝜏)151cos(4𝜏)+12cos(5𝜏),6=192𝛼13𝑎2𝛼+128𝛽+384𝜃+22𝑎2𝛼cos(𝜏)+3𝑎2𝐾𝛼cos(2𝜏),7𝐾=15963+17226cos(𝜏)+4846cos(2𝜏)+1216cos(3𝜏)+211cos(4𝜏)+18cos(5𝜏),8𝐾=1910cos(𝜏)+3cos(2𝜏),9𝐾=1345+1510cos(𝜏)+1334cos(2𝜏)+738cos(3𝜏)+274cos(4𝜏)+36cos(5𝜏),10𝐾=4714cos(𝜏)+cos(2𝜏),11=42+49cos(𝜏)+18cos(2𝜏)+3cos(3𝜏).(3.45) The higher order approximations for 𝜔𝑚1 and 𝑥𝑚(𝜏) when 𝑚>1 can be similarly calculated.

Two sets of results are presented. The first one is for free oscillations (Figures 17) and the second one is for driven oscillations (Figures 811).

The following characteristic scales were chosen to study the effect of varying both the initial radius of the MEB (by changing 𝜉) and the quantity of electrons in its surface (by changing 𝜃): [𝑅]=𝑅𝑐,[𝑡]=𝜌𝑅3𝑐,2𝜎(3.46) thus, the parameters (3.4) become 𝑒𝛼=1,𝜃=2𝑁264𝜋2𝜀𝜀0𝜎𝑅3𝑐,𝐴=𝜉.(3.47) Figure 1 shows the curves of 𝜔 as function of for different order approximations. The curves indicate that the valid regions of (i.e., the interval where 𝜔 versus is a horizontal line) are 3<<2 and 3<<1 for the 2nd and 5th-order approximations, respectively. Clearly, the convergence region of the series of the frequency 𝜔 increases with increasing the order approximation.

The 2nd-order approximation for a proper value of chosen in the valid interval 3<<2 is compared with the numerical solution by Runge-Kutta scheme as illustrated in Figure 2.

Solutions in Figure 2 correspond to a MEB containing 104 electrons, as that numerically simulated by Salomaa and Williams [11], accordingly its equilibrium Coulumb radius 𝑅𝑐 is about 1.064 microns. The initial radius of the MEB is 1.063 (precisely the value of 𝜉) times bigger than the equilibrium radius 𝑅𝑐; thus, after time zero, the bubble begins to oscillate with frequency around MHz. Clearly the 2nd-order HAM approximation is in excellent agreement with numerical solutions. This is quantified by computing the relative error for the 𝑀th-order approximation as follows ||RelativeError=𝑟(𝜏)𝑁𝑀𝑚=0𝑟𝑚||(𝜏)𝑟(𝜏)𝑁×100,(3.48) where 𝑟(𝜏)𝑁 is the solution by Runge-Kutta method.

The maximum relative error, indicated with a filled triangle in Figure 3, between 2nd-order approximation and numerical solution is about 0.129%. When increasing the order of the approximation, the maximum relative error becomes smaller; accordingly, the maximum relative error between 5th-order HAM approximation and numerical solution is reduced to 0.063% (indicated in Figure 3 with a hollow triangle), which is acceptably low.

The effect of incrementing the initial radius of the MEB is depicted in Figure 4. As before, a bubble containing 104 electrons is considered but now the initial dimensionless radius is fixed to 𝜉=1.15. Seemingly, the 2nd-order approximation in Figure 4(a) describes well the behavior of both the radius and the bubble wall velocity. Nevertheless, some deviation between HAM approximation and numerical solutions for the bubble wall acceleration is observed. This is overcome by increasing the order of approximation; consequently, the 5th-order HAM solution exhibits good agreement with the numerical solution as shown in Figure 4(b).

Figure 5 illustrates the first, second and fifth-order approximations of 𝜔 versus numerical results from Runge-Kutta method. As before, higher-order HAM solutions agree better with numerical results. The comparisons indicate that the 5th-order approximation agrees well with the numerical results even for a large initial dimensionless radius 𝜉=1.2. In fact, the maximum relative error between 5th-order approximation and numerical solution is about 0.066%; instead, the maximum error between the 2nd-order HAM solution and the numerical one is around 0.948%, as shown in Figure 6. It is interesting to note that the exponential growth tendency of the error curve for the 2nd-order approximation is dramatically reduced by increasing the order of HAM series. Clearly, the error curve for the 5th-order approximation shows slight upward trend with regard to that for the 2nd-order.

Figure 7 shows the validity of the 5th-order HAM solutions for 𝑟(𝜏) when compared with numerical integrations. The effect of varying the quantity of electrons in the MEB while fixing its initial radius 𝜉 is analyzed via changing the value of 𝜃. Figure 7(a) illustrates that smaller bubbles collapse more violently implying both smaller minimum radiuses and higher oscillation frequencies, for example, a small bubble containing 10 electrons attains at the first main collapse a minimum dimensionless radius 0.23 and develops an oscillation dimensionless frequency 𝜔0.33, whereas a big bubble containing 108 electrons attains at the first main collapse a minimum dimensionless radius 0.72 and develops a dimensionless frequency 𝜔0.27. Clearly, the analytical approximations for different values of 𝜃 closely approach the numerical results. Accordingly, in quantitative terms, the maximum relative errors (0.418, 0.154, 0.152, and 3.068%) for the respective values of 𝜃 (0.617, 0.833, 0.999, and 1.05) are acceptably low; even for the biggest bubble which contains 108 electrons.

In Figure 7(b), the robustness of the 5th-order HAM solution is analyzed for high values of the initial radius 𝜉 while fixing the quantity of electrons in the MEB. For 𝜉1.07, the bubble implodes reaching both high acceleration (≥1 × 106 m/s2) and a minimum radius more and more close to zero. This strong collapse conditions might involve sonoluminescence phenomenon as suggested by Tempere et al. [13]. Even for 𝜉 as high as 1.4 the analytical solution agrees to some extent with the numerical results, the maximum relative error in this case is around 5%. It is worthy of note that the parameter was varied to ensure the convergence of the analytical solution, the freedom of changing conveniently is precisely one of the virtues of the HAM [1, 6, 9].

In order to study the effect of varying the driving pressure 𝑝 in (3.25), the following scales where chosen: [𝑅]=143𝑒2𝑁2𝜋2𝜀𝜀0𝜎,[𝑡]=𝑒𝑁8𝜋𝜎𝜌2𝜀𝜀0,(3.49) then the parameters (3.30), become 𝑝𝛼=1,𝛽=83𝑒2𝑁2𝜋2𝜀𝜀0𝜎4,𝜃=1,𝑎=4𝑅𝑐3𝜋2𝜀𝜀0𝜎𝑒2𝑁2.(3.50) Figure 8 depicts comparisons between the 4th-order HAM solutions and numerical integrations of: 𝑥(𝜏), 𝑥(𝜏), and 𝑥(𝜏) for two values of the pressure step. Excellent agreement between analytical and numerical solutions is clearly seen for 𝑝=300 Pa (𝛽=0.443495); accordingly, the maximum relative error between numerical and analytical results for 𝑥(𝜏) is 0.304%, indicated by a hollow triangle in Figure 9(a). Nevertheless, when the pressure step is increased to 𝑝=1000 Pa (𝛽=1.47832), some disagreement between HAM approximations and numerical solutions of: 𝑥(𝜏), 𝑥(𝜏), and 𝑥(𝜏) is observed. In this case the maximum relative error for 𝑥(𝜏) is 4.664%, which might be acceptable if rough estimations are required. Even higher disagreement is clearly observed between HAM approximation and numerical integration for 𝑥(𝜏). This discrepancy might be reduced by incrementing the order of the HAM approximation; but this would imply larger and larger analytical solutions, and therefore more and more difficult to handle with. For the current example, an explicit expression for the 4th-order 𝑥(𝜏) solution would require several pages to write upon. As pointed out by Liao [3], in present time of computer with huge data storage capacities and high-speed CPU, an analytical expression with many terms might be accepted by most researchers. However, an important factor to be taken into account is the CPU time for evaluating an analytical solution with many terms (of course this time tends to increase together with the number of terms). The CPU time consumed in evaluating the 4th-order approximation for 𝑥(𝜏) is 12.933 s (computed by Timing command of Mathematica 5.1 in a Windows PC with Intel Core i7 processor at 3.4 GHz), which is around 275 times longer than time for computing a numerical solution of (3.27) and (3.28) for 𝑥(𝜏) using the fourth order Runge-Kutta algorithm (this time is computed by Timing and NDSolve commands of Mathematica 5.1 in a Windows PC with Intel Core i7 processor at 3.4 GHz). Of course, this result does not invalidate the 4th-order analytical approximations but it makes them unsuitable for certain situations. For instance, an analysis of shape instabilities of bubbles using the 4th-order approximations for 𝑥(𝜏), 𝑥(𝜏) and 𝑥(𝜏) would imply more computing time and less accuracy than those obtained by a numerical scheme.

Figure 10 illustrates the variation of dimensionless bubble frequency 𝜔 versus 𝛽. Clearly, the analytical approximation for 𝜔 (solid line) is very close to the numerical calculations (symbols). Even though the relative error has a tendency to grow exponentially as shown in Figure 11, its values are acceptably low for large magnitude of the step forcing, for example, for 𝛽=2 (𝑝=1360 Pa) the relative error is just about 1.15%.

4. Conclusions

In this study the HAM has been used to obtain analytical solutions of the Rayleigh equation for the radial oscillations of a MEB in liquid helium. The small-order HAM approximations for freely oscillating bubbles agree very well with numerical solutions even for bubbles with initial radial amplitudes as high as 1.6 times the equilibrium Coulomb radius. The analytical solutions for radius, velocity and acceleration of the freely oscillating bubble wall are accurate enough to accomplish surface stability studies (both parametric and Rayleigh-Taylor instabilities could be computed) with the possibility of both saving calculations and giving a bigger understanding of bubble shape instabilities when compared to solutions from a numerical scheme.

In the case of forced oscillations, the fourth-order HAM solutions for displacement and velocity of the bubble wall agree well with those computed numerically. Nevertheless, when the magnitude of pressure step is large enough (𝛽=1.47832, 𝑝=1000 Pa), noticeable differences are observed between analytical and numerical curves for acceleration. It was shown that higher order approximations can be adapted to increase the convergence of the solution but at the expense of a huge number of terms in it, which might be impractical. Even this inconvenience, current analytical approximations show the great potential of the HAM for complex problems with strong nonlinearities. HAM offers the possibility of controlling in a convenient way the convergence of approximation series; this fundamental characteristic is what makes the HAM more powerful than other nonperturbation techniques such as Adomian decomposition method [15] and the homotopy perturbation method (HPM) [16]. It was confirmed here that convergence of solutions is ensure by choosing parameter in an appropriate way.

Acknowledgments

The authors give thanks to Instituto de Ingeniería-UNAM for providing computing resources. This work was supported in part by DGAPA-PAPIIT-UNAM under Grant IN107509 and by II-FI UNAM under Grant 1135. The authors also wish to thank Julian Espinosa García for his useful comments and suggestions to improve the quality of the paper.