Basins of attraction in a harmonically excited spherical bubble model .

Basins of the periodic attractors of a harmonically excited single spherical gas/vapour bubble were examined numerically. As cavitation occurs in the low pressure level regions in engineering applications, the ambient pressure was set slightly below the vapour pressure. In this case the system is not strictly dissipative and the bubble can grow infinitely for sufficiently high pressure amplitudes and/or starting from large initial bubble radii, consequently, the stable bubble motion is not guaranteed. For moderate excitation pressure amplitudes the exact basins of attraction were determined via the computation of the invariant manifolds of the unstable solutions. At sufficiently large amplitudes transversal intersection of the manifolds can take place, indicating the presence of a Smale horseshoe map and the chaotic behaviour of system. The incidence of this kind of chaotic motion was predicted by the small parameter perturbation method of Melnikov.


Introduction
In many hydraulic systems, in engineering applications cavitation bubbles may form in the low pressure level regions.Propagating toward the high pressure fluid domain, these bubbles can collapse violently causing extensive erosion of the surfaces of the hydraulic components, see e.g.Chan [1], Escaler [2,3].One way of gaining better understanding of the physics of cavitation is to study the dynamics of a single gas/vapour bubble exposed to harmonically varying pressure field.This phenomenon has been extensively studied in the last century, both experimentally and numerically, by applying high amplitude and high frequency sound field to the liquid domain, see the comprehensive review of Lauterborn [4].
The response of this type of harmonically excited bubble is very feature-rich from dynamical point of view.Due to the time varying pressure field, equilibrium points do not exist, consequently, the simplest structure in the system is the periodic solution, whose period T p is an integral multiple of the period of the excitation T 0 , that is, T p = kT 0 , k = 1, 2, • • • .Under parameter variations these orbits can change its stability via bifurcations, such as period doubling (PD) or saddle-node/fold (FL) bifurcations.The simplest way to obtain information about the existing bifurcation curves (BC) is to integrate the system forward in time using a simple initial value problem (IVP) solver, meanwhile, continuously monitoring the possible convergence to a stable periodic/chaotic solution (attractor).One severe drawback of this method is that unstable structures cannot be computed, which are essentially important in the understanding of the dynamics, for instance to locate the basins of attractions.Although, there are already well developed numerical techniques, which are capable to compute the unstable orbits too treating the problem as a boundary value problem (BVP), the majority of the recent papers use the simple IVP method to study the bubble behaviour via bifurcation diagrams, see Akhatov [5], Behnia [6][7][8],Kafiabad [9], hence they often miss the unstable structures.
This study intends to compute the stable periodic orbits and its basins of attraction at constant excitation frequency ω and with varying pressure amplitude p A , employing the simplest and still widely used bubble model the Rayleigh-Plesset (RP) equation, for details see Plesset [10].In order to exploit the benefits of the different solvers, the IVP and BVP methods were combined during the investigation.For gaining a rough global picture about the coexisting attractors of the system, numerous IVP computations were performed in a wide range of the control parameter p A .Initiating the BVP solver (AUTO continuation and bifurcation software, Doedel [11]) from these results, we shall see that complete BCs of periodic orbits can be obtained under parameter variation including the unstable solutions and the detection of bifurcation points as well.
The resulted unstable structures play significant role in the determination of the basins of attraction related to the domains enclosed by the invariant manifolds of the saddle-type solutions.These are particularly important in the present study as our system is not strictly dissipative due to the very low ambient pressure, therefore, the bubble can grow infinitely for sufficiently high pressure amplitudes or for large initial bubble radii, thus, stable bubble motion is not guaranteed.The main aim of this study is to determine basins of attraction of the most significant stable solutions for moderate pressure amplitudes.
An additional application of the invariant manifolds is the determination of the incidence of chaotic motion.Their transversal intersection indicates the presence of a Smale horseshoe map.However, this map is an evidence of chaos it is usually unstable and impossible to compute even with sophisticated numerical algorithms.The parameter value at which the intersection occurs can be predicted by the small parameter perturbation method of Melnikov, see Guckenheimer [12].

Mathematical modelling
The radial motion of the bubble wall is governed by the RP eq.(for detailed description see Plesset [10]) defined as where R = R (t) is the bubble radius at time t, p ∞ (t) is the pressure far from the bubble consisting of a static and a periodic component: where p A is the pressure amplitude and is the angular frequency.The material properties of the liquid and the vapour phase were computed with the Haar-Gallagher-Kell equation of state [13] at temperature T ∞ = 30 • C and pressure P ∞ = 3768Pa (this value will be explained later in this chapter) resulted 995.61kg/m 3 liquid density, 7.973 • 10 −4 N s/m 2 liquid dynamic viscosity, 0.0712N/m surface tension and p V = 4242.7 Pa vapour pressure.In the bubble interior the gas content exhibits in general polytrophic state of change now with exponent n = 1 for simplicity (isothermal behaviour).The initial pressure p g0 and radius R 0 determine the mass of the gas within the bubble: m = 4p g0 R 3 0 /3RT ∞ , where R is the specific gas constant.
In the absence of excitation, the equilibrium radius R E of the bubble can be computed from the following algebraic equation (all time derivatives in the RP eq. are zero) A typical equilibrium radius curve is presented in Fig. 1 for a prescribed mass of gas (given p g0 and R 0 ).It contains a turning point which is usually referred to as Blake critical threshold, see Blake [14].At the critical point the derivative of the tension (p V − P ∞ ) with respect to the equilibrium radius R E is If the initial gas radius R 0 is chosen to be the critical radius R c then, according to equation ( 1), the initial gas pressure has to obey the following relation: p g0 = 2/3nR c .In this case merely the critical radius determines the amount of gas inside the bubble.In our computations the critical radius was R c = 10 −4 m, which is the upper bound of the typical nuclei size (Brennen [15]).
Fig. 1.The equilibrium radius R E curve as a function of the tension p V − P ∞ for a given amount of gas content.The vertical line denotes the applied tension on the system.R s E , R u E and R c are the stable, unstable and critical radii, respectively.The vertical line denote the applied tension to the system.
As this study focuses on the cases when the tension is between zero and the critical value, the tension was set to p V − P ∞ = (p V − P ∞ ) c /2 = 474.7 Pa, marked by the solid vertical line in Fig. 1, from which the static pressure is P ∞ = 3768 Pa.Observe that we have two equilibrium radii; the upper one R u E is unstable while the lower one R s E is stable.In what follows the stable equilibrium radius will be denoted by R E , which, after substituting the previously set numerical values, turns out to be R E = 0.6527R c .Introducing the dimensionless bubble radius y 1 = R/R c , the dimensionless time τ = ωt and defining a dimensionless velocity as y 2 = y 1 , where stands for the derivative with respect to τ , the governing equations can be written as where the parameters are where p re f = ρR 2 c ω 2 is a reference pressure.The linear resonant frequency of the bubble is (Brennen [15]) Due to the explicit time dependence of the system, the problem must be treated in the whole 3 dimensional phase space (y 1 , y 2 , τ).In spite of this difficulty the problem can be restricted to a 2 dimensional iterated map since the phase space is periodic in time with period T 0 , where T 0 = 2π is the period of the excitation, by defining a Poincaré map P : S → S (S ∈ R 2 : (y 1 ; y 2 ) is the Poincaré plane.The dynamics on the Poincaré plane S is clearly seen in Fig. 2 in which a trajectory was computed between the dimensionless time τ = 0 − 3T 0 with initial condition y 0 = y 0 1 ; y 0 2 .Observe that as the phase space is periodic in time, the trajectory was projected back to τ = 0 at the time instant τ = T 0 .Now the Poincaré map can be constructed by sampling the continuous solution at time instants τ = kT 0 , where k ∈ N, thus, the map of an arbitrary point y 0 can be obtained by integrating the system by one period T 0 initiating from the specific point.The end point of the solution is the map of y 0 denoted by P (y 0 ), see Fig. 2 left.If a trajectory starting from y 0 returns exactly at the same point after N iterations (P N (y 0 ) = y 0 ) then the solution is a periodic orbit whose period is exactly T p = NT 0 .For instance, Fig. 2 left and right presents the trajectory of a period 3 solution and its corresponding Poincaré map, respectively.
In order to obtain a global picture about the coexisting attractors, computations were performed with IVP solver at constant ω = 2ω 0 excitation frequency (first subharmonic) by varying the pressure amplitude p A between 10 Pa and 5000 Pa as control parameter.The solver was a standard 5 th order Runge -Kutta scheme with 4 th order embedded error estimation.At each pressure amplitude 10 simulations were performed in order to reveal the coexisting attractors.After the convergence of a solution 64 points were recorded from the Poincaré plane.In Fig. 3 the P (y 1 ) values can be seen as a function of the control parameter p A .The arabic numbers denote the periods of the attractors.Observe, that there are no stable solutions above the pressure amplitude 1700 Pa, although the maximum applied value was Fig. 4. Bifurcation curves of AUTO computations initiated from the solution marked by asterisks in Fig. 3 The black solid and red dashed curves are the stable and unstable orbits, respectively.The black dots denote the fold (saddlenode) bifurcations, while the crosses are period doubling bifurcation points.In case of zero pressure amplitude as a limit case of the uppermost and lowermost bifurcation curves one can obtain original equilibrium radii R s E and R u E .The relevant bifurcation curves are labelled by bc ps , where the value of p is the period of the solution and s is a suitable serial number.
5000 Pa.As we mentioned in the Introduction this is due to the dynamical nature of the system since it is non-strictly dissipative and the bubble can escape from stable domains.It is clear from Fig. 3 that with increasing pressure amplitude the basins of attraction gradually decrease and the chance of finding a stable solution becomes very difficult.This was the reason forcing us to examine the basins, however, only for moderate pressure amplitudes.
As we intend to compute the domains of attraction via the stable and unstable invariant manifolds of the saddle-type orbits, finding these unstable solutions is very important.Such kind of solution cannot be found applying a simple IVP solver even if one integrates the system backward in time.To overcome this problem the AUTO continuation software was employed which is capable of computing whole bifurcation curves, including unstable solutions and bifurcation detection, under parameter variation treating the mathematical problem as a boundary value problem.Because AUTO can handle only autonomous systems, equation (2) has to be extended with two additional decoupled ODEs as follows: marked by asterisks in Fig. 3, a series of complete BCs could be computed including some period doubled curves, see Fig. 4 Moreover, the unstable BC corresponding to R u E was also computed because of its dominant role in global basin of attraction, see the uppermost BC.In Fig. 4 the maximum of the solutions y max 1 is presented instead of the y 1 part of the Poincaré plane due to the special data storage mechanism of AUTO.The black solid and red dashed curves are the stable and unstable solutions, respectively.During the computations the bifurcation points were detected from which the FL bifurcations are denoted by black dots, while the PD bifurcations are marked by crosses.Each BC which is relevant in the next discussion is labelled by bc ps , where the value of p is the period of the solution and s is a suitable serial number.Observe that the FL bifurcation divides a BC into two parts.For instance, both bc 31 and bc 32 belong to the same period 3 family but indicate the upper and the lower branch, respectively.The PD bifurcation has no such effect, see e.g.bc 12 or bc 31 which have stable and unstable parts as well.With the aid of this notation we can refer to any solution by identifying the control parameter, for example, bc 12 (150) means the solution at pressure amplitude 150 Pa.
The above presented AUTO computations provide useful results to obtain the domains of attraction, inasmuch as they are areas enclosed by a pair of stable invariant manifolds of a saddletype unstable solution.By definition the stable (unstable) manifold W s (W u ) consists of points in the Poincaré plane from which the solution tends to a fix point y E as the time goes to infinity (minus infinity), more precisely, W s (y E ) : (S : y → y E , t → ∞), W u (y E ) : (S : y → y E , t → −∞).
In Fig. 5 left the stable manifold W s (black curve) and the unstable manifold W u (red curve) of the saddle-type fix point bc 11 (50) (red cross) can be seen.Due to the rather low p A = 50Pa pressure amplitude, the only other existing structure is the period 1 attractor labelled by bc 12 (50) (black dot).The basin of this stable solution is the light blue area enclosed by the stable manifolds, which is also the global domain of attraction of the system.Any trajectory started out of this region will result in an infinite growth of the bubble.The above mentioned notations in the parentheses hold for all the figures in the present study.
In case of pressure amplitude p A = 150Pa the originally stable bc 12 (50) fix point loses its stability via PD bifurcation resulted in another saddle-type structure marked by bc 12 (150), see Fig. 5 right.The basin of the bifurcated bc 22 (150) period 2 solution is defined by the stable manifolds of the fix point bc 11 (150) corresponding again to the bifurcation curve bc 11 , see Fig. 4 It should be noted that the unstable manifolds of bc 12 (150) have particularly complex shape, moreover, the gradually increasing oscillations of the manifolds belong to bc 11 (150) indicating the forthcoming transversal intersections of these manifolds.
At pressure amplitude p A = 219.6Pa a stable bc 31 and an unstable bc 32 period 3 solution appear via a FL bifurcation, see Fig. 4 In Fig. 6 the basin of this new attractor is represented by the green areas at p A = 250Pa as closed regions of the stable manifolds of its counterpart bc 32 (250).Observe that the three distinct domains contain black dots which are the fix points of the P 3 map and periodically alternate under the influence of the P 1 map.Only one unstable manifold was computed, which tends to the period 2 attractor of bc 22 resulting in a peculiar shape of the curve, see the red curve starting from bc 32 (250) in Fig. 6 The main consequence of the existence of this new structure is that the global basin is the sum of the basin of the period 3 attractor (green area) and the period 2 attractor (light blue area).Further increasing the pressure amplitude, the stable and unstable manifolds of bc 11 will intersect each other.The tangency occurs at approximately p A = 450Pa.One transversal intersection implies an infinite number of intersections resulting in a Poincaré homoclinic structure.The most drastic consequence is the presence of several Smale horseshoe maps, each containing infinite numbers of periodic points with arbitrary high periods.This dynamical behaviour is the evidence of the existence of chaotic motion which is usually unstable explaining the absence of chaotic attractor in Fig. 3 However, after this point the global basin of attraction cannot be computed exactly due to its fractal boundary, it may be possible to approximate via so called basin cells which is beyond the scope of this paper, see the details in (Nusse [16]).
Tab. 1 summarizes the computed basins with respect to the pressure amplitude.Assuming that the global basin will not change significantly from pressure amplitude 450Pa to 650Pa, the second largest stable structure corresponding to bc 31 has the area of basin of approximately 1/20 times the global area of basin.This implies that solution with higher periods have even less domain of attraction.
As we mentioned previously the exact global domain of attraction cannot be computed above p A = 450Pa due to the homoclinic tangency of the stable and unstable manifolds of bc 11 .Because the computation of these invariant manifolds are very resource demanding, it would be very convenient and useful to predict somehow the parameter value at which the intersection takes place.In the following we shall use the method of Melnikov which is based on the perturbation of a planar homoclinic orbit of a Hamiltonian system.If n = 1 (isothermal case) and in the absence of excitation and viscosity equation Eq. (2) reduces  to which now form a Hamiltonian vector field with the Hamiltonian function defined as In this case the system (3) can be written as where Q = 1/6y 3 1 .The stable and unstable manifolds of the unstable equilibrium point R u E of this unexcited time continuous system coincide forming a homoclinic orbit y 0 1 , y 0 2 .Next, the excitation and the viscosity are included as small perturbation: where Finally, we can define the Melnikov function as We should note, that M (τ 0 ) is T 0 = 2π-periodic and it provides a good measure of the separation of the manifolds.Thus, if M (τ 0 ) oscillates about zero, that is, M (τ 0 ) = 0 for some τ 0 and dM(τ 0 )/dτ 0 0, then it follows that there are infinite transversal intersections of the invariant manifolds.Varying the pressure amplitude p A and thus K 3A , we could determine the pressure amplitude at which the intersection occurs by continuously monitoring the values of M (τ 0 ).
However, ε must be sufficiently small, we set it to unity in order to be consistent with the original system (2).Observe, that it is not a real restriction since it is enough that εK 3A and εK 1 be small.The parameter K 1 = 6.1 × 10 −3 was constant, whereas the parameter K 3A was varied between 3.6 × 10 −5 and 1.6 × 10 −2 in the pressure range of 1 Pa and 450 Pa.Although, these parameter values are rather small the Melnikov method provides poor estimation as it approximates the tangency at pressure amplitude p A = 840 Pa.Therefore, in our system, this method is useless for predicting the homoclinic tangency.

Conclusion
In this paper the exact domains of attraction of the stable solutions of a harmonically excited bubble oscillator, the classical Rayleigh-Plesset equation, was examined.As these domains were obtained as an enclosed area of the stable invariant manifolds of the unstable saddle-type solutions, we used the AUTO continuation software to compute bifurcation curves including the unstable solutions.These simulations were initiated from results obtained from a simple initial value problem solver.The ambient pressure was set slightly below the vapour pressure since cavitation occurs at low pressure level regions.In this case the system is not strictly dissipative and the bubble can grow unlimited under certain initial conditions, thus the computation of the basins of attraction is very important.
A comprehensive analysis with a simple initial value problem solver revealed that above the pressure amplitude p A = 1700 Pa the finding of any stable structure becomes very difficult indicating that the area of the global basin is very small.Moreover, the exact domains of attractions showed that the size of the basin of the second largest period 3 structure is less than one-twentieth of the global basin.
The global basin could be computed up to the pressure amplitude p A = 450 Pa.Above this point the stable and unstable invariant manifolds intersect each other and the boundary of the basin becomes fractal.The inception of the tangency was tried to be predicted by the method of Melnikov based on the perturbation of a planar homoclinic orbit.Unfortunately it gave an poor estimation (p A = 840Pa instead of p A = 450Pa).

Fig. 2 .Fig. 3 .
Fig. 2. The construction of the Poincaré map (right) by sampling the continuous solution (left) at time instants τ = kT 0 , where k ∈ N. The present trajectory is a period 3 solution as the Poincaré map returns to itself after 3 iterates, that is, P 3 (y 0 ) = y 0 .

Fig. 5 .
Fig. 5.The global domain of attraction of the system (light blue area) at pressure amplitude p A = 50 Pa (left) and p A = 150 Pa (right).The stable W s and unstable W u invariant manifolds are denoted by black and red curves, respectively.The attractors are marked by black dots, while the unstable saddletype solutions denoted by red crosses.

Fig. 6 .
Fig.6.Basins of the period 3 attractor (green area) and the period 2 attractor (light blue area).The black dots and red crosses are the stable and unstable fix points of a P N map, respectively, where N is the period of the corresponding solution.The global domain of attraction is the sum of these two basins.
Tab. 1.The area of the domain of attractions as a function of the pressure amplitude p A .