Numerical investigation of the translational motion of bubbles: The comparison of capabilities of the time-resolved and the time-averaged methods

Highlights • In the present study, the accuracies of the time-resolved and time-averaged approaches are compared for modelling the translational motion of acoustic cavitational bubbles in a standing acoustic field.• The investigations are carried out in the parameter space of the driving frequency, pressure amplitude and equilibrium radius by utilizing the high processing power of GPUs.• For positionally unstable bubbles, at sets of parameters, the difference in translational frequency may be more than three and the amplitude of translational motion obtained by the time-averaged method is roughly 1.5 times higher than the time-resolved solution.• When the transient translational motion is important, the time-resolved approach is the proper choice for reliable results.


Introduction
A keen interest in sonochemistry is to utilize ultrasound in chemical reactions. Its physical background is the acoustic cavitation. Due to the high amplitude acoustic waves, bubble clusters consisting of thousands of radially pulsating bubbles are formed in the liquid domain. If the intensity of the irradiation is sufficiently high, the bubbles expand many time larger than their equilibrium size, then violently collapse. The collapse is stopped by the dramatically increased pressure inside the bubble. Near the minimum radius, the temperature can also reach thousands of degrees of Kelvin, which induce chemical reactions in the bubble interior [1][2][3][4][5][6][7] and shock waves are generated [8][9][10][11].
The free radicals produced via the bubble collapse are used in wastewater treatment [12][13][14][15][16], or in degradation and oxidation of pollutants [17][18][19]. For example, to degrade the hazardous chemical species, such as Perfluorooctanesulfonic acid (PFOS), sonochemistry seems to be the only option [20]. Some novel applications of sonochemistry are the production of hydrogen (Green fuel) [21][22][23] and the production of nanoalloys and nanoparticles [24][25][26][27]. In addition, ultrasound can be used to enhance the synthesis of organic compounds or pharmaceuticals by using fewer toxic or toxic-free solvents [28] (Green Chemistry).
The biggest challenge in sonochemistry is to scale up the applications feasible for industrial/commercial size. The main barriers of scalability are 1) the dependence on operating parameters, 2) the high number of involved parameters, 3) the non-uniform cavitational activity, and 4) the attenuation of sound waves [29][30][31]. Many papers deal with the dependence on operating parameters by investigating the nonlinear nature of single bubble [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46]. As the number of the involved parameters is large, the computing power of GPUs is utilized to investigate the wide range of parameter space [47]. To face with the challenge of non-uniform cavitational activity and the attenuation of sound waves, one has to deal with many-bubble systems. One of the key phenomena related to multi-bubble systems is the translational motion of bubbles besides the volume and shape oscillation. Nevertheless, translational motion is a less researched subject, while the spatial location of individual bubbles is of great practical importance to minimize the attenuation of sound waves caused by densely packed bubbles (acoustic shielding [48]) or decrease the volume of passive zones formed in the reactor [49]. From linear theory, it is known that bubbles with size below the resonance size are attracted by the pressure antinode, while bubbles larger than the resonance size are pushed to the pressure node in a standing acoustic field. However, in the nonlinear regime, the repulsion of bubbles from the pressure antinode can be observed. Thus; the bubble accumulation in parallel layers is observed [50][51][52][53][54]. As the bubbles behave as mechanical oscillators, wave dispersion is expected. The increasing void fraction in bubbly media changes the acoustic properties of the liquid, such as the speed of sound and the nonlinear attenuation [31,[55][56][57][58]. For close bubbles, the effect of bubble-bubble interaction is not negligible. Even for a small void fraction, the increasing pressure amplitude and the bubble-bubble interaction has huge influence on the properties of the acoustic field compared to the non-bubbly liquid [59]. The numerical simulations of Louisnard showed that the acoustic energy is dissipated in thin layer of oscillating bubbles [60]. A possible option to moderate the shielding effect is to control the behaviour of the bubble cluster. Control techniques such as a simple ON/OFF control [61,62], the application of bi-frequency driving to avoid the spatial instability [63,64] or the optimization of the geometry [65] already exist in the literature. A possible approach is to apply acoustic manipulation methods to control the position of individual bubbles. Acoustic manipulation methods are already extensively used for solid particles [66,28,67]. However, their application for cavitation bubbles is not yet widespread. Exceptions exist for single bubbles [68], where bubble position control was achieved by shifting the acoustic wave field using opposite transducers with a variable phase shift difference attached at the end of the liquid column. The manipulation of bubbles can be highly challenging due to the nonlinear behaviour of the bubbles and the spatially complex acoustic field generated in a sonochemical reactor [69][70][71][72]. The authors believe that the success of developing a feasible control method for real applications highly depends on numerical simulations. Due to the huge amount of involved parameters, efficient computational methods are required to exploit the capabilities of high-performance computing to investigate the dynamics in detail. For example, in a recent study, employing a six-dimensional parameter space even with a moderate resolution, the number of different parameter combinations was nearly 2 billion [47]. Besides being computationally efficient, an accurate modelling approach is required without oversimplification to obtain reliable results. Thus; the present paper aims to compare the two qualitatively different modelling approaches available in the literature, which are used to describe the translational motion of a radially pulsating bubble.
The translational motion of bubbles is governed by the acoustic radiation forces [73][74][75][76], which are usually classified as primary-and secondary Bjerknes forces. The primary Bjerknes force describes the effect of acoustic pressure gradient on single bubbles, and it gathers the bubbles in certain areas, such as pressure nodes or antinodes in the case of a standing sound field [77][78][79]78]. The bubble-bubble interaction via the emitted pressure is described by the secondary Bjerknes force that can attract or repel neighbouring bubbles [80][81][82]. In addition, the bubble-bubble interaction has an influence on the radial dynamics as well, which was observed numerically [83,59,84] and experimentally [85]. The effect of these radiation forces leads to the formation of stable bubble structures [52,54].
To obtain the individual bubble paths in an acoustic wave field, two kinds of numerical approaches exist. One possible candidate is to describe the translational behaviour of bubbles in terms of timeaveraged approach, which treats every bubble as a single particle accelerated by external forces such as acoustic radiation forces (discussed above), drag force and added mass force. These instantaneous forces exerted on the bubble by the surrounding liquid are averaged over one acoustic period (1/f). Then, the bubble positions and velocities are updated by solving the equation of motion. Although, this approach is justified for weak acoustic fields with small radial oscillation and translation motion; it is able to qualitatively capture the features of caviation clusters [50,52,53].
The second approach is to solve the differential equations describing the instantaneous bubble translation and oscillations. Watanabe and Kukita [86] investigated the translational motion by applying Newton's Second Law directly to a bubble immersed in the liquid domain. However, their model does not involve the direct feedback effect of the translational motion on the radial pulsation. Doinikov [79,87] refined the model, by using the Lagrangian formalism, which includes the direct coupling between the translational and radial motions. The interaction of two bubbles is taken into account via coupled equations of radial and translational motions. Doinikov [80] and Harkin et al. [88] derived coupling terms with an accuracy up to third order and fourth in the inverse distance, respectively. The theoretical model was refined for arbitrary distances between the bubble by taking into account the coupling terms up to tenth order in [89]. The time-resolved approach was generalized for multi-bubble systems with arbitrary spatial arrangement with coupling terms up to third order in inverse distance [90]. Note that the detailed theoretical review of the bubble translation is provided in [81].
The comparison of the numerical techniques has already been given in the literature for a limited amount of parameter sets. Reddy and Szeri [91] found that the decoupling results in an error of approximately 35% and 120% in terms of the averaged velocity in the case of strong collapse in sub-resonant driving and mild collapse in super-resonant driving. Kreftig et al. [92] investigated a broader range of driving frequencies and bubble sizes. An error of 30% in terms of translational speed obtained by the decoupling approximation relative to the translation speed obtained by the full equations was observed even for a weak acoustic field. However, the local error observed in the translational velocity does not indicate the qualitatively different translational dynamics of the bubble. Thus, the present study investigates the accuracy of both approaches numerically on a wide range of parameter combinations, which is not limited to weak forcing.
The numerical results provided in the present paper aid in finding a compromise between accuracy and computational effort. The computations were carried out in plane standing waves. For simplicity, the translational motion of the bubbles were restricted to one spatial dimension. The results showed that both modelling approaches provide reliable results in seeking stable equilibrium. On the contrary, in the case of transient bubble translational motions, the results showed poor agreement. This observation implies that in the case of cluster simulations, where the transient behaviour is inevitable, the time-resolved approach is the proper choice especially when position control is the main aim.

Mathematical model
The radial dynamics of the bubble is described by the Keller-Miksis equation [93]. The translational motion is modelled by two different approaches. The simplest one is to decouple the equation of motion from the differential equation describing the radial oscillation and solve it separately by considering the time-averaged forces acting on the bubble by averaging over one cycle of acoustic period. The approach of averaging and decoupling to calculate the averaged bubble path was applied to qualitatively explain structure formation of acoustic cavitation cluster [50]. In the present paper, this approach is referred as Time-Averaged method. The second approach is to solve the full Time-Resolved equation of motion coupled to the radial oscillation [78,80,87]. The governing equations are summarized below. Note that in the present paper the translational motion is restricted to one dimensional translation; thus, there is only one non-zero space coordinate.

Radial oscillations
The radial dynamics of a spherical bubble in viscous, compressible liquid is described by the Keller-Miksis equation [93] reads as where the extension of the last term takes into account the feedback from the translational motion (ẋ 2 /4). Note that this last term is omitted in the time-averaging approach. In the above equation, R and x are the instantaneous bubble radius and position; respectively. Furthermore, ρ L = 998 kg/m 3 ,c L = 1500 m/s, and p L are the liquid density, the sound speed and the pressure in the liquid at the bubble wall, respectively. Moreover, p ∞ (x, t) denotes the far field pressure, which depends on the acoustic field, see Section 2.4. Note that the dots stand for derivatives with respect to time.
The connection between the inner and outer pressure is described by a mechanical balance of normal stresses at the bubble wall where σ = 0.0725 N/m and μ L = 0.001 Pa s are the surface tension, and liquid dynamic viscosity, respectively. The internal pressure is the sum of the partial pressures of non-condensable gas p Gn and vapour p V = 0 Pa (pure gas bubble). The gas pressure inside the bubbles is calculated according to a polytrophic relationship where γ = 1.4 is the polytropic exponent (adiabatic behaviour), and R 0 is the equilibrium radius of the bubble.

Translation modelled by the time-resolved approach
The translational motion of the bubble is described as [78,86] where m b is the mass of the gas and vapour inside the bubble, F ex (x, t) is the instantaneous external force acting on the bubble and V(t) = 4πR 3 /3 is the bubble volume. After the derivative of the product is expressed, and the equation is reorganized the translational motion of the bubble is described as The mass of the bubble m b is negligible compared to the virtual mass term 1/2ρ L V; thus, after substituting the bubble volume V = 4πR 3 /3 and its time derivative V = 4R 2 πṘ into the equation above, the equation of translational motion is written as [79,81] Rẍ The external force F ex (x, t) is the sum of the instantaneous primary Bjerknes force and the viscous drag force. The primary Bjerknes force, which is induced by the acoustic pressure gradient, is calculated as where ∇p A (x, t) is the pressure gradient at the centre of the bubble. The drag force is obtained according to the Levich formula [94]: Hereby, v ex (x, t) is the liquid velocity induced by the acoustic irradiation calculated at the centre of bubble.

Translation modelled by the time-averaged approach
By neglecting the mass of the bubble m b , and integrating Eq. (4) for one acoustic period T = 1/f, the averaged equation of translational motion is obtained as [92] ρ L where 〈 and F ex (x, t) is the time-averaged external force that is the sum of the time-averaged primary Bjerknes force and the drag force. These forces are calculated as The first term on the right hand side in Eq. (9) is the averaged-virtual mass term. By decoupling this term the equation of motion in averaged and the decoupled form is The significant benefit of the decoupling of volumetric and translational motions is that once the average radial motion is known, it becomes straightforward to calculate the averaged forces. For its simplicity, it is the common technique to numerically simulate bubble clusters. It is worth to be mentioned, in most of the papers, the drag force is calculated according to an empirically obtained formula [50,51,95]; however, in the present paper, for the direct comparison of the numerical models, the averaged drag force is also calculated according to the Levich formula [94]; compare Eqs. (8) and (12). It is worth mentioning that for the adequate treatment of the transient cycles, instead of precalculated tabulated forces as a function of bubble position [51], the averaging was carried out for each acoustic cycle simultaneously during the solution of the Keller-Miksis equation.

Properties of the standing wave field
The acoustic field is assumed to be a standing wave; thus, the far field pressure at the bubble position x is where P ∞ = 1 bar is the ambient pressure and p A (x, t) is the harmonic forcing. This term is written as where P A is the pressure amplitude, f is the driving frequency, k = 2π/λ is the (angular) wavenumber and λ = c L /f is the wavelength. The derivative of the local pressure with respect the time is As the motion is restricted to one dimensional motion, the pressure gradient has only one component that is and the liquid velocity induced by the acoustic irradiation at the center of the bubble is v ex

Numerical technique
For numerical solution, the above equations were transformed into systems of ordinary differential equations by introducing dimensionless variables; namely, the dimensionless bubble radius y 1 = R/R 0 , the dimensionless position ξ = x/λ, the dimensionless time τ = t/(2π/ω) and the dimensionless bubble wall velocity y ′ 1 = y 2 . The dimensionless system of equations in the case of time-resolved and time-averaged approaches are given in A.1 and A.2, respectively. The dimensionless system of ordinary differential equations was solved by using an initial value problem solver.
For calculating simple time-series curves, the models are implemented in Python and solved by using the initial value problem solver called solve_ivp included in SciPy computational library. Depending on the stiffness of the problem, the built-in 'RK45' or 'BDF' methods were used. The high resolution, bi-parametric maps were calculated by using the MPGOS program packages [96][97][98]. MPGOS is an open-source, general-purpose program package, which is developed in C++ and CUDA C and capable to exploit the computing power of graphical processing units (GPUs). Its performance has already been demonstrated in [47,99]. From the available initial value problem solvers, the RKCK45 method was chosen, which is based on the fourth-order Runge-Kutta-Cash-Karp method with fifth-order error estimation. The high processing power of GPUs enables to carry out a considerable amount of numerical simulations within a reasonable time [100][101][102].

Quick overview of dynamical features, linear theory and model comparison
Depending on the bubble size, the excitation frequency and the driving amplitude the typical translational motion of the bubble can be classified [78]. The present section summarizes these different kinds of bubble motions and numerical examples are provided to compare the numerical models. In a weak standing wave, bubbles with size larger than the linear resonance translate to the pressure amplitude minima (node); on the contrary, bubbles with a size below the resonance size are pushed to the pressure maximum (antinode). At the antinodes, the pressure gradient is zero; thus, the primary Bjerknes force vanishes and the bubble position is in equilibrium. The node is also in an equilibrium position because of the zero pressure amplitude; therefore, the bubble is at rest and the time average of the Bjerknes force is zero. As the intensity of the acoustic irradiation is higher, the translational motion of the bubble is more complicated due to the nonlinear behaviour of the bubble. In this case, the large bubble volume expansion and the shift of the collapse phase relative to the driving pressure affects the absolute value and the sign of the Bjerknes force. Thus; the Bjerknes force may point away from the antinode even for small bubbles. In this case, stable positional equilibrium may exist between the antinode and node, where the acoustic pressure causes oscillation of small bubbles that lead to zero mean primary Bjerknes force. Besides stable equilibrium positions, the bubbles may demonstrate positional instability, for two reasons. In the first one, the bubble shape oscillation is parametrically excited by the volume oscillation. This causes the erratic small amplitude dancing motion of the bubble which was observed in [103][104][105][106][107]. In the present paper, this type of instability is not investigated as the models do not take into account the coupled shape mode oscillations [108][109][110][111]. The second type of instability is caused by the change in the sign of the primary Bjerknes force close to the resonance bubble size, where hysteresis of the radial dynamics sets in. Hence, the bubble performs large amplitude oscillation in the space. Fig. 1 demonstrates the different kinds of the aforementioned translational motions obtained in a standing acoustic field. Panels A to E depict the bubble position in terms of the normalized distance from the pressure antinode (x = 0) as a function of dimensionless time at bubble radii R 0 = 20, 60, 100, 110 and 140 μm. In all the cases, the driving frequency and the pressure amplitude are constant, i.e., f = 25 kHz, and p A = 60 kPa. The solid blue and the dashed red curves denote the bubble path obtained by using the time-resolved and the time-averaged method, respectively.
The results presented on the left-hand side of Fig. 1 obey the linear theory [95]. The bubbles with a size below the resonance size of R 0,res = 131.15 μm corresponding to the excitation frequency f = 25 kHz, which is calculated from are attracted by the antinode (panels A and C), while in the case of an equilibrium bubble radius higher than the resonance size, the bubble is attracted by the node (panel E). Although the results obtained by the different models show qualitatively similar behaviour; the dynamics at This is more than a factor of three. In addition, the amplitude of oscillation obtained by the time-averaged model x/λ = 0.126 is roughly 1.6 times higher than the amplitude of translational oscillation x/λ = 0.076 provided by the time-resolved model. From Fig. 1, one can observe that both models are capable to seek for equilibrium bubble positions; there are only small differences in the speed of convergence. The huge quantitative difference observed in the translational frequency implies that the time-resolved model is the proper choice in computations, where continuous translational dynamics plays an important role.

Model comparison in the parameter space of pressure amplitude, equilibrium radius and driving frequency
Investigating simple time-series curves is not feasible to conclude general observation; therefore, numerical investigations were carried out in a wide range of parameters. To ensure the convergence of the numerical solution, the first 25 000 acoustic cycles were treated as the transient solution and discarded. The high number of acoustic periods is justified by the slow translational velocity compared to the radial oscillation, especially for week acoustic driving and small bubble radius as the Bjerknes force is proportional to R 3 , see again Fig. 1A. The results obtained during the next 5000 acoustic cycles were treated as converged solutions and evaluated. In this case, the smoothed bubble trajectory, i. e., the bubble position and velocity at the end of each acoustic cycle was stored. These converged solutions were classified as stable and oscillating translational motion. The convergence to an equilibrium position was tested as the differences between the consecutive position and velocity values were less than pre-defined tolerances The tolerance was set as ε = 10 − 6 . In the case of translational oscillation, the frequency and the amplitude of oscillation were calculated by the Fast Fourier Transform similar to Fig. 1.
To compare the capabilities of the numerical techniques on a wide range of parameters, bi-parametric simulations were carried out in the parameter plane of the equilibrium bubble size R 0 and pressure amplitude p A at fixed driving frequency values f. The investigated set of parameters and their resolutions are summarized in Table 1.
The first column contains the excitation frequency and the second column contains the corresponding resonance radius R 0,r . The third and fourth columns show the parameter limits, and the applied resolution is depicted in the fifth column. At each parameter combination, initial value problem computations were carried out. The initial bubble size was prescribed as the equilibrium radius and the initial bubble wall velocity was set to zero. The velocity of the translational bubble motion was also set to zero. The initial bubble position was varied between the antinode and the next node. In terms of relative displacement, N IC = 5 equally distributed initial bubble positions were prescribed between x/ λ = 0.01 and x/λ = 0.249 to reveal the possible coexisting stable equilibrium positions. The applied limits of the initial displacement is given in the sixth column of Table 1 The detailed comparisons of the models are presented in the following subsections. First, the results are processed in terms of stable equilibrium positions and their stability threshold. Here, only the equilibrium positions (Section 5.1) and the number of equilibrium positions (multi-stability, Section 5.2) are compared; that is, only the final state of the bubble is important. The stable equilibrium positions give no information about the transient behaviour of the bubble; thus, in the case of translating bubbles, the frequency and the amplitude of translational oscillation are compared in Section 5.3. Fig. 2 shows the stable equilibrium positions obtained by the different methods as a function of the equilibrium radius R 0 and the pressure amplitude p A . The first, second, third and fourth rows correspond to excitation frequencies f = 25, 50, 100, and 200 kHz, respectively. The parameter maps on the left-hand side are obtained by the time-resolved approach, while the results on the right-hand side are calculated by means of the time-averaged method. For simplicity, only the equilibrium positions obtained by using x/λ = 0.01 as the initial displacement of the bubble are presented in Fig. 2. The co-existing equilibrium positions revealed by using all the initial conditions are presented in the next subsection. The colourmap denotes the obtained stable equilibrium position in terms of the relative displacement (x/λ). The blue and red regions correspond to the antinodes located at x/λ = 0 and x/λ = 0.5, while the grey domain denotes the pressure minimum (node). In the green domain, the bubble exhibits translational oscillation Table 1 The investigated parameter set, and the range of the applied initial conditions. discussed in details in Section 5.3. The equilibrium positions obtained by the different methods show a good agreement and the gross features in terms of translational stability are captured by both numerical models. Above the resonance size, the bubble path converges to the closest node (grey domain) for each excitation frequency. Well beyond the resonance size, the bubble is attracted by the antinode even for higher pressure amplitude values. From the resonance bubble size, with increasing pressure amplitude, a  Note that the numerical instability may be resolved by using a higher-order numerical scheme to solve the decoupled equation of motion. As usually the semi-implicit Euler method is applied in the literature [50,51], the present paper employed the same method and does not investigate the aforementioned numerical instability in detail. Although, such instability arises at low bubble sizes; it is limited to pressure amplitude values below the ambient pressure (1 bar), see Fig. 2 again. The typical sonochemical applications require pressure amplitude higher than the ambient pressure; thus, in the parameter range of real applications, the averaging method provides satisfactory results for a single bubble, i.e., capturing the attraction of the pressure antinode. In addition, the averaging method was used in particle simulations to capture key features of bubble clusters [50,51,113] and it has been validated for two bubbles experimentally [114].

Multi-stability
To present the co-existence of multiple stable equilibrium positions, the non-oscillating solutions are plotted in 3D in Fig. 4. Keep in mind that the number of initial conditions at each set of parameters is five. In the diagrams, the stable bubble positions in terms the of relative displacement are plotted as a function of the equilibrium radii and pressure amplitude.

Translating bubbles
In the case of translating bubbles, the main oscillation frequency and the corresponding amplitude of the bubble motion obtained by FFT are plotted in Figs. 5 and 6, respectively. Note that the colourmap for the oscillation frequency of the translational motion shows the main frequency in Hz (Fig. 5), while the oscillation amplitude (Fig. 6) is plotted in terms of relative displacement (x/λ). Again, the first, second, third and fourth rows correspond to excitation frequencies f = 25, 50, 100, The side-by-side comparison of the figures shows that the oscillation frequency of the translational motion is 2-3 times higher for pressure amplitude above p A = 0.5 bar in the case of f = 25 kHz and f = 50 kHz excitation frequencies (first and second rows). As the excitation frequency is increased, the difference in the frequency of the translational motion is getting smaller, see Fig. 5. Also with increasing excitation frequency the amplitude of translational motion (Fig. 6) decreases for both methods and the results are getting closer. To sum up, for translating bubbles, the application of the time-resolved method is mandatory to capture correctly the bubble path for excitation frequencies below f = 50 kHz.
An example of the difference in the oscillating solution is presented in Fig. 7, where the dimensionless bubble radius R/R 0 , the relative displacement x/λ and the average Bjerknes force F B1 are plotted as a function of the dimensionless time. The red and blue curves denote the results obtained by the time-averaged and the time-resolve methods, respectively. Note that in the case of the time-resolved method, the instantaneous forces were calculated during the integration, and then the results were averaged for each acoustic cycle. The results on the leftand right-hand sides are obtained for bubble radii R 0 = 80 μm (below resonance) and R 0 = 100 μm (near the resonance), respectively. The excitation frequency is f = 25 kHz and the pressure amplitude is p A = 0.75 bar.
By comparing the diagrams, one can see that below the resonance size of the bubble, the radius-time curves show a relatively good agreement. The primary Bjerknes force is a negative number; thus, the bubble approaches the closest antinode at x/λ = 0. Only a small difference in the speed of convergence can be observed. This is in accordance with the results presented in Section 5.1 and in Fig. 2. Near the resonance; however, the fluctuating pressure amplitude (due to the alteration of the bubble position) induces complicated bubble oscillation and translational motion [78]. Such behaviour are not captured by the time-averaged models; thus, the radius-time curves and the Bjerknesforce time curves show significant differences.

Summary and discussion
The translational motion of cavitation bubbles is responsible for various forms of bubble clusters. Depending on the acoustic field (standing wave, travelling wave) different distributions of bubbles have been observed experimentally, such as the cone at the top of the transducer [113], or streamers [50][51][52]54]. The common property of the observed clusters is that the bubbles do not fill the liquid container evenly; thus, it lets a massive amount of passive zones in a sonochemical reactor. These zones results in non-uniform cavitational activity, which is one of the main barriers to the scalability of sonochemical applications. Another limitation is the attenuation of the sound waves (acoustic shielding [48]) due to the densely packed bubble clusters near the irradiation surface, such as at the top of an immersed transducer. Due to the high energy dissipation near the irradiating surface, the penetration depth of the irradiation is approximately a few centimetres [49].
The simple idea of increasing the irradiation intensity to face the scalability of sonochemical reactors fails due to the aforementioned reasons. As the intensity of the irradiation is proportional to the square of the pressure amplitude (I∝p 2 A ); a small increase in the pressure amplitude decreases the energy efficiency of the reactor. The high pressure amplitude induces strong collapse; however, there is an optimal range of bubble collapse strength, which depends on the application [3]. In addition, the higher the pressure amplitude the more the number of the generated bubbles, and the resulting higher density of bubbles enhances the acoustic shielding effect as well. The possible solution is to optimize a sonochemical reactor by manipulating the acoustic field to achieve homogeneous bubble distribution. As the number of involved parameters is high, and the dynamics of a bubble show nonlinear behaviour, the numerical simulations are important tools in the field of sonochemistry. In a high dimensional parameters space by applying a moderate resolution for each parameter, the number of parameter combinations can be extremely high. An example for the investigation of a single bubble in dual frequency covered nearly 2 billion parameter combinations [47]. Therefore; to develop reliable control techniques, an efficient but reliable numerical model is required. The accurate models can be widely used to properly describe bubble dynamics in applications where translational motion is essential. Besides the sonochemical applications, translational motion is important in medical applications as well. In the case of clot lysis [115,116], the ultrasound-stimulated microbubbles are pushed close to the clot.
In general, the modelling of the translating and oscillating bubble requires solving the governing equations describing the radial oscillation and the equation of motion for the translation. There are two main approaches exits in the literature. The trivial technique is to solve the system of the coupled ordinary differential equation describing the translational motion and radial oscillations, which is referred to as the Time-Resolved method in the present paper. The second approach is to decouple the radial oscillation and the translational motion. In this case, the forces acting on the bubble are averaged for one acoustic cycle and the translation of motion is solved by using the averaged values. Although, studies have already been published investigating the validity of the decoupling approximation [91,92], they are limited to a few parameter combinations or weak acoustic forcing. The present study compares the results provided by both modelling approaches on a wider range of parameter combinations, which is not limited to small amplitude oscillation. Therefore, qualitatively different dynamics were identified in the oscillatory regime.
The results showed that in terms of equilibrium positions, both models give the same results. Below the resonance size, the bubble is attracted by the antinode, and above the resonance size, the bubble tends to approach to the node. This behaviour of the bubble is well known for weak standing acoustic fields, which is described by the linear theory [81]. In addition, it was observed that close to the higher harmonic resonances, intermediate equilibrium positions between the node and the antinode may exist. Contrary to the good agreement in terms of translational equilibrium, there are significant differences for continuously translating bubbles. In this case, the bubble oscillates if the bubble size is close to the resonance size. It must be emphasized that the models contains simplifications. In the case of the low frequency and high amplitude region (Giant Response region [44]), the bubble wall velocity calculated by the Keller-Miksis equation may exceed the speed of sound The graphs on the left-and right-hand sides are calculated for equilibrium radii R 0 = 80 μm and R 0 = 100 μm, respectively. in the liquid. Yasui [117] proved that based on the fundamental theory of fluid dynamics, the bubble wall Mach number never exceeds M = 1; thus, in this case, the correction of the model describing the radial dynamics would be necessary [118]. In addition, the equation of motion describing the translational motion neglect the effect of the complex flow around the bubble. Corrections are available to take into account the effect of the drag force acting on the bubble by adding the history force depending on the value of the translational Reynolds number Re t = R|u|/ν and the radial Reynolds number Re r = R ⃒ ⃒ ⃒Ṙ ⃒ ⃒ ⃒ν [119][120][121]. Reddy and Szeri [91] proved that for violent bubble collapse, when Re t > 1 or Re r > 1, the history force is negligible. Further improvement in the model accuracy can be achieved by including the history force.
In the case of a single bubble, the complexity of both models are quite simple. However, for bubble clusters containing thousands of bubbles, the interaction of individual bubbles via their emitted pressures (e.g., via the secondary Bjerknes force [122,113,50,114,51]) has to be taken into account. The secondary Bjerknes force can be attractive or repulsive for a pair of bubbles depending on the equilibrium radii and the properties of the acoustic field that complicates the prediction of the cluster behaviour. The calculation of the secondary Bjerknes force for bubble clusters by using the time-averaged approach includes simplifications such as equally varying bubble volume; thus, the complexity of the model does not increase significantly. On the contrary, in the case of the time-resolved model, the differential equations include the coupling terms up a given orders [80,89,88] in terms of inverse separation distance. These equations contain implicit terms that emerge numerical difficulties, i.e., at every timestep the evaluation of the right-hand side of the ODE requires solving a system of linear equations, which makes the efficient implementation of the models challenging (for example in CUDA to exploit the processing of power of GPUs). An important advantage of solving the full ODE instead of averaging is to take into account the time delay caused by the compressibility of the liquid in the bubble-bubble interaction [123,124]. The time delay may have a huge influence on the radial dynamics and thereby the translational motion as well, especially in the transient regimes. The explicit inclusion of the time delay makes the simulations computationally intensive as the delay differential equations have large memory requirement, which limits the size of the bubble cluster in simulations [83]. A promising approach is to decompose a model into subclusters and handle the interaction between subclusters [125]. The present results obtained in the case of a single bubble imply that the precise prediction of bubble-cluster behaviour requires the solution of the time-resolved approach [90], in spite of the arising numerical difficulties.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 1 = y 2 , the dimensionless position ξ 1 = x/λ and the dimensionless translation velocity ξ 2 = ξ ′ 1 , where the prime denotes the derivative with respect to the dimensionless time τ = t/(2π/ω). The dimensionless system in the case of the time-resolved approach is written as The local pressure p A and its time derivative ṗ A are the functions of the dimensionless variables p A = P A ⋅cos(2πξ 1 )⋅sin(2πτ), (A.7) ṗ A = P A ωcos(2πξ 1 )⋅cos(2πτ).

(A.8)
The external force is the sum of the instantaneous primary Bjerknes force F B1 and the drag force F D . These forces are calculated as Note that in the above equations C 0 − C 24 are precomputed constants, which are independent of the time or position of the bubble. These constants are summarized in A.3.

A.2. Time averaged method
The time averaged forces were calculated by changing the variable (A.13) By using the above notation, the averaged forces are calculated as The bubble position and velocity is updated according to the semi-implicit Euler method. Note that the timestep of the translational motion is one acoustic cycle which is one unit in the case of the dimensionless system; thus, the translational motion is updated as ξ 2,j = ξ 2,j− 1 + ξ ′ 2,j , (A.17) ξ 1,j = ξ 1,j− 1 + ξ 2,j , (A. 18) where j counts the number of acoustic cycles. The integral terms appearing in the above equations are calculated numerically by integrating the functions simultaneously with the dimensionless Keller-Miksis equation. Therefore; the following set of dimensionless ordinary differential equations are solved for each acoustic cycle were calculated according to the formulas described in the previous section, except N. In the case of N the feedback term from the translational motion (ξ 2 2 C 9 /4) was neglected.