Pattern Identification on the Nonlinear Radial Motion of an Oscillating Spherical Bubble Using Neural Networks

Abstract

The main goal of this article is to study the oscillatory motion of a spherical gas bubble immersed in a Newtonian liquid subjected to a harmonic pressure excitation. We use the classical Rayleigh-Plesset equation to study the radial motion of the bubble undergoing a forcing acoustic pressure field. The second order nonlinear ordinary differential equation that governs the bubble motion is solved through a robust fifth order Runge-Kutta scheme with adaptive time-step. Several interesting patterns are identified. First we develop an asymptotic solution for low amplitudes of excitation pressure to validate our numerical code. Then we develop a bifurcation diagram in order to show how the parameters of the flow modify the vibrational patterns of the bubble. We also train a neural network to identify the vibrational pattern through its FFT data. The combination of neural networks with a bifurcation diagram could be useful for the identification of the flow physical parameters in practical applications. For each pattern we also provide an analysis of the motion of the bubble on the phase-space and interpret physically the system behavior with its FFT. In addition, we analyze nonlinear patterns using standard tools of dynamical systems such as Poincaré sections and calculate the Lyapunov exponents of the system. Based on that, we have identified topological transitions in phase plane using for instance the analysis of Poincaré sections and the solution in the frequency spectrum. We have seen that the mechanisms that dominate the dynamics of the oscillating bubble is the competition of the acoustic field excitation with surface tension forces and momentum diffusion by the action of the surrounding fluid viscosity.

Keywords:
Bubble dynamics; nonlinear dynamics; neural networks; acoustic forcing; phase space; motion transition; spherical bubble; radial motion

1 INTRODUCTION

The classical equation describing the oscillatory motion of a spherical bubble subjected to a harmonic acoustic excitation field is known as Rayleigh-Plesset equation (Rayleight, 1917). This highly nonlinear ODE can be modified or extended to study several phenomena in the context of bubble dynamics. For decades, engineers and scientists have productively employed single-degree-of-freedom models to understand and explore the behavior of violently collapsing bubbles (Leighton, 1994; Young, 1989; Brenner et al., 2002). These studies have been applied in many areas such as: sonochemistry, bioengineering, cavitation erosion, underwater explosions and sonoluminescence (Geers et al., 2012).

The shape of the bubble can vary from spherical to nearly polyhedral, forming a complex geometrical structure insensitive to details of the liquid composition or the average bubble size (Gopal et al., 1995). However, when a small bubble is subjected to a small oscillating amplitude, we can assume that its shape remains spherical. This assumption is valid even in highly nonlinear situations. A bubble immersed in a non-Newtonian fluid, such as a viscoelastic polymeric solution (Street, 1968; Fogler, 1969; Albernaz and Cunha, 2011; Albernaz and Cunha, 2013) and even a magnetic fluid (Cunha et al. 2002; Malvar et al. 2016) may have several degrees of freedom, presenting very nonlinear oscillations.

This nonlinearity is related to a more complex control of the bubble's oscillation. When used in biomedical applications, such as contrast in ultrasonography (Chang et al., 1999), it is necessary to understand how the bubble is oscillating. In this case microbubbles are injected into the body and take part in image contrast. The ultrasound pulses are applied with the bubbles' resonance frequency. The bubbles respond to pulses increasing and decreasing their radius, generating echoes in their neighborhoods. The difference between these echoes and those shown by tissues generates a contrast. Furthermore, bubbles may be used to transport drugs (Marmottant, 2004), another example in which the control of its degrees of freedom may be very important.

Considering the relevance of understanding a bubble´s oscillating patterns and identifying them in order to control its motion, a nonlinear tool such as neural networks can be used. In practical applications, the identification of these parameters can become somewhat difficult due to the nonlinearity of the equation that controls the oscillatory regime of the bubble. However, it is a simple problem of pattern identification. Recent works (Sor, 2012; Schem, 2012, Behnia et al., 2013) have analyzed the nonlinear motion of a spherical bubble from a nonlinear dynamical system perspective, but none of them have used neural networks in order to identify some of its nonlinear patterns of motion.

The physical and mathematical theory of neural networks have been developed rapidly during the past 25 years. It is a theory whose diversity and complexity reflects the multifaceted organization of the brain in processes that it sets out to explain and identify (Grossberg, 1988). A multilayer neural network can be used in order to identify defects on metal beams and other dynamic systems (Genovese, 2001). The vibrational parameters of the beam and its response in the frequency domain are used to determine the involvement thereof.

In the last 25 years, multilayer perceptrons (MLPs) have been massively used in the area of pattern recognition (Gori, 1998). The experimental results have been impressive in some applications where we know in advance that the patterns belong to a small number of classes. In those cases, because of their strong discrimination capabilities, MLPs exhibit excellent performance (Waibel, 1989; Lecun 1989).

Considering the microbubble dynamic parameters, such as the system non-dimensional numbers, several vibrational patterns may be identified. Those patterns are classified and recognized by the neural network in order to predict the bubble's oscillatory motion considering its known parameters and vice versa.

It is important to highlight that recent works have explored the unsteady motion of an oscillating bubble immersed in complex fluids. For instance, Cunha et al. (2002) explored the motion of a spherical bubble immersed in a magnetic fluid. They have used a constant applied magnetic field and found out that magnetic effects may lead to a stabilization of the bubble´s motion. In other recent works, Albernaz and Cunha (2013) and Cunha and Albernaz (2013) have explored the radial motion of an oscillating bubble in a viscoelastic fluid. They investigated the effects of the fluid's elasticity on the radial motion of the bubble. In a recent work (submitted for publication) Malvar et al. (2016) have improved the model used by Cunha et al. (2002) in order to include spatial variations of the applied field in the framework of an oscillating bubble immersed in a magnetic fluid. It is important to highlight that in all of these works the physics was interpreted purely by phase space diagrams and time signal response analysis. In none of them the authors have used the set of tools proposed in the present work, to explore the complexity of the nonlinear time series that describes the radial motion of the bubble. In this sense, the present work may serve as an inspiration, in terms of a methodological approach, to future works exploring problems in the framework of bubble dynamics in complex fluids, since we use a larger set of tools to physically explain the behavior of our complex system.

In the present study, we use a sinusoidal excitation with variable amplitude and frequency as the pressure forcing field. The main goal is to identify the non-dimensional parameters and excitations that generate different vibrational modes, characterizing the system. Those vibrational patterns are used as a training input in a multilayer backpropagation neural network in order to identify them from the flow system non-dimensional physical parameters. This is a first step towards the understanding of the nonlinear response of the bubble from a dynamical system perspective.

From our knowledge this work seems to be a first step towards understanding and controlling the nonlinear response of bubbles in oscillatory motions from a dynamical system point of view. We identify topological transitions in phase plane using for instance the analysis of Poincaré sections and the solution in the Fourier space.

When identifying the history of bubble patterns and oscillations, the present analysis can be used to predict the collapse time and also to control the set of dynamic parameters that could avoid it. The presence of microbubbles enormously enhances delivery of genetic material, proteins and smaller chemical agents. Delivery of genetic material is greatly enhanced by ultrasound in the presence of microbubbles (Pitts, 2004). In this particular case, after characterizing the vibrational pattern of the bubble nonlinear oscillations, new investigations and studies could be performed in order to extract the full details of the bubble dynamic´s response. Moreover, determining which vibrational pattern one desires, the ultrasound may be adjusted likewise. In this context, using neural networks to identify cavitation vibrational modes seems to be a quite effective process.

2 MATHEMATICAL MODELLING

2.1 Rayleigh-Plesset Equation

The study of the dynamic oscillatory movement of a bubble presented in this work is based on the analysis of its radial motion when immersed in an incompressible fluid subjected to a harmonic excitation pressure field. The main governing equations used in the mathematical modeling of the problem are based on the principles of mass conservation and linear momentum balance. An appropriate constitutive model to represent the traction jump (i.e. boundary condition for the stresses) on the bubble surface is also used. The velocities are continuous on the surface.

Consider a spherical bubble immersed in a Newtonian incompressible fluid of viscosity μ, and density ρ. The inner side of the bubble is composed by a mixture of air and contaminant gases (which develop a polytrophic process) and steam. We assume that the bubble develops only radial motion due to its surface tension that resists to other non-radial deformational modes, simplifying our analysis to an unidimensional motion, maintaining its spherical shape. This assumption is valid for small excitation amplitudes or for small bubbles. Therefore, effects like pressure gradients in the liquid or surfactants in the fluid are neglected. Changes in the vapor temperature may modify its density modifying the bubble dynamics (Sog, 1996). However, for a small equilibrium vapor density the isothermal process is valid (Brujan, 1999).

The mass conservation equation for a compressible fluid is given by Batchelor (1967)

(1)

in which ρ is the fluid density, t represents time and u is the Eulerian velocity field on the liquid side. In this model spherical coordinates are used. The velocity components in the directions θ and φ are zero. Now, the mass conservation equation in spherical coordinates for the case of an incompressible fluid reduces to

(2)

where r represents the physical distance of the center of the bubble to an arbitrary point inside the liquid phase, u_r is the radial component of the velocity field. Integrating equation (2) from the surface of the bubble of radius R(t) to an arbitrary point in the liquid (distance r), after some algebraic manipulation, we have

(3)

The equation that expresses the balance of linear momentum is given by Batchelor (1967)

(4)

where Σ = -pI + 2μD is the stress tensor for a Newtonian incompressible fluid, written in terms of the pressure field p, the identity tensor I and the rate of strain tensor D = 1/2[∇u + (∇u)^T ]. Considering the radial component of Navier-Stokes equation in spherical coordinates, equation (4) reduces to

(5)

Integrating the resulting equation in the flow field (the region between the surface of the bubble and the infinite quiescent liquid) after some algebraic manipulation, we have

(6)

Where p _∞ is the ambient pressure and p _ℓ denotes the pressure on the liquid phase at the bubble surface. As the bubble is considered in the present context a clean curved interface between two immiscible fluids, there is a normal stress discontinuity in the interface due to its surface tension. This jump of normal stresses is given by Young-Laplace equation:

(7)

where _rr,ℓ and Σ_rr,b denote the radial component of the traction (σ) at the liquid phase and inside the bubble phase, respectively. The pressure inside the bubble is given by p_b (t) = p_v + p_g (t) where p_v is the liquid vapor pressure given by the Classius-Clapeyron relation and p_g (t) refers to the pressure of the gas, assumed as perfect and subjected to a polytropic process: p_{g ^Vn} = constant, where n is the polytropic coefficient, that will be considered unitary in this work, representing an isothermal process. The value of the polytropic index provides important information regarding the thermodynamic processes. Considering an equilibrium pressure that corresponds to the bubble condition when its radius is R_E (an equilibrium radius), we may assume that the internal pressure is given by represents the surface tension coefficient, Σ

(8)

Assuming that in the equilibrium condition the bubble radius is R_E , the pressure discontinuity in the interface assumes the form given bellow, base don the Young-Laplace condition:

(9)

where represents the internal pressure of equilibrium and represents the static pressure of equilibrium applied by the fluid in the initial time. Substituting equation (9) in , one reaches:

(10)

Thus, returning the term _rr,b = - p_b (t) we have to its orinal form, and defining Σ

(11)

The fluid velocity field induced by the bubble oscillations are given by equation (3). Hence, Σ_rr,ℓ results in

(12)

Substituting equations (11) and (12) in (7) we obtain

(13)

where p _∞ (t) = . The forcing acoustic pressure field imposed on the bubble is taken as being ∞ [1 + ε sin(ωt)], where ε represents the pressure amplitude and ω is the forcing pressure frequency. Now, substituting equation (13) in (6) the Rayleigh-Plesset equation is written as

(14)

2.2 Non-Dimensional Formulation

We shall make use of non-dimensional variables indicated by asterisks, namely

(15)

Typical scales of velocity and time are represented by

(16)

respectively. The excitation amplitude may also be written as ε* = ε(). Using these typical scales we can write the non-dimensional form of Rayleigh-Plesset equation (14) as being:

(17)

here, Re = ρU_cR_E /μ corresponds to the Reynolds number which measures the relative importance between inertial and viscous forces and We = ρU_c ² R_E / is the Weber number, which denotes the ratio of inertial and surface tension forces. From now on, all the variables will be presented in their nondimensional form and the upper asterisks will be suppressed in order to use a lighter and clear mathematical notation.

Equation (14) represents a classical equation in the framework of bubble dynamics. In practical applications involving cavitating bubbles the identification of these parameters and their dynamic control in nonlinear bubble oscillatory motion is a quite important feature for understanding and controlling the instability and cavitation in such a nonlinear system.

2.3 Neural Networks

As mentioned before, the main goal of this work is to show that it is possible to conceive, design and train a Neural Network used for pattern identification purposes in the specific dynamical problem of an oscillating micro bubble immersed in a Newtonian liquid. Multilayer neural networks are usually composed of an input layer, a hidden layer and an output layer. Those layers can be formed by activation functions and different training methods. Based on the network's input and random initial conditions, a transfer function is then used to calculate the output layer (Kang, 1996). The output y_k of the kth neuron of the Network is given by

(18)

where y_k is the output of the neural network, w_kj determines the weights given for each entry, x_j is the number of inputs and Ø is the transfer function used. In the present work the output of the network y_k is the number of the pattern found by the network, the input x_j is the information used in order to train the network. In the present work this information is composed by the amplitude and part of the frequency spectrum around the first harmonic. A sketch of the scheme used to train the network is shown in figure (1).

There are several different training methods available in the current literature (e.g. Sanaga 2006). In the present work the Levenberg-Marquardt (LM) backpropagation method (Chlouverakis, 2005) is explored. The Levenberg-Marquardt method seeks to minimize a nonlinear function. In this case the error in the pattern identification process is defined as

(19)

where ∊_φ,o is the error in identifying a given pattern φ for a given output o of the network, d _φ,o is the desired value of the pattern φ that the network should identify in the output o and n_φ,o is the current output that the network provides for a known pattern φ in a current output o.

Figure 1:
Schematic of a typical Neural Network explored in the present work.

The idea of the Levenberg-Marquardt method is to minimize this error ∊_φ,o for different patterns φ in different outputs o of the network. In this sense the network must adjust its weights w_kj for each input and each known pattern during the training period in order to maximize its chances of predicting the correct pattern. Basically the Levenberg-Marquardt algorithm uses the Jacobian matrix containing the first derivatives of the Network errors with respect to its weights and biases, given by

(20)

where n is the number of weights used for each neuron of the network. The minimization of the error in the identification of each pattern is done through the following equation

(21)

where α is called the combination coefficient. Depending on the value of α the LM algorithm switches between the steepest descent and the Gauss-Newton algorithms during the training process. The use of the Jacobian matrix is simpler to compute in comparison with the Hessian matrix in minimization methods such as the Newton algorithm. In this direction, the Levenberg-Marquardt algorithm is characterized by a fast convergence rate in comparison to other algorithms used to minimize nonlinear functions.

The nonlinearity of the bubble's oscillatory motion may generate not only pure harmonics, but also some spectral spreading. This phenomenon can be identified in the frequency spectrum of the signal R(t). The amplitude of the first harmonic and its spreading can be used to diagnose the nonlinearity degree and the involved vibrational modes in the coupled motion of the bubble and the surrounding liquid. We will show that this dynamical system has spectra with broad bands, similar to chaotic systems rather than isolated peaks. In this context a training set of 23 different simulations containing 7 identified patterns is inserted into the neural network. The validation group also possesses the same amount of information, but in this case, the neural network should identify the pattern from the frequency signature. Seven tests are conducted with different random initial conditions of the neural network weights. Each test improves the weights and the identification becomes more effective, displaying a very good percentage of correct answers in the sample.

The training of a neural network is a slow process, in particular because it requires the use of different random weights in its initial conditions. On the other hand, its application is rapid and simple once it is trained in order to perform a specific task. The current literature recommends that the number of epochs, neurons and training algorithm is an important factor for the choice of the network weights (Sanaga, 2006). More details on the implementation, training and performance of the Neural Network will be given in subsequent sections.

3 NUMERICAL METHODOLOGY

Due to the nonlinearity of the governing differential equation we do not have an analytical solution for arbitrary input parameters. In this sense we must recur to the use of a numerical algorithm to integrate the governing equation of the bubble-liquid dynamical system. In this work we use a Runge-Kutta scheme (RKS). Our RKS algorithm uses the advanced predictor-corrector algorithm, which uses weighted averages of a function f calculated at the extremes and at intermediate points of several time intervals [t_j, t_j+1 ] (Malvar, Gontijo and Cunha, 2016).

The nonlinearity of our dynamical system requires a careful calibration of the numerical time-step. For the most unstable regions where the radius of the bubble presents a rapid variation in time, the time step must be very small. On the other hand, when the bubble's oscillatory motion is harmonic, the computational cost may be unnecessarily increased by using small time steps. In order to solve this numerical issue, an adaptive time scheme is incorporated to the dynamic simulations. Here we use the adaptive time step method scheme first proposed by Fehlberg (Press et al., 1992) that uses the difference between two R predictions. The first prediction of R using a fifth order scheme and the second prediction of R by means of a fourth order scheme. In this way we have a good prediction of the truncation error by dynamically adjusting the numerical time step. The formula written in terms of k_n constants for each time step is given by Cash et al. (1990):

(22)

The recurrence formula of a fourth order RKS is given by:

(23)

with the associated error given by

(24)

In bubble dynamics studies it is useful to establish a minimum radius R_min , associated with a collapse criterion, to stop the simulation. In order to evaluate the minimum radius at the eminence of the collapse, an asymptotic theory (Albernaz and Cunha, 2013) is used as criteria to stop our simulations. In this asymptotic theory the most adverse flow condition is considered, that is Re → ∞, and a constant ambient pressure p _∞(t) = R_min results in that compresses the bubble. Under this condition we can solve the Rayleigh-Plesset equation using the integrating factor method. The asymptotic expression for

(25)

where R_min is the minimum radius, R_o is the initial radius and We is the Weber number.

It should be noted that the collapse in the present context is defined as an abrupt release of pressure that occurs outside the spectrum of continuum mechanics time scales. This sudden expansion occurs after significant rapid shape deformations. During these extremely small time intervals the bubble cannot be considered as an isotropic sphere. It is important to say that we are not focused in exploring bubble collapse in the present work. Instead, we are interested to explore the vibrational patterns associated with the nonlinear motion of a spherical bubble and to identifying them through the use of Neural Networks.

3.1 Numerical Tools

In this work we have used a series of numerical tools. The code that computes the motion of the bubble was written in Fortran95 for Linux platforms. For the post-processing of the numerical results, the Softwares Matlab and Tecplot have been used. The full analysis of the neural network is done by the neural networks toolbox of Matlab, which is based on the FFT responses acquired by a Scilab computational script. The neural network is based on the Levenberg-Marquardt method. We also used the gradient descent method, which allows us to recalculate the delta function based on the weights updated using the decrease of the error and determining the minimum value of the error function (Kazemi et al., 2012). The training method used, trainoss (i.e. one step secant backpropagation), is a more efficient version of the Levenberg-Marquardt algorithm and has allowed calculations requiring a smaller computational cost with better accuracy. For computing the Lyapunov exponents of the system we have also used a Matlab Lyapunov Exponents Toolbox (Matlab LET code [Siu, 1998]).

3.2 Code Validation

In order to validate our numerical code, an asymptotic solution for ε ≪ 1 is proposed. For this purpose, a regular asymptotic expansion method is proposed (Hinch, 1991). Specifically, a regular disturbance parametric method was used. In order to develop the asymptotic solution for the nonlinear governing equation (13) in the limit ε ≪ 1, a convenient change of variable is proposed so that R = (1 + r). The disturbances r around an equilibrium radius R = 1 are assumed to have the form r(t) = εy ₁(t) + ε ² y ₂(t) + 𝒪 (ε ³⁾, here we truncate the regular expansion for terms 𝒪(ε ³⁾. After this expansion we obtain the following system of two linear differential equations 𝒪(ε) and 𝒪(ε ²⁾ governing the quasi-nonlinear dynamic system.

(26)

and

(27)

Note that equations (26) and (27) are both second order ordinary differential equations. The constants that appear in (26) and (27) are given by

(28)

Figure 2 shows a comparison between the 𝒪(ε) and 𝒪(ε ²⁾ asymptotic solutions with the numerical solution (RKS). It is seen that in the limit ε → 0 the asymptotic solutions are exact solutions of the problem.

The insert of figure 2 shows the behavior of each solution for the interval 0 ≤ ε ≤ 0.2. The error between the 𝒪(ε ²⁾ and the numerical solution exceeds 1% for ε ≅ 1 . In addition, the results presented in figure 2 serve as a validation of the proposed algorithm. Therefore, in the next sections we use the developed code in order to explore other regimes that cannot be captured by our asymptotic solution

Figure 2:
Comparison between the asymptotic solutions 𝒪(ε) (filled line), 𝒪(ε ²⁾ (dashed line) and the numerical RKS (filled black circles) for the maximum bubble radius during its oscillatory motion as a function of the forcing pressure amplitude ε.

4 RESULTS

In this section we present different vibrational patterns of the nonlinear motion of the bubble that have been organized in seven possible patterns observed within the range of parameters Re = [0,100], ε = [0,1] for constant values of We = 6, n = 1 and ω = 1. In previous works (Malvar, 2014), it has been proved by a collapse diagram in which Reynolds and Weber numbers were taken into account, that the bubble is more sensitive to Reynolds number when considering collapse and instability. Physically, the increased of the Reynolds number means that more inertia is present in the system when compared to viscous forces. The Weber number, related to interfacial tension, represents the restoring forces acting on the bubble trying to bring its behavior to an equilibrium point. It is easy to realize that for small excitation amplitudes, the bubble presents a harmonic oscillation. However, for a fixed Reynolds number, when ε is increased, the bubble movement responds in a nonlinear fashion even though a harmonic pressure field is applied. This sequence of bifurcations is completely controlled by Re and ε. It is interesting to note that the change of the applied pressure field frequency is related to the system response. However, this change is also coupled with the inertial effects. In the present context, we do not intend to characterize all possible vibrational modes of the nonlinear bubble oscillatory motion, as done in previous works (Malvar, 2014). Our goal is the conception, design and training of a neural network in order to be applied as a tool for identifying patterns and nonlinear responses in the specific flow problem of an oscillating bubble in a Newtonian fluid. This step is quite necessary to extend the approach for investigating bubbles in complex fluids such as non-Newtonian fluids and magnetic suspensions. For this end we propose the classification of different vibrational modes into seven identified patterns based on the time response of the system, on the frequency spectrum and on the phase plane diagram. Each pattern is examined and discussed in details. A bifurcation diagram is also presented as a function of the parameters Re and ε. A bifurcation diagram in these dynamical system analysis is useful in order to identify the oscillating pattern for a given pair of Re and ε. The concept of Poincaré sections and Lyapunov exponents commonly used for analyzing nonlinear responses of dynamical systems is also applied to the investigated bubble system. Actually with these tools, we investigate if a chaotic characteristic is presented in the bubble nonlinear motion. Next, a neural network is designed and trained based on information regarding the first harmonic on the frequency response of the bubble motion.

4.1 Different Patterns of Response

Four of the seven identified patterns are shown in figure 3. This figure shows a typical bifurcation diagram for a fixed Weber number. In this diagram we may see how a combination of Reynolds number and excitation amplitude leads to the emergence of different vibrational patterns. We have omitted patterns 5, 6 and 7 to avoid visual pollution in the image. Based on the distortions of the bubble response on the phase diagrams, the vibrational modes of the bubble system are identified in seven different patterns. 1) periodic with mild distortion in the lower peak; 2) periodic with moderate distortion in the lower peak; 3) periodic with mild distortion; 4) periodic with moderate distortion; 5) without harmonic distortion; 6) non-harmonic and non-periodic; 7) non-harmonic and non-periodic type 2. Here we have proposed that the concept of distortion should be based on a pattern deformation coefficient.

So we define this coefficient as the ratio between the two coexisting peaks shown in figure 4. This ratio determines the degree of distortion in the time response. While values bellow a unit correspond to a slight distortion, a null value means that the answer is harmonic and values greater than the unit represent a moderate signal distortion by our definition. The simulations all assume a polytrophic coefficient n = 1, a forcing pressure frequency ω = 1 and the Weber number set as We = 6. The input parameters that we have varied in order to produce these vibrational patterns are the Reynolds number and the amplitude of excitation. We have decided to adopt this methodology because they are the most relevant physical nondimensional parameters of the problem. The Reynolds number expresses a relation between a diffusive and a convective time scale. Here the term diffusive refers to the diffusion of linear momentum by the fluid's viscosity, while the convective time scale is related to a period of excitation. In this sense, the Reynolds number relates relevant time scales of our dynamical system. The amplitude of excitation was varied in order to increase the degree of nonlinearity of the system. Since we have developed an asymptotic solution for small amplitudes, we decided to explore higher values of ε to study some cases that cannot be solved analytically.

Figure 4 indicates that under the simulated conditions if Re < 5, the bubble response corresponds always to the harmonic pattern 1 regardless of the pressure forcing amplitude. It is also interesting to note that with ε ~ 1 all patterns tend to a harmonic configuration at low Re. So, we identify a critical Reynolds number for which the bubble motion seems to be always harmonic even for ε ~ 1. Now, the stability of the bubble motion is analyzed by examining the phase plane (where g represents ), frequency and time responses of the seven identified patterns.

Figure 3:
Bifurcation diagram. Parameter Re as a function of the pressure forcing amplitude ε, for We = 6.

Figure 5 shows that the phase space for pattern 1 approaches to a stable, periodic and harmonic response. In this case, the bubble response is represented by a figure similar to a circle in the phase diagram, despite a slight distortion in the second quadrant as a direct consequence of the velocity variations on the lower peak for smaller radius. The limit cycle, is also seen as a closed orbit, representing in the phase space a periodic trajectory with a finite and well defined period. In this case, we have an attracting limit cycle. The bubble response on the frequency domain shows that most part of the energy is located in frequencies around ω = 1. This implies a signal (or time series) pretty close to a harmonic response. A second harmonic is presented due to the small deformation in the lower peak of the time response signal.

Figure 4:
Identification of a typical time series (R versus t) of the bubble response.

In pattern 2, the shape anisotropy of the phase diagram is more significant due to higher distortions. A sudden change in speed not only increases the amplitude of the frequency response in the Fourier space, but produces a third harmonic to the system. On the other hand, the phase diagram distortion is not enough to produce an inner limit cycle. In pattern 2, shown in figure 6, the distortion on the phase space suggests the deviation from the harmonic response. This deviation is associated with a higher contribution of inertial forces on the liquid. As shown in figure 3 for a given value of ε, say ε = 0.1, the transition from pattern 1 to pattern 2 occurs for Re = 10 so for Re < 10 the bubble response is represented by pattern 1, while for Re > 10 the bubble responds according to pattern 2.

Figure 5:
Pattern 1 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 5, We = 6,n = 1ω = 1 and ε = 0,1.

The presence of a slight distortion on the time response (pattern 2) allows us to infer about the existence of different harmonics with different energetic levels. Pattern 2 is characterized by different states with nonlinear characteristics but still having one period. In contrast pattern 3 is seen clearly to have a time-series with double-periods that represents a more nonlinear motion of the bubble in comparison to the harmonic one of pattern 1. At higher Reynolds number, even for ε ~ 0.1, we can see a quite nonlinear response of the bubble motion tending to present a continuous spectrum of vibrational degrees of freedom or periods. In this sense, we may state that the oscillating motion of the bubble may be highly nonlinear, even when it is immersed in a Newtonian liquid. We can see similar speed changes on patterns 3 and 4, figures 7 and 8.

Figure 6:
Pattern 2 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 12,We = 6,n = 1, ω = 1 and ε = 0.3.

However, the change in position of pattern 4 is much more expressive, also increasing the amplitude of the fourth harmonic. In particular, the moderate distortion provides a new limit cycle, smaller than the main one. Thus, it is observed that the common generated distortion transforms the stable cycle and attenuates abrupt velocity variations of the bubble surface. In other words, we can say that large accelerations of the bubble surface induce internal extra cycles with different periods. The size of the main limit cycle is proportional to acceleration's variations.

Figure 7:
Pattern 3 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 20,We = 6,n = 1, ω = 1 and ε = 0.3.

The third pattern, shown in figure 7, does not present a distortion over the valley in the time response. The phase space becomes more symmetric, showing that the velocity changes of patterns 3 and 4 are quite similar. However, there is an expressive change in position in pattern 4, which increases the amplitude of the fourth harmonic.

Figure 8:
Pattern 4 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 4,We = 6,n = 1, ω = 1 and ε = 0.5.

Patterns 5 and 6 are shown in figures 9 and 10. It is observed that these patterns can be divided into two categories. Pattern 5 presents a harmonic behavior without distortion. Pattern 6 is not harmonic and non-periodic. Pattern 5 has a tendency of stability in the phase space and only two harmonics in the frequency response. The first harmonic contains almost all the energy of the system. On the other hand, pattern 6 shows a damping effect on the bubble time response and a continuous decrease in its energetic level in the phase diagram. We speculate that this energy loss is associated with an energy exchange between the bubble and the surrounding liquid. This damping effect presented by pattern 6 occurs due to the diffusion of linear momentum by the action of the fluid's viscosity. For a certain range of Reynolds number, the times the surrounding liquid takes to diffuse the energy of the oscillating bubble is much smaller than a period of oscillation. When this phenomenon occurs, the energy related to the oscillating motion of the bubble is quickly dissipated by the surrounding liquid's viscosity.

Figure 9:
Pattern 5 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 10,We = 6,n = 1, ω = 1 and ε = 0.1.

During this bubble-liquid interaction the viscosity of the liquid works as a damper, decreasing the overall energy of the oscillating bubble. In the time response, we note that in each oscillation, the distortion parameter changes, inducing different vibrational responses.

Figure 10:
Pattern 6 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 60,We = 6,n = 1, ω = 1 and ε = 0.3.

Finally, pattern 7, shown in figure 11, presents a low distortion of the peak and a weak aperiodicity. This anisotropic behavior in the shape of the phase diagram is observed, which may be interpreted as a superposition of several degrees of freedom or energy levels. This system is highly nonlinear, showing a different frequency spectrum and phase space according to the variation of the nondimensional parameters. This indicates that in the studies of dynamic systems stability, controlling the physical parameters of the problem is as important as the oscillation amplitude or even the type of external forcing. Viscous and restoring forces in relation to inertia forces are as important on the dynamics of the system as its forcing field.

Figure 11:
Pattern 7 characterized by a time series of R(t) (a), (b) phase plane and (c) Power spectrum given by FFT of the output. The parameters used are Re = 30,We = 6,n = 1, ω = 1 and ε = 0.1.

4.2 Poincaré Sections and Lyapunov Exponents

Another interesting approach used to explore the behavior of our nonlinear system is the Poincaré section of the phase plane. For instance, the construction of the Poincaré map (section) is extremely useful in order to convert the study of the flow in the vicinity of a closed orbit. Roughly speaking, the Poincaré map corresponds to the process of periodically strobing a phase portrait. In particular, we can see from figure 12 that different parameters can lead to an attractor if the collapse criterion is turned off.

Figure 12:
Poincaré section for different patterns. (a) the phase plane and the Poincaré section of an harmonic pattern with low Reynolds and Weber numbers. The transient response is represented by the dashed line. (b) a nonlinear pattern obtained with high Reynolds number and excitation amplitude when the collapse criteria is turned off. Those different parameters leads to an attractor.

However, the results indicate that even a bubble describing nonlinear oscillations its motion will not bifurcate for a chaotic motion if it is immersed in a linear liquid. This type of response is seen even under a high Re and We numbers condition. We might argue, however, that for a chaotic configuration to appear, the bubble should develop oscillations with a period in the same time scale of a complex nonlinear liquid such as a viscoelastic and magnetic fluid (Cunha et al., 2002; Albernaz and Cunha, 2013). In another words, the motion of the bubble should be coupled with the mechanism of structure relaxation (e.g. macromolecules and chains present on these liquids). As the relaxation time of a Newtonian linear liquid is zero compared with any period of the bubble oscillation, the effect of the liquid is related with dissipation and inertia. These mechanisms are not sufficient to produce a chaotic route. These findings have currently motivated us to examine in full details the bubble dynamics in a magnetic anisotropic suspension (Malvar, Gontijo and Cunha, 2016). Actually, in practical application the bubble would collapse before its motion displays a chaotic motion. In order to produce the attractor shown in figure 12 the collapse criteria suggested in Albernaz and Cunha (2013) based on an asymptotic critical radius of collapse was turned off in this specific simulation.

Now, we continue examining the nonlinear response of the bubble dynamics through the concept of Lyapunov exponents of this system. This study allows to observe the sensitivity of the bubble system to initial conditions and also the divergence and convergence rates of nearby trajectories in the phase space. In other words, the Lyapunov exponents (LE) herein called λ are a measure of the exponential separation (eλ) of the neighboring trajectories over all points of a trajectory around an attractor. For stable cycles λ < 0 and neighboring trajectories converge, whereas λ > 0 corresponds to a chaotic attractor. In terms of the Kolmogorov entropy (entropic metric) the LE may be interpreted as a measure of the disorder in the bubble motion response due to the acoustic forcing pressure and its nonlinear interaction with the surrounding liquid. The Lyapunov exponents are also considered to be a dynamic measure of the complexity of delays and may be used for the characterization of chaos and bifurcations. These are standard consequences of a dynamical system with a high degree of nonlinearity and instabilities as can be seen in the unsteady oscillatory motion of a bubble in a complex liquid. LE are found in the present work by numerical computations. It can only be evaluated analytically in some simple cases. Most of them without practical importance in the scenery of nonlinear dynamic systems.

We shall now explore the behavior of Lyapunov exponents for an equivalent autonomous system representing the nonlinear motion of the bubble. Since we are considering the radial motion of a non-deformable bubble in a three dimensional quiescent liquid in the absence of any rotational motion, we have a dynamical system with three translational degrees of freedom. In this sense we can use equation (14) in order to express the primary system as being an autonomous system. In this particular case, in which the bubble is immersed in a Newtonian liquid, Rayleigh-Plesset equation can be just written as an equivalent autonomous system as

(29)

The Jacobian of the system of nonlinear ordinary differential equations given in (29) is expressed by equation (30). Considering different values of Re, We, ε and the initial conditions it is possible to test whether the bubble presents a chaotic behavior or not. For this purpose, we have used the Matlab Lyapunov Exponents Toolbox (Siu, 1998). The output can present different behaviors which are characterized as periodic cycles or stable equilibrium, nonlinear chaos and pure random processes. In the Jacobian approach, the Jacobians are usually found by locally linear mapping in the neighborhoods near the reference trajectory to neighborhoods at a subsequent time (Ataei et al., 2003). Along with the Lyapunov exponents, the Lyapunov dimension is also given for each time step. This dimension, also known as Kaplan-York dimension can be defined as the fractal dimension in which a cluster of initial conditions will neither expand nor contract as it evolves in time (Chlouverakis et al., 2005). The rate of expansion is the sum of the Lyapunov exponents, and this sum will necessarily be negative for an attractor, for instance,

(30)

Figures 13 and 14 show that even though the bubble response in a Newtonian fluid is highly nonlinear, the formal requirements for reaching chaotic configurations of the bubble oscillations are not satisfied.

Actually all Lyapunov exponents are negatives and the Lyapunov dimensions of both systems are null. The Lyapunov dimension is zero when the system is a first order one-dimensional system, such as it occurs in the logistic map or when all the Lyapunov exponents are negative (or zero). All the tests were conducted with no collapse criteria. Under this condition, we can argue that none of the combination of parameters physically accepted would lead to a chaotic motion.

Figure 13:
(a) Lyapunov exponents for the bubble dynamic system, (b) phase space and (c) Poincaré section for Re = 150, We = 2, ε = 0.8 and ω = 8.

Figure 14:
(a) Lyapunov exponents for the bubble dynamic system, (b) phase space and (c) Poincaré section for Re = 200, We = 200, ε = 1.5 and ω = 15.

4.3 Neural Network Training Data from the Bubble Dynamical System

In this part of our dynamic analysis a neural network is designed for identifying the vibrational patterns of the nonlinear motion of the bubble through its FFT response. Previously, we know that the bubble motion is stable. The set of parameter ranges for these tests are Re = [0,100], ε = [0,1] under a condition of constant ω = 1, n = 1 and We = 6. More precisely the input data on the network consists on the frequency f spectrum ranging from 0.2 ≤ f ≤ 1.5. This band in the frequency spectrum includes information regarding the amplitude of the first harmonic together with the spectral spreading behavior on its vicinities.

The trained network uses a set of 36 input matrices called Patterns. Each one with 15 different pairs of information concerning a set of points that accounts for the spectral spreading around frequencies in the interval 0.2 ≤ f ≤ 1.5. The system uses 1000 epochs with a linear transfer function on the first layer of the network and a logarithmic transfer function on the second layer. The first layer is built with 30 neurons and the second with a single neuron. Each training set has an output, called target, referring to its vibrational pattern. The validation set is not used for training. In this sense we test the performance of the network by providing an input with information never seen before during its training. In order to validate the performance of the network and its ability in identifying the vibrational patterns of the bubble nonlinear dynamic we have carried out seven tests. In each test the network uses a different set of weights in its initial condition. For a given input data containing 15 points of the FFT signal in the interval 0.2 ≤ f ≤ 1.5, the output of the network is a real number. This real number corresponds to the pattern the network is supposed to identify. Figure 15 shows the plots of the network performance for identifying the particular patterns.

Figure 15:
(a) Network output as a function of the patterns considering seven tests with different initial weights for each neuron. The symbols represent the Network output for different tests. (b) Average values of the Network output for seven different tests, denoted by the filled black circles. In both plots the full line represents the desired output for each pattern.

Figure 15a shows the output for seven different tests for a given input containing all the possible patterns. The symbols denote the output of the network for different tests. The continuous line would be the expected answer of the network. Since the network uses different initial conditions its output changes for each test. In order to demonstrate the ability of the network in identifying the patterns we might perform an ensemble average over all validation tests. The result of this ensemble average is presented in figure 15b. Here, the expression Network output denotes a real number identified by the trained network that represents the guessed vibrational pattern. Since we are dealing with a nonlinear continuous ODE, in each test the Neural Network may think that a pattern could display some intermediary features of more than a single pattern. In this sense, different tests (with different initial conditions for the Network's neurons) may obtain different responses, as shown in figure 15a. It is seen that the resulting network is able to identify the correct pattern for all the defined patterns within the range of the physical parameters explored in this work. This result indicates that the network has the ability of reading a new set of data that it had never seen before and predict what pattern is responsible for that specific data. The tests point out that using a simple statistical analysis (by considering an ensemble average of 7 tests) the confidence of the network ability for recognizing the vibrational patterns is of 100%.

Identifying the pattern supports us for understanding the nonlinear bubble motion and the influence of the physical parameters on the vibrational modes of the oscillating bubble. Using the collapse diagram available from our previous work (Malvar, Gontijo and Cunha, 2014), we can identify the full range of Reynolds and Weber numbers that leads to the bubble collapse. As a promising application the control of a cavitating bubble by monitoring the bubble collapse time could be done by using neural networks.

5 CONCLUSIONS

The dynamic behavior of a single gas bubble immersed in a Newtonian liquid in the presence of an acoustic pressure forcing was successfully explored and identified with a neural network. The main conclusions are:

• The numerical simulations performed in this work have idendetified several anisotropic shapes in the phase diagram and different modes in the frequency response of the nonlinear bubble oscillatory motion.

• The existence of several vibrational modes has shown a rich dynamical behavior of the bubble-liquid dynamic system explored here. Surprisingly, for bubble nonlinear dynamics in a Newtonian fluid, a chaotic configuration is never reached.

• Attractors resulting from the nonlinear motion in the phase plane in the absence of bubble collapse was largely examined. This analysis was performed by using both Lyapunov exponents as well as Poincaré section methods. The Lyapunov exponents were always negative under all conditions simulated.

• A fractal dimension was not identified in this nonlinear system of bubble dynamics in a Newtonian liquid. However, the lack of a purely chaotic behavior contributes to the stability of the neural network.

• By determining the network parameters, the vibrational patterns were identified and related as being a function of the relevant non-dimensional physical parameters governing the bubble motion dynamics. The present study has successfully shown a rate of 100% guesses of the network in identifying the vibrational pattern of the bubbleresponse system by using numerical data which did not belong to the training set.

• The present work certainly can be considered a first step towards understanding and controlling the nonlinear response of bubbles in oscillatory motions from a dynamical system point of view. Topological transitions in phase plane using of Poincaré sections and the solution in the Fourier space have been identified. We can argue that a chaotic behavior related to a strange attractor could evolve if the bubble is immersed in a complex fluid like a viscoelastic liquid or a magnetic suspension. In this case a drastic change in the bubble oscillations may be possible since a typical period of the bubble oscilations occurs in the same time scale of the particle relaxation time of the ambient liquid. In a future work we plan to continue investigating bubble dynamics in complex fluids using the new approach presented here by incorporating the larger set of dynamical system tools explored successfully in this work.

References

• Albernaz, D. L., Cunha, F. R. 2011. Bubble dynamics in a Maxwell fluid with extensional viscosity, Mechanics Research Communications, Vol 38/3, 255.
• Albernaz, D. L., Cunha, F. R. 2013. Unsteady Motion of a Spherical Bubble in a Complex Fluid: Mathematical Modelling and Simulation, Applied Mathematical Modelling, Vol 37: 8972-8984.
• Ataei, M., Khaki-Sedigh, A., Lohmann, B., Lucas, C. 2003. Estimating the Lyapunov exponents of chaotic time series: A model based method, European Control Conf (ECC), Cambridge, UK.
• Batchelor G. K. 1967. An introduction to fluid dynamics, Cambridge University Press.
• Behnia, S., Mobadersani, F., Yahyavi, M. and Rezavand, A. 2013. Chaotic behavior of gas bubble in non-Newtonian fluid: a numerical study, Nonlinear Dyn (2013) 74:559570.
• Brenner, P. M., Hilgenfeldt, S., and Lohse, D. 2002. Single-bubble sonoluminescence, Rev. Mod. Phys. 74: 425-483.
• Brujan, E. A. 1999. A first-order model for bubble dynamics in a compressible viscoelastic liquid, J. Non-Newtonian Fluid Mech., 84, 83-103.
• Cash, J. R., Karp, A. H. 1990. A variable order RungeKutta method for initial value problems with rapidly varying right-hand sides., ACM Trans. Math. Softw. 16, 201-222.
• Chang, R., Alejos, J. C., Atkinson, D., Jesen, R., Drant, S., Galindo, A., Laks, H. 1999. Bubble contrast echocardiography in detecting pulmonary arteriovenous shunting in children with univentricular heart after cavopulmonary anastomosis, Journal of the American College of Cardiology, Vol 33/7:2052-2058.
• Chlouverakis, K. E., Sprott, J. C. 2005. A comparison of correlation and Lyapunov dimensions, Physica D, 156-164.
• Cunha, F. R., Albernaz, D. L. 2013, Oscillatory motion of a spherical bubble in a non-Newtonian fluid, Journal of Non-Newtonian Fluid Mechanics, Vol 191: 35-44.
• Cunha, F. R., Sousa, A. J., Morais, P. C. 2002, The Dynamic Behavior of a Collapsing Bubble in a Magnetic Fluid, Journal of Magnetism and Magnetic Materials, Journal of Non-Newtonian Fluid Mechanics , Vol 252: 271-275.
• Fogler, H. S., Goddard, J. D. 1969. On the Stability and Collapse of Spherical Cavities in Viscoelastic Fluids, paper number 63A presented at the 62nd. Ann. Meeting A. I. Ch. E., Washington, D.C.
• Geers, T. L., Lagumbay, R. S., and Valisyev, O. V. 2012. Acoustic-wave effects in violent bubble collapse, Journal of Applied Physics 112, 054910.
• Genovese, M., Bauchpiess, A., Brito, J. L. V. de, Doz, G. N. 2001. Damage Detection Using an Hibrid Formulation Between Changes in Curvature Mode Shapes and Neural Networks, CILAMCE - 22th. Iberian Latin-American Congress on Computational methods in Engineering, Campinas. XXII CILAMCE, 1-10.
• Gopal, A. D., Durian, D. J. 1995. Nonlinear Bubble Dynamics in a Slowly Driven Foam, Physical Review Letters, Vol 75, 13.
• Gori, M., Scarselli, F. 1998. Are Multilayer Perceptrons Adequate for Pattern Recognition and Verification?, IEEE Transations on Pattersn Analysis and Machine Intelligence, Vol 20, 11.
• Grossberg, S. 1988. Nonlinear Neural Networks: Principles, Mechanicsms, and Architectures, Neural Networks, Vol 1, 17-61.
• Hinch E.J. 1991. Pertubation methods, Cambridge University Press, Cambridge.
• Kang, K., Oh, J-H., Kwon. C., Park, Y. 1996. Generalization in a two-layer neural network with multiple outputs, Physical Review E, 54/2.
• Kazemi, P., Renka, R. J. 2012. A LevenbergMarquardt method based on Sobolev gradients, Nonlinear Analysis: Theory, Methods and Applications, 75/16, 6170- 6179.
• le Cun, Y., Boser, B., Denker, J., Henderson, D., Howard, R., Hubbard, W., Jackel, L. 1989. Backpropagation Applied to Handwritten Zip Code Recognition, Neural Computation, Vol 1, 541-551. Generalization in a two-layer neural network with multiple outputs.
• Leighton, T. G.. 1994. The Acoustic Bubble, Academic, New York. 20: 638-645.
• Malvar, S. 2014. Bubble Dynamics in Magnetic Fluids: Theory and Applications, Master Dissertation, University of Brasília.
• Malvar, S, Gontijo, R. G. and Cunha, F.R. 2016. On the Nonlinear Motion of an Oscillating Bubble Immersed in a Magnetic Fluid, In preparation.
• Malvar, S., Gontijo, R. G., Cunha, F. R. 2014. Vibrational modes on the nonlinear motion of an oscillating bubble in a Newtonian Fluid using Neural Networks, DINAME - 17th. International Symposium on Dynamic Problems of Mechanics, Natal - RN - Brazil.
• Marmottant, P., Hilgenfeldt, S. 2004. A bubble-driven microfluidic transport element for bioengineering, J. American Chemical Society, Vol 101, Issue 26.
• Pitt, W. G., Husseini, G. A., Staples, B. J. 2004. Ultrasonic Drug Delivery A General Review, Expert Opin Drug Deliv, 3756.
• Press, W. H., Teukolsky, S. A., Vetterling, layerW. T., Flannery, B. P. 1992. Numerical recipes in C: The art of scientific computing, Cambridge: Cambridge Univ. Press, second edition.
• Rayleigh, L. (1917). On the Pressure Develop in a Liquid During the Collapse of a Spherical Cavity, Phil. Magazine, 6/34: 94-98
• Sanaga, S., Jain, A. 2006. A comparative analysis of training methods for artificial neural network rainfallrunoff models, Applied Soft Computing, Vol 6/3, 295-306.
• Schembri, F., Sapuppo, F. and Bucolo, M. 2012. Experimental classification of nonlinear dynamics in microfluidic bubble's flow, Nonlinear Dyn (2012) 67:28072819
• Siu, S., 1998. Lyapunov Exponents Toolbox (LET).
• Soh, W. K., Karimi, A. A. 1996. On the calculation of heat transfer in a pulsating bubble, App. Math. Modell.
• Sorokin, V.S., Blekhman, I.I. and Vasilkov, V.B. 2012. Motion of a gas bubble in fluid under vibration, Nonlinear Dyn (2012) 67:147158
• Street, J. R. 1968. The rheology of phase growth in elastic liquids, Trans. Soc. Rheol., 12: 103.
• Waibel, A., Hanazawa, T., Hinton, G., Shikano, K., Lang, K. 1989. Phoneme Recognition Using Time-Delay Neural Networks, IEEE Trans. Acoustics, Speech, and Signal Processing, vol. 37, 3, 328-339.
• Young, F. R. 1989. Cavitation, Mc-Graw-Hill, London, Chap. 2.
• Young, F. R. 2005. Sonoluminescente, CRC, Boca Raton, Chap. 5.

Publication Dates

• Publication in this collection
  Dec 2016

History

• Received
  18 Mar 2016
• Reviewed
  07 June 2016
• Accepted
  27 June 2016