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 [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]: subject to the initial conditions 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, 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: 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: Now a homotopy in general form can be constructed as follows: where 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 to , and the solution varies from the initial approximation to the exact solution . Put differently, is the solution of , and is the solution of . Setting , the zero-order deformation equation is obtained as under the initial conditions: where is a constant. The functions and can be expanded as power series of using Taylor’s theorem as follows where and are called the mth-order deformation derivatives.
Differentiating the zero-order deformation equation with respect to and then setting , yields the first-order deformation equation which gives the first-order approximation of as follows: subject to the initial conditions: Differentiating Equations (2.4) and (2.5) times with respect to , then setting and finally divided them by m!, yields the so called th-order () deformation equation: subject to the initial conditions where , and are defined as follows:
3. Application of the HAM
3.1. Radial Free Oscillations
A nondimensional radius and time are defined as follows: where and are characteristic scales.
In terms of the dimensionless variables (3.1), the equation of motion (1.1) becomes and the initial conditions (1.2) become where α, , and are dimensionless parameters given by the following: Introducing the new independent variable , where is the dimensionless natural frequency, (3.2) and the initial conditions (3.3) can be rewritten as follows: In order to reduce the number of recursions and save computational efforts, this will be clarified later, an additional transformation is introduced then, (3.5) becomes under the initial conditions 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]: in the compact form With the purpose of satisfying the initial conditions (3.9), the initial guess of for zero-order deformation equation is chosen as follows: Under the Rule of Solution Expression denoted by (3.11), the linear operator is selected as: with the property where and are the integral constants.
From (3.8), the nonlinear operator is written as For complying with the general form of the base functions (3.10), the auxiliary function is chosen as follows: Substituting (3.15) into (2.11) and performing the homotopy-derivatives [3], yields: 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: 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 exists in , the so-called secular term will appear in the final solution disobeying the rule solution expression (3.11). To avoid the secular terms, the coefficient of in must be zero.
From (3.17), the first-order approximation () of HAM for is as follows: Thus, results in: Solving (2.7) subject to conditions (2.8) is obtained as: Similarly, from the coefficient of in , is obtained as follows: Solving (2.9) under conditions (2.10) for , results in: where The higher-order approximations of and , for , are calculated similarly.
3.2. Driven Oscillations
The Rayleigh equation (1.1) with a driving pressure step is written as follows: subject to the initial conditions A nondimensional radius and time are defined as follows: In terms of the dimensionless variables (3.27), the equation of motion (3.25) becomes and the initial conditions (3.26) become where , , , and are dimensionless parameters given by the following: Under the transformation , where is the dimensionless natural frequency, (3.28) and the initial conditions (3.29) can be rewritten as follows: Equation (3.31) can be integrated once [13], yielding: 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., , which is at the order of magnitude ) 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., ) is at the order of magnitude . Even though the term appears in (3.33), its homotopy-derivative (i.e., ) is only at the order of magnitude .
An additional simplification is found when (3.33) is substituted into (3.31), resulting in 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., ) 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: in the form It should be noted that the system (3.34) and (3.32) oscillates between the initial radius and the minimum radius , that is, when bubble wall velocity becomes zero ; then from (3.33): where is given by the following: To satisfy the initial conditions (3.32), the initial guess of for zero-order deformation equation is chosen as follows: 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: The solution must comply with the general form of the base functions (3.35). Hence, the auxiliary function must be chosen as follows: Substituting (3.40) into (2.11) and performing the homotopy-derivatives, yields: As before, to avoid the secular terms, the coefficient of in must be zero. Thus, the frequency is obtained by setting to zero the coefficient of in , yielding: Then solving (2.7) under the conditions (2.8), is obtained as follows: where The higher order approximations for and when can be similarly calculated.
Two sets of results are presented. The first one is for free oscillations (Figures 1–7) and the second one is for driven oscillations (Figures 8–11).
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
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 ): thus, the parameters (3.4) become 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 and 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 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 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 . 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 . 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 , whereas a big bubble containing 108 electrons attains at the first main collapse a minimum dimensionless radius ≈0.72 and develops a dimensionless frequency . 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 , 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: then the parameters (3.30), become 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 Pa (); 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 Pa (), 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 ( 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 (, 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.