High Dimensional Parameter Fitting of the Keller–Miksis Equation on an Experimentally Observed Dual-Frequency Driven Acoustic Bubble

A parameter identification technique of an underlying bubble model of an experimentally observed single bubble in a cluster under dual-frequency external forcing is presented. The measurements are carried out via high-speed camera recordings at a rate of 162750 frames per second. The used frequencies during the experiment are 25 kHz and 50 kHz. With a digital image processing technique, the measured bubble radius as a function of time is determined. The employed governing equation for the parameter fitting is the Keller–Miksis equation being a second order ordinary differential equation. The unknown four-dimensional parameter space is composed by the two pressure amplitudes, the phase shift of the dual-frequency driving and the equilibrium size of the bubble. In order to obtain an optimal parameter set within reasonable time, an in-house initial value problem solver is used running on a  graphics processing unit (GPU). The error function measuring the distance between the numerical simulations and the measurement is based on the identification of the maximum bubble radii during each subsequent period of the external forcing. The results show a consistent estimation of both pressure amplitudes. The optima of phase shift and equilibrium bubble size are less significant due to a valley-like shape of the error function. Nevertheless, reasonable values are found that lead to estimations of pressure and temperature peaks during bubble collapse.


Introduction
Acoustic cavitation is the phenomenon where gas bubbles are formed in a liquid by strong acoustic excitation [1]. After the nucleation of the bubbles, they start to oscillate radially and move in space. During their oscillation, these bubbles can collapse so violently, that the local pressure and temperature can become as high as 1000 bar and 8000 K [2], respectively. This phenomenon has attracted the attention of many researchers due to the possible applications in several fields of science. For example, the high pressure can be utilized in surface cleaning or in material erosion [3][4][5][6]. Chemical reactions and radical production can also be enhanced by the strong collapse of the bubbles [2,[7][8][9]. The dynamics of single frequency driven bubbles is mostly understood. During the last decades, it has been extensively studied both numerically [10][11][12][13][14][15] and experimentally [16,17]. Recently, the attention has turned towards dual-or multi-frequency driven systems due to their several positive effects [18,19]. However, the influence of another frequency on the oscillating bubbles or bubble clusters is still not well understood in terms of bubble dynamics. This prevents a coherent interpretation of the partly contradictory experimental results in the literature [20], as well as further optimization of dual-frequency driven sonochemical reactors.
Our group has made a significant step toward the theoretical understanding of dual-frequency driven single bubbles by employing a graphics processing unit (GPU) during the simulations [21]. This allows the investigation of much larger parameter space than it was possible with conventional approaches before (using CPUs).
The main aim of the present paper is to perform experiments to support the numerical simulation, and propose a technique to identify the underlying model of an individual oscillating bubble.
The experimental observation of a single gas bubble in a liquid is not an easy task. One full cycle of its oscillation takes only a couple of microseconds, and its size is in the μm range [22][23][24]. During their radial oscillation, bubbles also move in space [25] which prevents the capture of their motion for a long period. With a high-speed camera, a very high framerate can typically be reached only on the expense of the observable area, from which the bubbles then can move out very fast. In case of single bubbles, two methods exist to stabilize them in space: 1. capture and levitate a bubble in an acoustic trap [26] and 2. bubble production by focused laser light [27,28].
However, in real applications, bubbles usually appear in streamers or clusters; they can act on each other via their emitted pressure wave, merge or disintegrate into smaller bubbles. But still, bubble structures are built up from single bubbles, and their dynamics is the basic building block of the behavior in a larger ensemble.
During the present experimental investigation, the driving frequencies are set to f 1 25 = kHz and f 2 50 = kHz . Observe that the excitation is periodic, but not purely harmonic due to the rational fraction of the frequencies.
The period is defined by the smaller, first frequency component: t f p = = 1 40 1 µs . We do not use any methods to catch and levitate a bubble. Instead, we observe a bubble cluster and try to select and reconstruct the radial oscillation of an individual bubble in the cluster with digital image processing techniques. The bubble is chosen according to the criteria: pronounced collapse ("active" bubble), longer life time, sharp image, and minimum interaction with neighbor bubbles.
Although the above described dual-frequency experiments are mandatory in the understanding of the dynamics of bubble clusters, they cannot provide details on the cavitation activity (e.g. chemical activity in case of sonochemistry) of the individual bubbles. One way to obtain such information is to use numerical simulations fitted to data [7]. The extracted amount of information and their validity is dependent on the complexity of the model itself. The main technical difficulty is the determination of the emerging unknown parameters by comparing the simulated and the measured bubble radii as a function of time via the definition of a suitable error function. This paper proposes such a method employing the well-known Keller-Miksis bubble oscillator [29] that is a second order ordinary differential equation, in which the unknown parameters are the two amplitudes of the external forcing, the phase shift between the two sinusoidal frequency components and the equilibrium size of the observed bubble. The numerical calculations are carried out using an in-house code written in C++/ CUDA C to exploit the high computational capacities of GPUs, and accelerate the identification of the unknown system parameters. Supposing a potentially non-smooth and complicated structure of fitness landscape in the 4D-parameter space, we start with a simple "brute force" approach of test point scanning on a parameter grid.

System set-up and measurement
The measurements were carried out at Georg-August University, Göttingen in a rectangular water tank of inner dimensions 14 cm , in which the water level height was 11 cm. The walls of the container were made of 1 cm thick transparent Polymethyl methacrylate (PMMA) and the base was made of a 3 mm steel plate. The container was sonicated with two piston transducers glued symmetrically to the bottom. The transducers were driven by a function generator (Tektronix AFG 3022) and two power amplifiers (Electronics & Innovation RF-Power Amplifier 1040L and RF-Broadband Power Amplifier 1140LA) via an in-house built impedance matching device. For the observation of the bubble oscillation, a fast camera was used (Photron Fastcam SA5 with Infinity long distance microscope lense and in-house made LED background cw illumination). All the measurements were carried out using fresh tap water.
In order to reduce the dissolved gas content of the liquid, at the beginning of each measurement, only one transducer was operated at 90 kHz with modulating amplitude for eight minutes. Then both transducers were operated at different driving frequencies; namely, one at f 1 25 = kHz and the other at f 2 50 = kHz . The oscillation of the generated bubbles was observed with an exposure time of 1 µs at 162750 frames per second, which gives one frame at about every 6 μs. At this frame rate, the maximum width of the observable area was adjusted to 1.05 mm, which led to a resolution limit of 5.4 µm. Exact size calibration was done a posteriori after each recording with a syringe of 0.5 mm diameter placed in the focus of the camera. Fig. 1 shows a typical frame of the high-speed recordings. Sharp bubbles appear dark in front of the brighter background. Blurred grey spots correspond to bubbles that are out of focus of the camera. Almost all the individual bubbles can be considered as nearly spherical.

Dynamics of individual bubbles under two-frequency excitation
In order to extract the radius of individual bubbles, digital processing was applied using ImageJ software. In this study, only the oscillation of one bubble is presented. After choosing a suitable bubble (according to the criteria above), its close area was cut out from the video. Subtracting the background and applying a black and white threshold on the pictures, the bubble appears black on a white background. An example of the sequence of a bubble oscillation is presented in Fig. 2.
Here the frame size of each subfigure is 130 120 µm µm × . Again, the time difference between every frame is 6 μs.
This sequence lacks the fine details of the oscillation from the collapse to the after bounces that could be observed with an acoustically trapped bubble with higher framerate [30].
The bubble radius on each frame was determined based on the pixel size and the longest distance between any two points along the boundary of the bubble. Fig. 3 a) and b) shows the resulting bubble radius curve as the function of time. Fig. 3 a) presents 351 points of the oscillation (2.12 ms); after that, the bubble moved out of the focus of the camera. The period of the excitation at the used frequency combination is determined by the smaller frequency component f 1 25 = kHz , that is t p = 40 µs . Therefore, 53 periods of the excitation were observed. The solid blue line connects the "period maxima" of the bubble radii, i.e. the overall maxima assumed during every period of the excitation (note that further maxima of smaller value can occur in each period). Fig. 3 b) is a magnification of the measured points during five periods of the external driving. In every period, six or seven data points are obtained. Typical nonlinear bubble oscillations occur in a way that the minima are assumed during very short time (fast bubble wall), while the maxima are approached on a longer time scale (slower bubble wall). Due to this effect and the resolution limit, in Fig. 3 a) and b),   the smallest bubble radii along an oscillation period cannot be determined precisely (or even cannot be observed at all). However, the periodic structure of the oscillation corresponding to the maximum radii at each period is still well reconstructed. Therefore, the comparison of this measured signal with numerical simulations is carried out by taking into account only the respective "period maxima". Fig. 4 depicts the Fourier spectrum of the bubble-radius curve plotted in Fig. 3 a). Clear peaks appear at the excitation frequencies, i.e. at 25 kHz and 50 kHz, above a small noise floor. This indicates a high periodic part of the bubble oscillation, and rather small aperiodic or chaotic parts. The third harmonics of the main driving frequency at 75 kHz appears as well.

The bubble model and the numerical method
For computation, a slightly modified Keller-Miksis equation was used [22]. This second order nonlinear differential equation describes the time evolution of the bubble radius R t ( ) in time: It assumes that the bubble is spherical, and by the compressibility of the liquid, it incorporates sound radiation to first order. The parameters of the liquid are the density ρ L and the sound speed c L . p L is the pressure at the bubble wall and p ∞ is the pressure far away from the bubble which contains the ambient pressure P 0 , and the dual-frequency excitation: Here p Ai and ω i ( , ) i = 1 2 are the amplitudes and frequencies of the excitations, respectively. Their phase difference is denoted by θ . The connection between the pressures inside and outside the bubble at the bubble wall can be written as where the total pressure inside the bubble (left hand side) is the sum of the partial pressures of the non-condensable gas content p G and the vapor pressure p V . The surface tension is σ and the liquid kinematic viscosity is μ L . The gas inside the bubble obeys a simple polytropic relationship [22] where R E is the equilibrium radius of the unexcited system and n = 1 4 . is the polytropic exponent. For the numerical calculations, Eqs. (1)-(4) were rewritten into a dimensionless first order equation system. The dimensionless time is defined as τ ω π = t 1 2 ( ). The dimensionless bubble radius is y R R E 1 = , while the dimensionless bubble wall velocity is y R R E ). The parameters of the water were calculated from the Haar-Gallagher-Kell equation of state [31] at ambient temperature T ∞ = 25 C° and at ambient pressure P 0 1 = bar . The angular frequencies of the excitation were set according to the experiment: ω π 1 2 25000 = ⋅ rad/s and ω π 2 2 50000 = ⋅ rad/s . Table 1 shows the list of the known and unknown parameters of the system.
To find the optimal parameter set corresponding to the best fit of the calculated radius-time curves to the measured data shown in Fig. 3 (in terms of a suitable error function discussed in more details later), several simulations were computed with different parameter combinations. To this end, a 4D parameter scan was performed. Again, the unknown parameter space includes the two pressure amplitudes p A1 and p A2 , the phase shift between the two sine waves θ and the equilibrium radius of the bubble R E . The list of the four parameters, their ranges and resolutions (equidistant distribution) are summarized in Table 2. The dimensionless equation system was treated as an initial value problem and the employed numerical method was a fourth order Runge-Kutta-Cash-Carp method with fifth order embedded error estimation. The algorithm is adapted from [32].
All the possible parameter combinations mean approximately a number of 1.2 million number of initial value problems (IVPs) altogether. In order to reduce the required computational times, the high processing power of a graphics processing unit (GPUs) were exploited. Since the applied numerical algorithm uses only function evaluations to solve the millions of differential equations which are independent of each other, our problem is well suited for parallelization in GPUs. The program code was implemented in C++ and in CUDA C software environment. The interested reader is kindly directed to the website www.gpuode.com, where the program code is freely downloadable. The GPU was a Titan Black card (Kepler architecture, 1707 GFLOPS double precision processing power). Since in bubble dynamics large amplitude oscillation in the collapse-like response region of the system is inevitable, the application of double precision floating point arithmetic was necessary. The final, optimized code is approximately 50 times faster on the aforementioned card than on a four core Intel Core i7-4790 CPU. For further details on GPU accelerated ODE solvers, the interested reader is referred to the publications [10,21,[32][33][34][35][36].
At each parameter combination, the IVPs were solved with initial conditions y R R E representing the equilibrium condition of the unexcited system. After integrating 1024 cycles of oscillation, the properties (maximum bubble radii) of the subsequent 64 cycles of the converged solution were recorded. One cycle means the integration of the system forward in time by τ ω π p p t = 1 2 ( ) . Again, the period of the external forcing is t p = 40 s µ and the first frequency component of the excitation is f 1 1 2 25 = = ω π kHz . With the aid of the GPU accelerated code, the calculations took approximately four hours. Since the minimum values of the measurement do not carry valuable information (see the discussion in the previous section), and the fine structure of the oscillation cannot be reconstructed from the recordings (due to the relatively small frame rate), general methods of parameter fitting of nonlinear ordinary differential equations are not suitable here [37]. Therefore, the comparison presented in this paper only relies on the maximum values of the oscillation.
In the following, let us denote the measured and simulated maximum radii at each cycle of the excitation by  Fig. 3. It must be emphasized that the radial oscillation of the bubble can be chaotic (not to mention the unknown starting point of the simulation in time). In this situation, exact match between the measurement and the simulations cannot be achieved. Therefore, the inclusion of statistics in the determination of the optimal parameter set is necessary.

Results and discussion
As an initial attempt, three error functions are defined as where Δ max , Δ min and Δ avg try to minimize the difference between the measured and the simulated largest, smallest and averaged local maxima in the ith period of the oscillation, respectively. For each error function, the optimal parameter set from the 4D scanned space (where each error function has a minimum) is given in the first three rows of Table 3.
The pressure amplitudes in all the three cases are not significantly different. The values of p A1 and p A2 are  around 0.73 and 0.45 bar, respectively. Therefore, they can be estimated relatively confidently. The equilibrium radius R E , on the other hand, is significantly different for all the error functions. This suggests that the precise identification of the bubble size is a cumbersome task. The phase difference θ has -roughly speaking -two different cases: 0 rad (inphase) and 2.9 rad (antiphase). Fig. 5 a)-c) visualizes the obtained optimal solutions (red curves) together with the measurement (black curves). Each figure shows only a short period of the overall time domain in order to avoid overcrowded subpanels. Moreover, the simulation is also shifted in time so that to match the local maxima as much as possible calculating the cross-correlation function between the two curves. During these calculations the simulations were fixed in time and the measurements were shifted. In each period of the oscillation, after the bubble size reaches its maximum size and shrinks suddenly, the red curve presents the wellknown afterbounces [30]. These high frequency oscillations could not be resolved by the camera at the given frame rate. Nevertheless, the sample rate is high enough to compare at least the local maxima of bubble radius-time curves. Keep in mind again, that with increasing frame rate, the size of the observable area shrinks; thus, the possible measured time domain shortens as well due to the displacement of the bubbles in space. That is, there is always a compromise between the frame rate and the number of the observed bubbles and their recorded time interval. This is another reason why only the local maxima are involved in the definition of the error functions (Eqs. (5)-(8)). Fig. 5 a) corresponds to MIN ∆ max { } in Table 3.
As expected, the largest local maximum of the oscillation fits well to the measured data. However, mostly they show relatively large differences. Here the R E = 12 5 . m µ bubble size seems to be a good estimation as the minimum values of the measurement are lying near the middle in the afterbounces. Fig. 5 Table 3.
The amplitudes in this case are mostly smaller. The parameter set of this solution differs from the previous case only in the equilibrium radius R E . This confirms the observation that this parameter has a strong effect on the oscillation [38,39]. Comparing the average of the maxima and looking for the smallest difference MIN ∆ avg { } , one gets the solution presented in Fig. 5 c). This corresponds to the parameters in the third row of Table 3. Here the phase difference is much Table 3 Optimal parameter values of the error functions (Eqs. (5)-(8)).   Eq. (7) has a minimum at this parameter set.
higher than in the previous two cases. Moreover, the equilibrium radius is much smaller. This curve is periodic, and thus does not fit well in amplitudes to the measurement.
Neither of the three previously discussed error functions can give an acceptable fit to the measurement. Therefore, let us now combine Eqs. (5)- (7). The simplest and straightforward combination would be to look at the weighted average of the three equations. The definition of the new error function Δ w is: The parameter set, where Δ w has a minimum is in the fourth row of Table 3, and its related solution (red) is presented in Fig. 6 a) together with the measurement (black). Fig. 6 b) plots the comparison where the simulation is sampled at the rate of measurement. In both Fig. 6 a) and b), the simulation shows a relatively good agreement (compared to the other three cases) with the measured signal.
Another way to compare the measurement data and the simulation is to re-sort the time series curves according to the period of the driving: that is, the results on the subsequent time domain and plotted to the same diagram. Here i = 0 1 52 , , ,  and t p = 40 s µ is the period of the excitation. Fig. 6 c) shows such a diagram, where the black crosses are the measured values, while the red curves are again the results of the simulation correspondent to Δ w . In this comparison, the simulation still shows a relatively good agreement with the measured signal.
Additionally, one may estimate the maximum pressure p max and temperature T max during the oscillation of the bubble. Based on the adiabatic compression law, these estimates can be calculated with Eqs. (9) and (10): where the subscript 0 means the values at t = 0 , hence R R E 0 = is the equilibrium radius, P 0 1 = bar and T 0 25 = C°, and n = 1 4 . is the polytrophic exponent. In the above case, where Δ w has a minimum, the values are between . These values meet the expectations that approximately below 1 bar excitation amplitude, the bubble oscillation lacks really strong collapses [30,38]. This means that the specific bubble oscillation studied here is probably not able to support effects of strong collapses, such as sonochemical applications. However, the main objective of this paper is to present a technique to estimate the parameters of bubble oscillations.
According to Table 3, the pressure amplitudes can be estimated with high confidence. However, as it was   already highlighted before, the biggest uncertainty in the optimal parameter set is the phase difference of the two excitation frequencies and the equilibrium radius of the bubble. Therefore, it is interesting to look into the effect of these two parameters. Fig. 7 shows four different phase diagrams at constant pressure amplitudes; namely, at p A1 0 70 = . bar and p A2 0 48 = . bar . The horizontal axis is the phase difference θ , and the vertical axis is the equilibrium radius R E . The four diagrams show the error functions Δ calculated via Eqs. (5)-(8), respectively.
The magnitudes of the errors are shown by the color bar next to the subfigures. Darker blue means smaller errors, while green and yellow means higher errors of the respected error function. The red dot represents the minimum place of Δ w defined by Eq. (8), see also the fourth row in Table 3.
Naturally, the best solution that fits to the measured signal would be at a parameter set, where all the four equations are minimal. In Fig. 7, in all the four figures, the structure of the dark blue region (small error) is very similar and appears almost like a sine wave. With changing R E , the phase difference θ also varies along these minimal values.
Consequently, there is a very large set of parameter combinations which almost equally well represent a good fit in terms of the corresponding error function. This is definitely the major difficulty of the present situation; observe that the bubble size can vary by a factor of two (8 μm to 16 μm) without altering the error functions significantly.

Summary and conclusion
The presented work demonstrates the measurement of dual-frequency excited bubble clouds in water and the reconstruction of the model of a single bubble oscillation using the Keller-Miksis oscillator. The free parameters of the fitting of this equation were the two pressure amplitudes of the excitation p A1 and p A2 , the phase difference of the two signals θ and the equilibrium bubble radius R E . With the help of the high processing power of a graphics processing unit (GPU) the numerical solution of more than one million initial value problems (IVPs) took only a couple of hours. This efficiency opens up the door for high dimensional parameter fitting in case of strongly nonlinear differential equations. For parameter fitting, one needs to choose a suitable error function that needs to be minimized. In the present work, the employed error functions were based on the comparison of the difference between the local maxima of the simulated and the measured bubble radii in time. Investigating the results, it has turned out that the pressure amplitudes can be estimated confidently, while the other two parameters (phase shift and bubble size) have high uncertainty.
It must be emphasized that the aim of the present paper is not to find a perfect fit to the measurements but to demonstrate that with small frame rate (compared to the bubble collapse time scale) the investigation of the local maxima of the bubble radii together with the usage of GPUs can be a feasible way to identify the optimal parameter set of the Keller-Miksis equation.
Naturally, there are still many problems that need to be handled. For instance, other type of error functions should be tested by varying the weights in Eq. (8) or including e.g. the standard deviation in the optimization process. Alternatively, measures on back folded radius-time data could be employed. In addition, extracting more bubble radius curves from the same video recording, the fitting could be compared to each other and maybe consistent parameter regions could be selected in case of large, nearly optimal domains, see again the blue area in Fig. 7.