# Homotopy Analysis Method for the Rayleigh Equation Governing the Radial Dynamics of a Multielectron Bubble

Authors: F. A. Godínez; M. A. Escobedo; M. Navarrete
Journal: Journal of Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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 [1–3]. 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]: (1.1)R5d2Rdt2+32R4(dRdt)2=-1ρ[2σR3-e2Z232π2εε0],
subject to the initial conditions
(1.2)R(0)=ξRc,dR(0)dt=0,
where R is the bubble radius, ρ and σ are the mass density and surface tension of the liquid Helium, respectively, e is the electron charge, Z is the number of electrons in the bubble, ε is the dielectric constant of helium, ε0 is thevacuum permittivity, Rc 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:
(2.1)N[r(t)]=0,
where N is a nonlinear differential operator, t denotes the independent variable; r(t) is an unknown function.Applying the transformationτ=ωt, (2.1) can be expressed by the following:
(2.2)N[r(τ),ω]=0.
Now a homotopy in general form can be constructed as follows:
(2.3)H-(ϕ,q,h,H(τ))=(1-q)L[ϕ(τ,q)-r0(τ)]-qhH(τ)N[ϕ(τ,q),ω(q)],
where q∈[0,1] is an embedding parameter, ϕ is a function of τ and q, h is a nonzero auxiliary parameter, H(τ) is a nonzero auxiliary function, and L denotes an auxiliary linear operator. When the parameter q increases from 0 to 1, ω(q) varies from ω0 to ω, and the solution ϕ(τ,q) varies from the initial approximation r0(τ) to the exact solution r(τ). Put differently, ϕ(τ,0)=r0(τ) is the solution of H-(ϕ,q,h,H(τ))|q=0=0, and ϕ(τ,1)=r(τ) is the solution of H-(ϕ,q,h,H(τ))|q=1=0. Setting H-(ϕ,q,h,H(τ))|=0, the zero-order deformation equation is obtained as
(2.4)(1-q)L[ϕ(τ,q)-r0(τ)]=qhH(τ)N[ϕ(τ,q),ω(q)],
under the initial conditions:
(2.5)ϕ(0,q)=k,dϕ(0,q)dτ=0,
where k is a constant. The functions ϕ(τ,q) and ω(q) can be expanded as power series of q using Taylor’s theorem as follows
(2.6)ϕ(τ,q)=ϕ(τ,0)+∑m=1∞1m!∂mϕ(τ,q)∂qm|q=0qm=r0(τ)+∑q=1∞rm(τ)qm,ω(q)=ω0+∑m=1∞1m!∂mω(q)∂qm|q=0qm=ω0+∑q=1∞ωmqm,
where rm(τ) and ωm are called the mth-order deformation derivatives.Differentiating the zero-order deformation equation with respect toq and then setting q=0, yields the first-order deformation equation which gives the first-order approximation of r(τ) as follows:
(2.7)L[r1(τ)]=hH(τ)N[r0(τ),ω0],
subject to the initial conditions:
(2.8)r1(0)=0,dr1(0)dτ=0.
Differentiating Equations (2.4) and (2.5) m times with respect to q, then setting q=0 and finally divided them by m!, yields the so called mth-order (m>1) deformation equation:
(2.9)L[rm(τ)-χmrm-1(τ)]=hH(τ)Dm(r→m-1,ω→m-1),
subject to the initial conditions
(2.10)rm(0)=0,drm(0)dτ=0,
where Dm(r→m-1,ω→m-1), r→m-1 and ω→m-1 are defined as follows:
(2.11)Dm(r→m-1,ω→m-1)=1(m-1)!dm-1N[ϕ(τ,q),ω(q)]dqm-1|q=0,(2.12)r→m-1={r0,r1,r2,…,rm-1},(2.13)ω→m-1={ω0,ω1,ω2,…,ωm-1},(2.14)χm={0,m≤1,1,m>1.


## 3. Application of the HAM

### 3.1. Radial Free Oscillations

A nondimensional radiusR- and time t- are defined as follows:
(3.1)R-=R[R],t-=t[t],
where [R] and [t] are characteristic scales.In terms of the dimensionless variables (3.1), the equation of motion (1.1) becomes
(3.2)R-5d2R-dt-2+32R-4(dR-dt-)2=θ-αR-3,
and the initial conditions (1.2) become
(3.3)R-(0)=A,dR-(0)dt-=0,
where α, θ, and A are dimensionless parameters given by the following:
(3.4)α=2σ[t]2ρ[R]3,θ=e2Z2[t]232π2εε0ρ[R]6,A=ξRc[R].
Introducing the new independent variable τ=ωt-, where ω is the dimensionless natural frequency, (3.2) and the initial conditions (3.3) can be rewritten as follows:
(3.5)ω2R-5d2R-dτ2+32ω2R-4(dR-dτ)2=θ-αR-3,(3.6)R-(0)=ξ,dR-(0)dτ=0.
In order to reduce the number of recursions and save computational efforts, this will be clarified later, an additional transformation is introduced
(3.7)R-(τ)=r(τ)3,
then, (3.5) becomes
(3.8)13ω2rd2rdτ2-118ω2(drdτ)2+αr-θ=0,
under the initial conditions
(3.9)r(0)=ξ3,dr(0)dτ=0.
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]:
(3.10){cos⁡(mτ)∣m≥0},
in the compact form
(3.11)r(τ)=∑m=0+∞cmcos⁡(mτ).
With the purpose of satisfying the initial conditions (3.9), the initial guess of r(τ) for zero-order deformation equation is chosen as follows:
(3.12)r0(τ)=(ξ3-1)cos⁡τ+1.
Under the Rule of Solution Expression denoted by (3.11), the linear operator is selected as:
(3.13)L[r(τ;q)]=ω02(∂2r(τ;q)∂r2+r(τ;q)),
with the property
(3.14)L[C1sin⁡τ+C2cos⁡τ]=0,
where C1 and C2 are the integral constants.From (3.8), the nonlinear operator is written as
(3.15)N[r(τ;q),ω(q)]=13ω(q)2r(τ;q)∂2r(τ;q)∂τ2-118ω(q)2(∂r(τ;q)∂τ)2+r(τ;q)-θ.
For complying with the general form of the base functions (3.10), the auxiliary function H(τ) is chosen as follows:
(3.16)H(τ)=1.
Substituting (3.15) into (2.11) and performing the homotopy-derivatives [3], yields:
(3.17)Dm(τ)=13∑i=0m-1rm-1-i′′(τ)∑j=0iri-j(τ)∑l=0jωj-lωl-118∑i=0m-1rm-1-i′(τ)∑j=0iri-j′(τ)∑l=0jωj-lωl+αrm-1(τ)-θ(1-χm),
where a prime indicates first derivative with respect to τ.At this point, it is worth comparing (3.17) with the mth-order homotopy derivative corresponding with (3.5), which reads as follows:
(3.18)Dm(τ)=∑i=0m-1R-m-1-i′′∑j=0iR-i-j∑l=0jR-j-l∑m=0lR-l-m∑p=0mR-m-p∑q=0pR-p-q∑r=0qωq-rωr+32∑i=0m-1R-m-1-i′∑j=0iR-i-j′∑l=0jR-j-l∑m=0lR-l-m∑p=0mR-m-p∑q=0pR-p-q∑r=0qωq-rωr+α∑i=0m-1R-m-1-i∑j=0iR-i-jR-j-θ(1-χm),
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 L, if the term cos(τ) exists in Dm(τ), 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 Dm(τ) must be zero.From (3.17), the first-order approximation (m=1) of HAM for D1(τ) is as follows:
(3.19)D1=(α-θ-736ω02+1436ξ3ω02-736ξ6ω02)+(ξ3α-α+13ω02-13ξ3ω02)cos⁡⁡(τ)+(1036ξ3ω02-536ω02-536ξ6ω02)cos⁡(2τ).
Thus, ω0 results in:
(3.20)ω0=3α.
Solving (2.7) subject to conditions (2.8) r1(τ) is obtained as:
(3.21)r1(τ)=-hsin2(τ/2)[54(θ-α)+ω02(ξ3-1)2(13+5cos⁡(τ))]27ω02.
Similarly, from the coefficient of cos(τ) in R2(τ), ω1 is obtained as follows:
(3.22)ω1=h1296ω03(ξ3-1)(1944αθ-1944α2+1584αω02-1224ξ3αω02+288ξ6αω02-1296θω02+648ξ3θω02-137ω04+315ξ3ω04-219ξ6ω04+41ξ9ω04).
Solving (2.9) under conditions (2.10) for m=2, r2(τ) results in:
(3.23)r2(τ)=-111664ω04×{hsin2(τ2)[-23328hα(α-θ)+ω02(144Γ8+Γ12ω0(Γ7ω0+864Γ6ω1))]},
where
(3.24)Γ1=(ξ3-1),Γ2=(Γ1h(29ξ3-107)-162),Γ3=(27+13Γ1h),Γ4=α(5+ξ3)-6θ,Γ5=(216-7h+103ξ3h),Γ6=(13+5cos⁡(τ)),Γ7=(5616-1379h+1859ξ3h+10Γ5cos⁡(τ)+195Γ1hcos⁡(2τ)),Γ8=(Γ2α+6Γ3θ-5hΓ1Γ4cos⁡(τ)).
The higher-order approximations of rm(τ) and ωm-1(τ), for m>2, are calculated similarly.


### 3.2. Driven Oscillations

The Rayleigh equation (1.1) with a driving pressure step p is written as follows:
(3.25)R5d2Rdt2+32R4(dRdt)2=-1ρ[2σR3+pR4-e2Z232π2εε0],
subject to the initial conditions
(3.26)R(0)=Rc,dR(0)dt=0.
A nondimensional radius x and time t- are defined as follows:
(3.27)x=R[R],t-=t[t].
In terms of the dimensionless variables (3.27), the equation of motion (3.25) becomes
(3.28)x5d2xdt-2+32x4(dxdt-)2=θ-αx3-βx4,
and the initial conditions (3.26) become
(3.29)x(0)=a,dx(0)dt-=0,
where α, β, θ, and a are dimensionless parameters given by the following:
(3.30)α=2σ[t]2ρ[R]3,β=p[t]2ρ[R]2,θ=e2Z2[t]232π2εε0ρ[R]6,a=Rc[R].
Under the transformation τ=ωt-, where ω is the dimensionless natural frequency, (3.28) and the initial conditions (3.29) can be rewritten as follows:
(3.31)ω2x5d2xdτ2+32ω2x4(dxdτ)2+βx4+αx3-θ=0,(3.32)x(0)=a,dx(0)dτ=0.
Equation (3.31) can be integrated once [13], yielding:
(3.33)ω2x4(dxdτ)2+23βx4+αx3-(α+23β+2θ)x+2θ=0.
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., ∑i=0m-1xm-1-i′′∑j=0ixi-j∑ℓ=0jxj-ℓ∑m=0ℓxℓ-m∑p=0mxm-p∑q=0pxp-q∑r=0qωq-rωr, which is at the order of magnitude n8) 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., ∑i=0m-1xm-1-i′∑j=0ixi-j′∑ℓ=0jxj-ℓ∑m=0ℓxℓ-m∑p=0mxm-p∑q=0pxp-q∑r=0qωq-rωr) is at the order of magnitude n8. Even though the term (α+(2/3)β+2θ)x appears in (3.33), its homotopy-derivative (i.e., (α+(2/3)β+2θ)xm-1) is only at the order of magnitude n1.An additional simplification is found when (3.33) is substituted into (3.31), resulting in
(3.34)ω2x5d2xdτ2-12αx3+(32α+β+3θ)x-4θ=0.
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)β(∑i=0m-1xm-1-i∑j=0ixi-j∑ℓ=0jxj-ℓxℓ)) 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:
(3.35){cos⁡(mτ)∣m≥0},
in the form
(3.36)x(τ)=∑m=0+∞cmcos⁡(mτ).
It should be noted that the system (3.34) and (3.32) oscillates between the initial radius x(0)=a and the minimum radius xmin, that is, when bubble wall velocity becomes zero (dx/dτ)=0; then from (3.33):
(3.37)xmin⁡=-13-α2β+(3α-4β)(3α+2β)6βΦ+Φ6β,
where Φ is given by the following:
(3.38)Φ={[8θ(3α+2β)2(7β-3α)-(α-2β)(3α+2β)3+1296β2θ2]1/2324β2θ-(3α-7β)(3α+2β)2+9β[8θ(3α+2β)2(7β-3α)-(α-2β)(3α+2β)3+1296β2θ2]1/2}1/3.
To satisfy the initial conditions (3.32), the initial guess of x(τ) for zero-order deformation equation is chosen as follows:
(3.39)x0(τ)=(a-xmin⁡)2cos⁡τ+(a+xmin⁡)2.
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:
(3.40)N[x(τ;q),ω(q)]=ω2(q)x5(τ;q)∂2x(τ;q)∂τ2-12αx3(τ;q)+(32α+β+3θ)x(τ;q)-4θ.
The solution must comply with the general form of the base functions (3.35). Hence, the auxiliary function H(τ) must be chosen as follows:
(3.41)H(τ)=1.
Substituting (3.40) into (2.11) and performing the homotopy-derivatives, yields:
(3.42)Dm(τ)=∑i=0m-1xm-1-i′′∑j=0ixi-j∑l=0jxj-l∑m=0lxl-m∑p=0mxm-p∑q=0pxp-q∑r=0qωq-rωr-12α∑i=0m-1xm-1-i∑j=0ixi-jxj+(32α+β+3θ)xm-1-4θ(1-χm).
As before, to avoid the secular terms, the coefficient of cos(τ) in Dm(τ) must be zero. Thus, the frequency ω0 is obtained by setting to zero the coefficient of cos(τ) in D1(τ), yielding:
(3.43)ω0=22-(a-xmin⁡)(3(5a2-16)α-32(β+3θ)+3αxmin⁡(6a+5xmin⁡))93a6-68a5xmin⁡-15a4xmin⁡2+15a2xmin⁡4+68axmin⁡5-93xmin⁡6.
Then solving (2.7) under the conditions (2.8), x1(τ) is obtained as follows:
(3.44)x1(τ)=-hsin2(τ/2)215040ω02[(-21αK10-16a3K11ω02sin4(τ2))1680K1+a6K2ω02+5a2K3xmin⁡4ω02-4aK4xmin⁡5ω02+K5xmin⁡6ω02-4xmin⁡×(420K6+a5K7ω02)+5axmin⁡2(336αK8+a3K9ω02)+80xmin⁡3(-21αK10-16a3K11ω02sin4(τ2))],
where
(3.45)K1=-192aα+49a3α-128aβ+1024θ-384aθ+18a3αcos⁡(τ)+a3αcos⁡(2τ),K2=71087+53974cos⁡(τ)+9974cos⁡(2τ)+1674cos⁡(3τ)+199cos⁡(4τ)+12cos⁡(5τ),K3=-2255-310cos⁡(τ)+74cos⁡(2τ)+38cos⁡(3τ)-103cos⁡(4τ)+36cos⁡(5τ),K4=10013+5326cos⁡(τ)-1454cos⁡(2τ)+516cos⁡(3τ)-139cos⁡(4τ)+18cos⁡(5τ),K5=56737+25274cos⁡(τ)-4726cos⁡(2τ)+974cos⁡(3τ)-151cos⁡(4τ)+12cos⁡(5τ),K6=192α-13a2α+128β+384θ+22a2αcos⁡(τ)+3a2αcos⁡(2τ),K7=15963+17226cos⁡(τ)+4846cos⁡(2τ)+1216cos⁡(3τ)+211cos⁡(4τ)+18cos⁡(5τ),K8=19-10cos⁡(τ)+3cos⁡(2τ),K9=-1345+1510cos⁡(τ)+1334cos⁡(2τ)+738cos⁡(3τ)+274cos⁡(4τ)+36cos⁡(5τ),K10=-47-14cos⁡(τ)+cos⁡(2τ),K11=42+49cos⁡(τ)+18cos⁡(2τ)+3cos⁡(3τ).
The higher order approximations for ωm-1 and xm(τ) when m>1 can be similarly calculated.Two sets of results are presented. The first one is for free oscillations (Figures1–7) and the second one is for driven oscillations (Figures 8–11).Figure 1
Theh-curves for the frequency ω. Dash-dotted curve: 2nd-order approximation; solid line: 5th-order approximation (ξ=1.063, α=1, θ=0.999939).Figure 2
Comparison between 2nd-order HAM approximations (solid lines) and Runge-Kutta numerical solutions (open symbols) for:r(τ) (circles), r′(τ) (squares), and r′′(τ) (triangles); ξ=1.063, α=1, θ=0.999939, h=-2.3.Figure 3
The relative error between HAM approximations and numerical solutions for the radiusr(τ) (ξ=1.063, α=1, θ=0.999939, h=-2.3). Dotted curve: 2nd-order approximation; solid curve: 5th-order approximation.Comparison between HAM approximations (hollow symbols) and Runge-Kutta numerical solutions (solid lines) forr(τ) (circles), r′(τ) (squares), and r′′(τ) (stars); (a) 2nd-order approximation, (b) 5th-order approximation. The parameters for both the HAM and numerical solutions were set to: ξ=1.15, α=1, θ=0.999939, h=-2.3.
(a)
(b)Figure 5
Comparison of the 1st, 2nd, and 5th-order HAM approximations ofω with the numerical results (α=1, θ=0.999939, h=-2.3). Symbols: Runge-Kutta numerical results; dotted line: 1st-order; dash-dotted line: 2nd-order; solid line: 5th-order.Figure 6
The relative error for the frequencyω (α=1, θ=0.999939, h=-2.3). Dash-dotted line with symbols: relative error between 2nd-order HAM approximation and numerical results; solid line with symbols: relative error between 5th-order HAM approximation and numerical results.(a) Comparison between 5th-order HAM solution (hollow symbols) and Runge-Kutta numerical solution (solid lines) ofr(τ) for different values of θ with the following fixed values: ξ=1.15, α=1, h=-2.3. Square symbols: Z=10 electrons, Rc=125 Å and θ=0.617405; triangular symbols: for Z=100 electrons, Rc=525 Å and θ=0.83334; diamond symbols: for Z=104 electrons, Rc=1.064 μm and θ=0.999939; circular symbol: for Z=108 electrons, Rc=486 μm and θ=1.05049. (b) Effect of varying the dimensionless initial radius ξ. Comparison between 5th-order HAM solution (hollow symbols) and Runge-Kutta numerical solution (solid lines) for r(τ) with: ξ=1.2 and h=-1.36 (square symbols), ξ=1.3 and h=-1.25 (circular symbols), ξ=1.4 and h=-0.97 (triangular symbols). The remaining parameters are set to: α=1, θ=0.999939.
(a)
(b)Comparison between 4th-order HAM solutions (hollow symbols) and Runge-Kutta numerical solutions (solid lines) for:x(τ), x′(τ) and x′′(τ) at different values of β. (a) α=1, β=0.443495 (p=300 Pa), θ=1, a=1.00002, h=-1.4; (b) α=1, β=1.47832 (p=1000 Pa), θ=1, a=1.00002, h=-1.7. All parameter values are for a MEB containing 104 electrons with an equilibrium radius Rc=1.06441 μm.
(a)
(b)Relative error between 4th-order HAM solution and the corresponding numerical integration for displacementx(τ); (a) and (b) parts of this figure match with those in Figure 8.
(a)
(b)Figure 10
The curve ofω changing with β. Symbols: numerical results; solid line: 4th-order HAM approximation.Figure 11
Relative error between 4th-order HAM approximation and numerical solutions for the frequencyω.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.46)[R]=Rc,[t]=ρRc32σ,
thus, the parameters (3.4) become
(3.47)α=1,θ=e2N264π2εε0σRc3,A=ξ.
Figure 1 shows the curves of ω as function of h for different order approximations. The curves indicate that the valid regions of h (i.e., the interval where ω versus h is a horizontal line) are -3<h<-2 and -3<h<-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 ofh chosen in the valid interval -3<h<-2 is compared with the numerical solution by Runge-Kutta scheme as illustrated in Figure 2.Solutions in Figure2 correspond to a MEB containing 104 electrons, as that numerically simulated by Salomaa and Williams [11], accordingly its equilibrium Coulumb radius Rc is about 1.064 microns. The initial radius of the MEB is 1.063 (precisely the value of ξ) times bigger than the equilibrium radius Rc; 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 Mth-order approximation as follows
(3.48)Relative  Error=|r(τ)N-∑m=0Mrm(τ)|r(τ)N×100,
where r(τ)N is the solution by Runge-Kutta method.The maximum relative error, indicated with a filled triangle in Figure3, 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 Figure4. 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).Figure5 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.Figure7 shows the validity of the 5th-order HAM solutions for r(τ) 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 Figure7(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 h was varied to ensure the convergence of the analytical solution, the freedom of changing conveniently h is precisely one of the virtues of the HAM [1, 6, 9].In order to study the effect of varying the driving pressurep in (3.25), the following scales where chosen:
(3.49)[R]=14e2N2π2εε0σ3,[t]=eN8πσρ2εε0,
then the parameters (3.30), become
(3.50)α=1,β=p8e2N2π2εε0σ43,θ=1,a=4Rcπ2εε0σe2N23.
Figure 8 depicts comparisons between the 4th-order HAM solutions and numerical integrations of: x(τ), x′(τ), and x′′(τ) for two values of the pressure step. Excellent agreement between analytical and numerical solutions is clearly seen for p=300 Pa (β=0.443495); accordingly, the maximum relative error between numerical and analytical results for x(τ) is 0.304%, indicated by a hollow triangle in Figure 9(a). Nevertheless, when the pressure step is increased to p=1000 Pa (β=1.47832), some disagreement between HAM approximations and numerical solutions of: x(τ), x′(τ), and x′′(τ) is observed. In this case the maximum relative error for x(τ) 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 x′′(τ). 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 x(τ) 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 x(τ) 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 x(τ) 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 x(τ), x′(τ) and x′′(τ) would imply more computing time and less accuracy than those obtained by a numerical scheme.Figure10 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 (p=1360 Pa) the relative error is just about 1.15%.


## 3.1. Radial Free Oscillations

A nondimensional radiusR- and time t- are defined as follows:
(3.1)R-=R[R],t-=t[t],
where [R] and [t] are characteristic scales.In terms of the dimensionless variables (3.1), the equation of motion (1.1) becomes
(3.2)R-5d2R-dt-2+32R-4(dR-dt-)2=θ-αR-3,
and the initial conditions (1.2) become
(3.3)R-(0)=A,dR-(0)dt-=0,
where α, θ, and A are dimensionless parameters given by the following:
(3.4)α=2σ[t]2ρ[R]3,θ=e2Z2[t]232π2εε0ρ[R]6,A=ξRc[R].
Introducing the new independent variable τ=ωt-, where ω is the dimensionless natural frequency, (3.2) and the initial conditions (3.3) can be rewritten as follows:
(3.5)ω2R-5d2R-dτ2+32ω2R-4(dR-dτ)2=θ-αR-3,(3.6)R-(0)=ξ,dR-(0)dτ=0.
In order to reduce the number of recursions and save computational efforts, this will be clarified later, an additional transformation is introduced
(3.7)R-(τ)=r(τ)3,
then, (3.5) becomes
(3.8)13ω2rd2rdτ2-118ω2(drdτ)2+αr-θ=0,
under the initial conditions
(3.9)r(0)=ξ3,dr(0)dτ=0.
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]:
(3.10){cos⁡(mτ)∣m≥0},
in the compact form
(3.11)r(τ)=∑m=0+∞cmcos⁡(mτ).
With the purpose of satisfying the initial conditions (3.9), the initial guess of r(τ) for zero-order deformation equation is chosen as follows:
(3.12)r0(τ)=(ξ3-1)cos⁡τ+1.
Under the Rule of Solution Expression denoted by (3.11), the linear operator is selected as:
(3.13)L[r(τ;q)]=ω02(∂2r(τ;q)∂r2+r(τ;q)),
with the property
(3.14)L[C1sin⁡τ+C2cos⁡τ]=0,
where C1 and C2 are the integral constants.From (3.8), the nonlinear operator is written as
(3.15)N[r(τ;q),ω(q)]=13ω(q)2r(τ;q)∂2r(τ;q)∂τ2-118ω(q)2(∂r(τ;q)∂τ)2+r(τ;q)-θ.
For complying with the general form of the base functions (3.10), the auxiliary function H(τ) is chosen as follows:
(3.16)H(τ)=1.
Substituting (3.15) into (2.11) and performing the homotopy-derivatives [3], yields:
(3.17)Dm(τ)=13∑i=0m-1rm-1-i′′(τ)∑j=0iri-j(τ)∑l=0jωj-lωl-118∑i=0m-1rm-1-i′(τ)∑j=0iri-j′(τ)∑l=0jωj-lωl+αrm-1(τ)-θ(1-χm),
where a prime indicates first derivative with respect to τ.At this point, it is worth comparing (3.17) with the mth-order homotopy derivative corresponding with (3.5), which reads as follows:
(3.18)Dm(τ)=∑i=0m-1R-m-1-i′′∑j=0iR-i-j∑l=0jR-j-l∑m=0lR-l-m∑p=0mR-m-p∑q=0pR-p-q∑r=0qωq-rωr+32∑i=0m-1R-m-1-i′∑j=0iR-i-j′∑l=0jR-j-l∑m=0lR-l-m∑p=0mR-m-p∑q=0pR-p-q∑r=0qωq-rωr+α∑i=0m-1R-m-1-i∑j=0iR-i-jR-j-θ(1-χm),
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 L, if the term cos(τ) exists in Dm(τ), 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 Dm(τ) must be zero.From (3.17), the first-order approximation (m=1) of HAM for D1(τ) is as follows:
(3.19)D1=(α-θ-736ω02+1436ξ3ω02-736ξ6ω02)+(ξ3α-α+13ω02-13ξ3ω02)cos⁡⁡(τ)+(1036ξ3ω02-536ω02-536ξ6ω02)cos⁡(2τ).
Thus, ω0 results in:
(3.20)ω0=3α.
Solving (2.7) subject to conditions (2.8) r1(τ) is obtained as:
(3.21)r1(τ)=-hsin2(τ/2)[54(θ-α)+ω02(ξ3-1)2(13+5cos⁡(τ))]27ω02.
Similarly, from the coefficient of cos(τ) in R2(τ), ω1 is obtained as follows:
(3.22)ω1=h1296ω03(ξ3-1)(1944αθ-1944α2+1584αω02-1224ξ3αω02+288ξ6αω02-1296θω02+648ξ3θω02-137ω04+315ξ3ω04-219ξ6ω04+41ξ9ω04).
Solving (2.9) under conditions (2.10) for m=2, r2(τ) results in:
(3.23)r2(τ)=-111664ω04×{hsin2(τ2)[-23328hα(α-θ)+ω02(144Γ8+Γ12ω0(Γ7ω0+864Γ6ω1))]},
where
(3.24)Γ1=(ξ3-1),Γ2=(Γ1h(29ξ3-107)-162),Γ3=(27+13Γ1h),Γ4=α(5+ξ3)-6θ,Γ5=(216-7h+103ξ3h),Γ6=(13+5cos⁡(τ)),Γ7=(5616-1379h+1859ξ3h+10Γ5cos⁡(τ)+195Γ1hcos⁡(2τ)),Γ8=(Γ2α+6Γ3θ-5hΓ1Γ4cos⁡(τ)).
The higher-order approximations of rm(τ) and ωm-1(τ), for m>2, are calculated similarly.


## 3.2. Driven Oscillations

The Rayleigh equation (1.1) with a driving pressure step p is written as follows:
(3.25)R5d2Rdt2+32R4(dRdt)2=-1ρ[2σR3+pR4-e2Z232π2εε0],
subject to the initial conditions
(3.26)R(0)=Rc,dR(0)dt=0.
A nondimensional radius x and time t- are defined as follows:
(3.27)x=R[R],t-=t[t].
In terms of the dimensionless variables (3.27), the equation of motion (3.25) becomes
(3.28)x5d2xdt-2+32x4(dxdt-)2=θ-αx3-βx4,
and the initial conditions (3.26) become
(3.29)x(0)=a,dx(0)dt-=0,
where α, β, θ, and a are dimensionless parameters given by the following:
(3.30)α=2σ[t]2ρ[R]3,β=p[t]2ρ[R]2,θ=e2Z2[t]232π2εε0ρ[R]6,a=Rc[R].
Under the transformation τ=ωt-, where ω is the dimensionless natural frequency, (3.28) and the initial conditions (3.29) can be rewritten as follows:
(3.31)ω2x5d2xdτ2+32ω2x4(dxdτ)2+βx4+αx3-θ=0,(3.32)x(0)=a,dx(0)dτ=0.
Equation (3.31) can be integrated once [13], yielding:
(3.33)ω2x4(dxdτ)2+23βx4+αx3-(α+23β+2θ)x+2θ=0.
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., ∑i=0m-1xm-1-i′′∑j=0ixi-j∑ℓ=0jxj-ℓ∑m=0ℓxℓ-m∑p=0mxm-p∑q=0pxp-q∑r=0qωq-rωr, which is at the order of magnitude n8) 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., ∑i=0m-1xm-1-i′∑j=0ixi-j′∑ℓ=0jxj-ℓ∑m=0ℓxℓ-m∑p=0mxm-p∑q=0pxp-q∑r=0qωq-rωr) is at the order of magnitude n8. Even though the term (α+(2/3)β+2θ)x appears in (3.33), its homotopy-derivative (i.e., (α+(2/3)β+2θ)xm-1) is only at the order of magnitude n1.An additional simplification is found when (3.33) is substituted into (3.31), resulting in
(3.34)ω2x5d2xdτ2-12αx3+(32α+β+3θ)x-4θ=0.
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)β(∑i=0m-1xm-1-i∑j=0ixi-j∑ℓ=0jxj-ℓxℓ)) 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:
(3.35){cos⁡(mτ)∣m≥0},
in the form
(3.36)x(τ)=∑m=0+∞cmcos⁡(mτ).
It should be noted that the system (3.34) and (3.32) oscillates between the initial radius x(0)=a and the minimum radius xmin, that is, when bubble wall velocity becomes zero (dx/dτ)=0; then from (3.33):
(3.37)xmin⁡=-13-α2β+(3α-4β)(3α+2β)6βΦ+Φ6β,
where Φ is given by the following:
(3.38)Φ={[8θ(3α+2β)2(7β-3α)-(α-2β)(3α+2β)3+1296β2θ2]1/2324β2θ-(3α-7β)(3α+2β)2+9β[8θ(3α+2β)2(7β-3α)-(α-2β)(3α+2β)3+1296β2θ2]1/2}1/3.
To satisfy the initial conditions (3.32), the initial guess of x(τ) for zero-order deformation equation is chosen as follows:
(3.39)x0(τ)=(a-xmin⁡)2cos⁡τ+(a+xmin⁡)2.
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:
(3.40)N[x(τ;q),ω(q)]=ω2(q)x5(τ;q)∂2x(τ;q)∂τ2-12αx3(τ;q)+(32α+β+3θ)x(τ;q)-4θ.
The solution must comply with the general form of the base functions (3.35). Hence, the auxiliary function H(τ) must be chosen as follows:
(3.41)H(τ)=1.
Substituting (3.40) into (2.11) and performing the homotopy-derivatives, yields:
(3.42)Dm(τ)=∑i=0m-1xm-1-i′′∑j=0ixi-j∑l=0jxj-l∑m=0lxl-m∑p=0mxm-p∑q=0pxp-q∑r=0qωq-rωr-12α∑i=0m-1xm-1-i∑j=0ixi-jxj+(32α+β+3θ)xm-1-4θ(1-χm).
As before, to avoid the secular terms, the coefficient of cos(τ) in Dm(τ) must be zero. Thus, the frequency ω0 is obtained by setting to zero the coefficient of cos(τ) in D1(τ), yielding:
(3.43)ω0=22-(a-xmin⁡)(3(5a2-16)α-32(β+3θ)+3αxmin⁡(6a+5xmin⁡))93a6-68a5xmin⁡-15a4xmin⁡2+15a2xmin⁡4+68axmin⁡5-93xmin⁡6.
Then solving (2.7) under the conditions (2.8), x1(τ) is obtained as follows:
(3.44)x1(τ)=-hsin2(τ/2)215040ω02[(-21αK10-16a3K11ω02sin4(τ2))1680K1+a6K2ω02+5a2K3xmin⁡4ω02-4aK4xmin⁡5ω02+K5xmin⁡6ω02-4xmin⁡×(420K6+a5K7ω02)+5axmin⁡2(336αK8+a3K9ω02)+80xmin⁡3(-21αK10-16a3K11ω02sin4(τ2))],
where
(3.45)K1=-192aα+49a3α-128aβ+1024θ-384aθ+18a3αcos⁡(τ)+a3αcos⁡(2τ),K2=71087+53974cos⁡(τ)+9974cos⁡(2τ)+1674cos⁡(3τ)+199cos⁡(4τ)+12cos⁡(5τ),K3=-2255-310cos⁡(τ)+74cos⁡(2τ)+38cos⁡(3τ)-103cos⁡(4τ)+36cos⁡(5τ),K4=10013+5326cos⁡(τ)-1454cos⁡(2τ)+516cos⁡(3τ)-139cos⁡(4τ)+18cos⁡(5τ),K5=56737+25274cos⁡(τ)-4726cos⁡(2τ)+974cos⁡(3τ)-151cos⁡(4τ)+12cos⁡(5τ),K6=192α-13a2α+128β+384θ+22a2αcos⁡(τ)+3a2αcos⁡(2τ),K7=15963+17226cos⁡(τ)+4846cos⁡(2τ)+1216cos⁡(3τ)+211cos⁡(4τ)+18cos⁡(5τ),K8=19-10cos⁡(τ)+3cos⁡(2τ),K9=-1345+1510cos⁡(τ)+1334cos⁡(2τ)+738cos⁡(3τ)+274cos⁡(4τ)+36cos⁡(5τ),K10=-47-14cos⁡(τ)+cos⁡(2τ),K11=42+49cos⁡(τ)+18cos⁡(2τ)+3cos⁡(3τ).
The higher order approximations for ωm-1 and xm(τ) when m>1 can be similarly calculated.Two sets of results are presented. The first one is for free oscillations (Figures1–7) and the second one is for driven oscillations (Figures 8–11).Figure 1
Theh-curves for the frequency ω. Dash-dotted curve: 2nd-order approximation; solid line: 5th-order approximation (ξ=1.063, α=1, θ=0.999939).Figure 2
Comparison between 2nd-order HAM approximations (solid lines) and Runge-Kutta numerical solutions (open symbols) for:r(τ) (circles), r′(τ) (squares), and r′′(τ) (triangles); ξ=1.063, α=1, θ=0.999939, h=-2.3.Figure 3
The relative error between HAM approximations and numerical solutions for the radiusr(τ) (ξ=1.063, α=1, θ=0.999939, h=-2.3). Dotted curve: 2nd-order approximation; solid curve: 5th-order approximation.Comparison between HAM approximations (hollow symbols) and Runge-Kutta numerical solutions (solid lines) forr(τ) (circles), r′(τ) (squares), and r′′(τ) (stars); (a) 2nd-order approximation, (b) 5th-order approximation. The parameters for both the HAM and numerical solutions were set to: ξ=1.15, α=1, θ=0.999939, h=-2.3.
(a)
(b)Figure 5
Comparison of the 1st, 2nd, and 5th-order HAM approximations ofω with the numerical results (α=1, θ=0.999939, h=-2.3). Symbols: Runge-Kutta numerical results; dotted line: 1st-order; dash-dotted line: 2nd-order; solid line: 5th-order.Figure 6
The relative error for the frequencyω (α=1, θ=0.999939, h=-2.3). Dash-dotted line with symbols: relative error between 2nd-order HAM approximation and numerical results; solid line with symbols: relative error between 5th-order HAM approximation and numerical results.(a) Comparison between 5th-order HAM solution (hollow symbols) and Runge-Kutta numerical solution (solid lines) ofr(τ) for different values of θ with the following fixed values: ξ=1.15, α=1, h=-2.3. Square symbols: Z=10 electrons, Rc=125 Å and θ=0.617405; triangular symbols: for Z=100 electrons, Rc=525 Å and θ=0.83334; diamond symbols: for Z=104 electrons, Rc=1.064 μm and θ=0.999939; circular symbol: for Z=108 electrons, Rc=486 μm and θ=1.05049. (b) Effect of varying the dimensionless initial radius ξ. Comparison between 5th-order HAM solution (hollow symbols) and Runge-Kutta numerical solution (solid lines) for r(τ) with: ξ=1.2 and h=-1.36 (square symbols), ξ=1.3 and h=-1.25 (circular symbols), ξ=1.4 and h=-0.97 (triangular symbols). The remaining parameters are set to: α=1, θ=0.999939.
(a)
(b)Comparison between 4th-order HAM solutions (hollow symbols) and Runge-Kutta numerical solutions (solid lines) for:x(τ), x′(τ) and x′′(τ) at different values of β. (a) α=1, β=0.443495 (p=300 Pa), θ=1, a=1.00002, h=-1.4; (b) α=1, β=1.47832 (p=1000 Pa), θ=1, a=1.00002, h=-1.7. All parameter values are for a MEB containing 104 electrons with an equilibrium radius Rc=1.06441 μm.
(a)
(b)Relative error between 4th-order HAM solution and the corresponding numerical integration for displacementx(τ); (a) and (b) parts of this figure match with those in Figure 8.
(a)
(b)Figure 10
The curve ofω changing with β. Symbols: numerical results; solid line: 4th-order HAM approximation.Figure 11
Relative error between 4th-order HAM approximation and numerical solutions for the frequencyω.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.46)[R]=Rc,[t]=ρRc32σ,
thus, the parameters (3.4) become
(3.47)α=1,θ=e2N264π2εε0σRc3,A=ξ.
Figure 1 shows the curves of ω as function of h for different order approximations. The curves indicate that the valid regions of h (i.e., the interval where ω versus h is a horizontal line) are -3<h<-2 and -3<h<-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 ofh chosen in the valid interval -3<h<-2 is compared with the numerical solution by Runge-Kutta scheme as illustrated in Figure 2.Solutions in Figure2 correspond to a MEB containing 104 electrons, as that numerically simulated by Salomaa and Williams [11], accordingly its equilibrium Coulumb radius Rc is about 1.064 microns. The initial radius of the MEB is 1.063 (precisely the value of ξ) times bigger than the equilibrium radius Rc; 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 Mth-order approximation as follows
(3.48)Relative  Error=|r(τ)N-∑m=0Mrm(τ)|r(τ)N×100,
where r(τ)N is the solution by Runge-Kutta method.The maximum relative error, indicated with a filled triangle in Figure3, 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 Figure4. 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).Figure5 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.Figure7 shows the validity of the 5th-order HAM solutions for r(τ) 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 Figure7(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 h was varied to ensure the convergence of the analytical solution, the freedom of changing conveniently h is precisely one of the virtues of the HAM [1, 6, 9].In order to study the effect of varying the driving pressurep in (3.25), the following scales where chosen:
(3.49)[R]=14e2N2π2εε0σ3,[t]=eN8πσρ2εε0,
then the parameters (3.30), become
(3.50)α=1,β=p8e2N2π2εε0σ43,θ=1,a=4Rcπ2εε0σe2N23.
Figure 8 depicts comparisons between the 4th-order HAM solutions and numerical integrations of: x(τ), x′(τ), and x′′(τ) for two values of the pressure step. Excellent agreement between analytical and numerical solutions is clearly seen for p=300 Pa (β=0.443495); accordingly, the maximum relative error between numerical and analytical results for x(τ) is 0.304%, indicated by a hollow triangle in Figure 9(a). Nevertheless, when the pressure step is increased to p=1000 Pa (β=1.47832), some disagreement between HAM approximations and numerical solutions of: x(τ), x′(τ), and x′′(τ) is observed. In this case the maximum relative error for x(τ) 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 x′′(τ). 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 x(τ) 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 x(τ) 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 x(τ) 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 x(τ), x′(τ) and x′′(τ) would imply more computing time and less accuracy than those obtained by a numerical scheme.Figure10 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 (p=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, p=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 h parameter in an appropriate way.



---

*Source: 591058-2012-09-13.xml*