Simulations of a full sonoreactor accounting for cavitation

Highlights • A full sonoreactor, transducer/cavitating liquid/vessel walls, is simulated.• The US generator’s automatic frequency selection is accounted for.• The current feeding the transducer is the unique model input.• The dissipated electrical power is compared to experiments.• Transducer immersion depth strongly influences the power dissipated.

In the 1980s and 1990s, there was a rapid resurgence of interest in cavitation and sonochemistry applications, mainly at the laboratory scale. As any growing science, the difficulty its extrapolation to the industrial scale inevitably occurred, and raised critical issues. The phenomena occurring within a sonoreactor are very complex, and unfortunately, unlike more traditional branches of process engineering, there is currently hardly any tool to design sonochemical reactors. One reason for this is the great physical complexity of acoustic cavitation. The enormous range of spatial and temporal scales involved, for example, makes dimensional analysis impractical. On the other hand, because ultrasound is employed, acoustics should be one of the most significant ingredients of sonoreactors science. Nevertheless, it remains one of the most disregarded in prior research.
To shed new light on this issue, several groups have attempted to simulate the acoustic field inside sonoreactors in order to predict cavitation events within the reactor see [25] for a review.
Most simulations rely on linear acoustic [26][27][28][29][30][31][32][33][34][35][36][37][38][39][40], which takes the form of an Helmholtz equation and can be easily solved by the finiteelements method (FEM). The main drawback of such an approach is that the effect of cavitation on the acoustic field is not accounted for, apart for some studies using an empirical attenuation coefficient. In order to account for the presence of bubbles, nonlinear models of bubbly liquids acoustics should be used [41][42][43], but the latter are difficult to solve on large spatial scales, despite interesting results have been obtained on reduced scales by using involved numerical methods [44,45]. Several groups [46][47][48][49][50][51][52][53][54] have used the linearized version of the Caflisch model developed by Commander & Prosperetti [55]. As the latter assumes linear oscillations of the bubbles, this approach yields unrealistic values of the attenuation coefficient where inertial cavitation is involved [56] (see also Ref. [57] for an in-depth discussion).
Bubbles are the main dissipators in the liquid. Each bubble dissipates mechanical energy into heat and the sum of these contributions is the power measured by the calorimetric method. The above-mentioned studies cannot therefore compute correctly the calorimetric power, either because bubbles are neglected, or because the energy they dissipate is largely underestimated by linearized bubble models [58]. To get rid of the latter limitation, Louisnard [59,60] developed a reduced nonlinear model of sound propagation in cavitating liquids, based on Caflisch's model and accounting for fully nonlinear bubble dynamics. As the proposed nonlinear Helmholtz equation directly involves the power dissipated by each bubble, computing the power dissipated in the liquid becomes natural within this model. This model has been successfully used to predict commonly observed cavitation phenomena [60], such as conical bubble structures under ultrasonic horns [61] or flare structures that are very similar to those observed experimentally in ultrasonic cleaning tanks [62]. Several refinements have been proposed since then, to account for radiation losses [63,64]. Sojahrood and co-workers have extended the theory in various aspects [65], including extension to coated microbubbles [66,67]. The latter group also recently performed experiments with a layer of almost monodisperse coated microbubbles to assess experimentally the attenuation and sound-velocity in the bubbly liquid and found good agreement between their extension of our model and the experimental results [68]. Both attenuation and sound velocity were found to vary with acoustic pressure, and to our best knowledge, this constitutes the first experimental results supporting the concepts raised by Louisnard in Ref. [59]. Other groups have recently used Louisnard's model for different fluids other than water [69], in liquid aluminium [70], at frequencies larger than 20kHz [71]. Interesting results have been obtained by Lesnik and co-workers, who coupled our model with the equations governing the translational motion of bubbles and their effect on the fluid motion [72].
A sonoreactor consists of a volume of liquid enclosed within solid walls, in contact with one (or more) transducer, which in turn is driven by an electrical generator. Using the reasonably realistic model for the cavitating liquid of Ref. [59], simulating the full assembly consists in modelling the transducer and coupling its vibrations with the one of the liquid. Modelling transducers has now become an easy task thanks to the implementation of the piezo-electricity equations in most FEM codes, in our case COMSOL Multiphysics. Moreover, the equations of mechanical coupling between the solid vibrations and the acoustics of the liquid are classical [73] and can also be used to study the recipient wall vibrations [31,57]. Therefore, as far as the internal design of the transducer is precisely known, and circumventing the additional difficulty brought by the automatic frequency selection implemented on modern ultrasound generator, modelling a full sono-reactor under COMSOL becomes feasible. Demonstrating this assertion is the first objective of this paper.
In order to test our numerical predictions, we also performed a significantly large set of experiments by varying several controllable parameters: the input current amplitude, which is the quantity experimentally controlled by the generator level button, the geometry of the enclosing vessel, and the immersion depth of the transducer.
The paper is organized as follows: Section 2 describes the experimental setup, the set of experiments realized and the principles underlying the frequency control loop of the generator. Section 3 describes the model. We first recall the main lines of Louisnard's model [59], give a summary of the equations used for solid vibrations, including piezoelectrics, and how we account for the automatic frequency control. Section 4 shows the comparison between numerical and experimental results, for the two set of experiments, varied current and varied immersion depth. The agreement between the simulated and experimental results is discussed in Section 5, along with suggested future enhancements of Louisnard's model.

Experimental setup
A schematic of the experimental setup is shown in Fig. 1. The apparatus used in this study was a 20kHz homemade standard Langevintype transducer (SinapTec, Lezennes, France) with two piezoelectric rings pre-stressed by a steel screw between a mass and counter-mass, both made of titanium alloy. The exact resonance frequency of the unloaded transducer, f 0 = 20342Hz, was measured with an impedancebridge. The transducer was driven by a computer-controlled ultrasonic generator (SinapTec NexTgen Inside 500), by which several electrical quantities, such as power, frequency, voltage and current and their mutual phase, impedance etc, can be monitored and logged every 50ms during experiments.

Sets of experiments
The transducer was immersed in the center of three different commercial glass beakers (hereinafter referred to as A, B and C) filled with tap water. The dimensions of the beakers are shown in Table 1.
For each beaker, two additional geometrical parameters can be varied: the transducer immersion depth h T and the water level h liq (Fig. 1). In this study, two sets of experiments were carried out. In the first set, both the transducer immersion depth h T and the liquid height h liq were kept at constant values for each beaker, whereas the input current of the transducer was varied. For beaker A, the transducer immersion depth was arbitrarily fixed at h T = 3 cm and the liquid height h liq at 12 cm. For beakers B and C the values of h T and h liq were chosen so that the two ratios h liq /H and h T /H remain beaker independent. The corresponding values of h T and h liq can be found in Table 2. The input RMS current provided by the generator was varied from 0.08 to 0.8 A in all experiments of set 1.
In the second set of experiments, the current was fixed at a constant value of I 0,RMS = 0.6 A and h T was swept in a beaker-dependent range (Fig. 2). For each beaker, the liquid volume V liq was kept constant so that the liquid level h liq varied in function of immersion depth h T following:
where D int = D − 2d is the beaker internal diameter and D T is the transducer diameter. Finally, for each set of current/ transducer depth, the liquid was irradiated continuously for only for 30s to avoid heating, at ambient temperature and pressure, in all experiments. As the generator internal control loop was able to set the frequency in less than 1s, all the electrical data monitored during the first two seconds were removed. The remaining data were averaged for each experiment, and the standard deviation of the fluctuations was calculated, yielding an error bar for each point. In summary, all the measured electrical quantities presented hereinafter correspond to averages/standard deviations over 28s with data points logged every 50ms.

Frequency control loop
When the transducer is immersed, it operates against a mechanical load and is therefore in contact with a non-zero complex mechanical impedance. This has the effect to slightly shift the resonance frequency and possibly results in a reduction of the vibration amplitude [74]. To circumvent this problem, ultrasonic generators usually employ some control strategy to ensure a constant vibration amplitude of the ultrasonic tool by adjusting in real time the frequency of the system to the load.
In the vicinity of the resonance frequency, the vibration behavior of the transducer can be described by the equivalent circuit displayed in Fig. 2 [75]. The capacitive branch C b accounts for the purely dielectric properties of the piezo-electric rings, whereas the motional branch in blue (RLC circuit in series) reflects the mechanical properties of the vibrating element and the load. It is important to note that all the energy dissipation mechanisms are concentrated in the resistance R a = RZ BM . The latter include energy dissipated into heat in the transducer, but more importantly, all the energy dissipated in the load (for example by cavitation). The branch C b is classically compensated by a parallel matching inductance so that only the motional branch remains at resonance. The control strategy of the generator consists in estimating the complex motional impedance Z BM = V/I BM and to adjust the frequency so that its reactive part vanishes (IZ BM = 0). This ensures that the whole available input power is transferred to the load, whatever the fluctuations of the latter. The additional control of the input current amplitude allows a good stability in operating conditions.

Numerical study
The model detailed hereinafter was implemented in the commercial code COMSOL Multiphysics, based on the Finite Element Method (FEM). The numerical modeling was performed in the 2D-axisymmetric domain and the sketch of the modelled setup is presented in Fig. 3a. Similar computation could be easily extended to 3D geometries if required, at the expense of heavier computational resources.

Cavitating liquid
In order to account for acoustic energy dissipated by cavitation, the liquid was modelled by a previously published reduced nonlinear model [59], based on the rigorously derived Caflisch equations [43] describing acoustics of bubbly fluids. Our model results in a nonlinear Helmholtz equation, whose squared wavenumber is expressed in function of the energy dissipated by cavitation bubbles by thermal diffusion and by viscous friction in the liquid along their radial oscillations around the bubble. The latter energies are estimated by solving a fully nonlinear bubble dynamics equation. The nonlinear Helmholtz equation is used only in zones where inertial cavitation takes place, i.e. where the

Fig. 2.
Equivalent circuit for the transducer [75]. The capacitor C b originates from the pure dielectric character of piezo-electrics. The motional branch reflects the mechanical load of the transducer. acoustic pressure exceeds the Blake threshold. Assuming monoharmonic oscillations and setting the pressure field as p ( , the equations of the model can be summarized as: where k l = ω/c l , P B is the Blake threshold [76][77][78]: is the dimensionless Laplace tension. ω is the angular frequency, c l the sound velocity in pure liquid, p 0 the static pressure, σ the surface tension, R 0 bubble ambient radius. The nonlinear wavenumber k NL is defined by: where nance frequency, and Π v (|P|) the viscous dissipation function. The latter is pre-computed by solving a Keller equation augmented with a reduced model accounting for the heat and water transfer between the bubble and the surrounding liquid (see A), fitted, and injected in COMSOL as an analytical function [59,57]. The nonlinear Helmholtz equation is implemented within the Pressure Acoustics COMSOL physics module. For all the simulation results presented hereinafter, we considered air bubbles of ambient radius R 0 = 5μm in water at ambient temperature and pressure, and the bubble density was set to 50bubbles/mm 3 , unless otherwise specified. This choice is justified by experimental results which shows that at 20kHz, the ambient bubble sizes lie in a narrow range between 2μm and 10μm [79,80]. This is moreover consistent with theoretical results on bubble shape instabilities leading to bubble breakup above some critical size [81][82][83][84][85][86]. The order of magnitude of tens of bubbles per mm 3 for N is mentioned in [62] and the selected value N = 50bubbles/mm 3 yields void fractions β = 4/3πR 3 0 N ≃ 2.6 × 10 − 5 , which is of the order of magnitude of results reported in Ref. [79].
The computational domain and the associated boundary conditions (BC) are illustrated in Fig. 3b. On the liquid surface, a soft boundary condition (P = 0) was set. The conditions to be applied on the liquid boundaries in contact with solid parts are deferred in the next section.

Transducer and vessel wall
The motion of linear elastic materials (mass, screw, counter-mass, and vessel boundaries) can be represented by Newton's equation complemented with Hooke's law [31]. In time-harmonic form at frequency ω, the equations read: where U S is the mono-harmonic displacement field amplitude, T the stress tensor, S the strain tensor and c the elasticity tensor. The latter equations are available in the basic version of COMSOL within the Solid Mechanics physics module. Extension to the motion of piezo-electrics is classically performed by adding a constitutive electrostatics equation cross-coupled with Hooke's law, here in stress/charge form: where D is the electric displacement field, E = − ∇V the electric field. The system is closed by solving Gauss law ∇⋅D = 0 in the piezo-electrics. The latter equations are available in COMSOL additional Acoustics module: extension of Hooke's law for piezo-electrics is available under an optional sub-node of the Solid Mechanics physics module, Gauss law is implemented in a Electrostatics physics module, and both modules are coupled through a Piezoelectric Effect multi-physics node, implementing Eqs. (9) and (10). The Frequency domain solver computes the complex displacement field U S in all materials plus complex electric potential V in the piezoelectric material. The electric displacement vector field D can further be deduced from Eq. (10) and the input current at the junction between two piezo-electric rings can be computed as: where n 12 is the normal unit vector pointing from piezoelectric disc 1 towards piezoelectric disc 2.
Isotropic losses in materials are considered, in order to account for mechanical energy dissipated in the transducer. A simple way is to use a complex Young modulus Preliminary simulations showed that dissipation occurs mainly in the titanium mass and counter-mass, and in the steel pre-stress screw. Piezoelectric rings would also suffer losses by various mechanical, electrical or electro-mechanical mechanisms, but their contribution to the power dissipated was neglected in this study. The mechanical boundary conditions on all solid parts exposed to air are free conditions (T n = 0). For the electric boundary conditions of piezo-electrics, the lateral sides are assumed free of surface charge (D.n = 0), the faces in contact with mass and counter-mass are grounded (V = 0). The classical boundary condition on the interface between the two piezo-electric rings would be an imposed known complex voltage.
Here however, to mimic the generator behaviour, we need to impose the input current I. To do so, the boundary between the piezo-electrics is assigned to a voltage V, whose value is defined as the solution of a global equation constraining the complex input current I, computed from Eq (11), to the required value I 0,RMS ̅̅̅ 2 √ . Finally, Fluid/Structure Interface Conditions must be specified at boundaries separating vibrating solids and fluid (labelled as FSIC on Fig. 3b). These conditions classically express continuity of velocity and stress at the boundary [73,31] and are easily implementable in COMSOL by an Acoustics Structure Boundary multi-physics node.
The complex electrical impedance of the transducer Z Trans = V/I can be easily computed in any conditions. The impedance curve of the unloaded transducer can be obtained by performing a frequency sweep in a range around the resonance frequency. Since the mechanical properties of metals can slightly vary between different samples, they were fitted by using the experimental impedance curve of the unloaded transducer: the Young modulus of titanium was adjusted so that the computed unloaded resonance frequency matched the experimental one, and its loss angle tanδ so that the computed and experimental impedance module |Z Trans | at unloaded resonance matched. The values found by this method (Table 3) are in agreement with commonly reported values [87]. The loss factor for the steel pre-stress screw is set arbitrarily to tanδ = 5 × 10 − 4 , a value representative of results reported in the literature [88]. This procedure ensures that the simulated transducer shows exactly the same unloaded impedance curve as the one used in experiments.

Modelling of the control loop
Since the behavior of the piezoelectric transducer is accounted for in the model described above, the input parameters for a given experiments are the same as the experimental ones: frequency f and driving current amplitude I 0,RMS . However, as described in Section 2.3, frequency is not known by advance and is adjusted by the generator to cancel the imaginary part of the motional branch impedance IZ BM = 0. The working frequency is thus load-dependent and this frequency selection method must be reproduced to perform simulations as close as possible to the experimental conditions.
To do so, the static capacity C b , which is load independent, is first deduced once for all by a fit of the unloaded impedance curve. Then, for a given configuration (f, I 0,RMS , geometry of the liquid), the electrical impedance can be calculated and, knowing C b , the motional impedance Z BM can be deduced by: We thus proceed as follows: for each configuration, several simulations are performed for a range of step-wise increasing frequencies, close to, and slightly below the resonance frequency. All quantities of interest, especially Z BM , are computed each time. The frequency sought is the one for which IZ BM crosses zero and is determined by inverse interpolation on the (f, IZ BM ) curve. This allows an exact mimicking of the generator strategy, at the price of performing multiple simulations for a given configuration. More details and illustrations can be found in Ref. [57].

Variable current
In this set of experiments, the liquid geometry is kept fixed and the input current amplitude is varied. The comparison between simulations and experimental results is shown in Fig. 4. The operating frequency of the system chosen by the generator is displayed in Fig. 4a: it can be seen that the simulation correctly captures the frequency evolution with current, which asymptotes to a value about 30Hz below the unloaded resonance frequency, indicating that our model produces a reasonable estimate of the load seen by the transducer despite a slight difference of approximately 50Hz between numerical and experimental results (which is less than 0.25% error).
The phase of the transducer impedance is shown in Fig. 4b. The  numerical curve reasonably reproduces the experimental result, with a minimum value near 0.2 A. The descending part of the curve corresponds to a transition from a constant phase (leftmost part of the curve), when there is no cavitation and linear acoustics holds, up to the progressive buildup of a cavitation cloud at the transducer tip. Fig. (4b) also display the experimental and computed phase of the motional branch impedance. Apart from small fluctuations of about 1 • , Z BM is close to zero, which clearly demonstrates that the generator correctly ensures a purely real impedance of the motional branch. Fig. 4c shows that the shape of the transducer impedance modulus vs. current is also correctly predicted by the model, except for very low currents. So far, we don't have a clear explanation for this mismatch, but within this range of low currents, the acoustic impedance seen by the transducer varies considerably, rendering frequency control more difficult. Furthermore, the generator does not allow the use of very low input currents, making a further comparison in this area difficult. Finally, in Fig. 4d, the computed and experimental electrical power consumed are compared. Our simulation, as can be seen, overestimates the electrical power by a factor of about two. Furthermore, the experimental curve is almost linear with increasing currents, whereas the numerical curve bends upwards. This suggests that our model overestimates the energy dissipated by cavitation in the liquid.
In order to examine the sensitivity of the results to the model free parameters, computations were repeated for R 0 in the range 2μm to 10μm and N 0 in 20bubbles/mm 3 to 110bubbles/mm 3 . The same conclusions can be drawn, apart from very slight changes in the other output variables, as shown in Fig. 5 and 6, respectively.

Variable transducer immersion
The influence of the beaker geometry and the immersion depth of the transducer were investigated. The three different-shaped beakers formerly labeled as A, B, and C were used in experiments for various transducer immersion conditions, and the corresponding acoustic fields were computed with our model. The input current of the transducer was set to 0.6 A for all experiments.
The comparisons between numerical results and experimental data for the main electrical output quantities in function of immersion depth h T are shown in Figs. 7, 9 and 10, for beakers A, B and C, respectively.
For the 1L beaker A, Fig. 7e shows that the electrical power increases with immersion depth, i.e. as the tip of the transducer approaches the bottom of the vessel. Our model captures this feature reasonably well, and this reflects the progressive enlargement of the cavitation zone on the lateral sides of the transducer as immersion increases. This is evidenced on Fig. 8 which shows the acoustic field and the cavitation zones for various immersion depths. This increase of dissipated power with immersion is a general trend in such experimental configurations [89], but other acoustic resonance effects may qualify this conclusion. For example, the peaks observed for some immersion depths in the power curve of beaker A on Fig. 7e are not captured by our simulations. The situation becomes more complex for larger volumes, as will be seen for beakers B and C. On the other hand, it can also be seen on Fig. 7e that here again, our computations overestimate the active electrical power by a factor of about two. This conclusion also holds for beakers B and C (Figs. 9d and 10d).
The experimental and computed values of the phases of transducer and motional branch impedances are superimposed on Figs. 7b, 9b and 10b for the three beakers. It can be seen that in all cases, the generator correctly maintains IZ BM = 0 despite the variations of the load as immersion depth varies (pink dashed line), which demonstrate that the control-loop for automatic frequency locking works properly. On the other hand, for beaker A (Fig. 7c), our model reasonably captures the variations of IZ Trans , and our estimation of the working frequency differs from the experimental one by less than 100Hz (Fig. 7a).
It is interesting to compare Figs. 7b and 7e, which show that the variations of the impedance module |Z Trans | closely follows the ones of the dissipated power, both experimentally and in computations. This is a natural result since the phase of Z Trans is close to zero so that |Z Trans | ≃ R(Z Trans ). Since the electrical power reads P elec = R(Z Trans )I 2 0,RMS , both curves are therefore proportional as long as the input current remains constant. Our overestimation of P elec is therefore also naturally reflected on |Z Trans |.
A practical instructive conclusion can be drawn from the power curves in function of immersion depth. The imposed current value corresponds in fact to a given position of the "power level" button on the generator. It can be seen therefore that fixing this position does not ensure a given active power since as shown on Figs. 7e, 9d and 10d the latter is immersion depth-dependent. The case of beaker B is particularly interesting since the dissipated power exhibits a marked maximum for an immersion depth slightly above 6 cm and a minimum near 10 cm. Thus, since the liquid volume remains constant, two experimenters working with the same graduation level, but with these different immersion depths, would conduct sonochemistry experiments with very different dissipated powers.
Going further into this issue, Figs. 9d and 10d show that our model also yields extrema in the power curves of beakers B and C, but not at the same locations as for the experimental ones. We suspect that the presence of these extrema reflects some acoustic resonance effects of the vessel for specific immersion depths. Were the liquid governed by linear acoustics, the location of these peaks would be independent of the driving amplitude I 0,RMS . But as the liquid behaviour is nonlinear, their location may vary with amplitude. As shown so far, our model overestimates the energy dissipated in the liquid, possibly predicting cavitation zones larger than the experimental ones, and therefore more non-linearity. The discrepancy between computed and experimental results observed for the locations of the power extrema might therefore be a further consequence of our overestimation of the dissipated power.
In order to check whether the correct peak locations can be recov-ered for smaller drivings, the current amplitude was lowered to I eff = 0.3 A in the simulations of beaker C, and compared to the experimental curve at I eff = 0.6 A (Fig. 11): the shape of the computed power curve now fits much better the experimental one, and exhibit extrema at the same immersion depths, but the electrical power is now underestimated. The above considerations clearly show that the vessel and liquid  geometry noticeably influence the energy transfer from the transducer to the liquid, which is a manifestation of acoustic effects. Our model captures such effects qualitatively, but is yet unable to make quantitatively correct predictions. Similar results have been reported in the literature [90][91][92][93][94][95]89,96] and deserve further investigations.

Discussion
The orders of magnitudes and trends of our predictions are in reasonable agreement with experimental results. Owing to the numerous multi-physics aspects of our simulations, which mimic as close as possible a sonoreactor experiment, these results are rather encouraging. But the observed systematic overestimation of the dissipated electrical power by a factor of about two deserves further discussion. As the latter directly reflects the mechanical power dissipated by cavitation, this reflects a general overestimation of the dissipated power by cavitation bubbles in the model of Ref. [59], despite the latter is based on real bubble dynamics. Apart from the semi-arbitrary choice of parameters R 0 and N, whose parametric variations have shown negligible influence, several enhancements should be considered.
One well-known weakness of our model is the arbitrary setting of the real part of k 2 to its classical linear expression (4), which is independent of pressure amplitude. This issue has been recently re-examined independently by Trujillo [97,64] and Sojahrood and co-workers [66,68]. They showed that both real and imaginary parts of k 2 can be expressed as period-averaged expressions involving the pressure field and the local void fraction β = N4πR 3 /3 (with possible extension to poly-disperse bubble sizes). The authors finally obtain the same expression for I(k 2 ) as Eq. (5) but R(k 2 ) is also found to depend on the local acoustic pressure |P|. This not only affects the local sound velocity but also the attenuation. Calculations for uncoated 2μm bubbles showed that the values of R(k 2 ) can drastically fall as acoustic pressure increases, in an increasing range of frequencies around resonance [68]. Interestingly, the authors found that the drop of R(k 2 ) could yield values of attenuation falling to near 50% of our own estimation, which is the order of magnitude of the discrepancy observed in the present work. Whether similar conclusions could still be drawn for the low f/f res ratio used in the present study remains to be clarified, and including the correct value of R(k 2 ) in our model is part of an ongoing work.
Our model still disregards the radiation losses [63,64,67,65], despite we use the Keller equation in our numerical computations. The historical reason for that is that the Caflisch model, from which our model is a simplification, was rigorously derived with the Rayleigh-Plesset equations. To what extent the latter can be replaced by a compressible bubble dynamics in the fully nonlinear model remains to be explored, despite it is a well-known practice in linear versions [55,[98][99][100][101]. Nevertheless, accounting for the radiative dissipation term has now become common, and it will be included in a future version of our model. One should note however that adding this term would further amplify the discrepancy observed in the present study.
Another possible overestimation of attenuation may be due to the fact that the bubble cloud is very dense in some regions (especially near   Fig. 8. Computed acoustic fields for beaker A (Fig. 7), for different immersion depths (from left to right: h T = 0.52,2.89, 5.26,7.63,10 cm). The color-plot represents the acoustic field amplitude in bar. The white lines are the computed paths of bubbles in regions above the Blake threshold, assuming they follow the primary Bjerknes force field [60]. The red line are the deformed shape of the transducer and beaker's walls. the transducer), invalidating the basic scaling assumptions of Caflisch model, which is valid to O(β). This should be solved by accounting for bubble-bubble interactions more carefully than the way they are modelled in O(β) models, where bubbles interact only through their individual coupling with the mean field. Some extensions have been explored for example by [98] under the linear assumption, and in a more general context by [44] which yield better agreement with the wellknown experimental results of Silberman [102] once linearized. Another approach may consist in solving coupled radial dynamics of a given set of N randomly dispersed bubbles [103,104] extending the twobubbles model of Ref. [105]. The latter approach also yields better agreement with Silberman's experiments near resonance [68]. We think however that such subtle refinements would probably compromise the simplicity of our model, whose efficiency is guaranteed by the ability to compute the power dissipated by a single bubble upstream from the COMSOL model (see A).
Bubble size poly-dispersity may also alter the dissipated power evaluation, despite the parametric study proposed in the present paper shows little sensitivity of the global power consumed to the precise value of bubbles ambient radius. The recent study of Lesnik and co-workers [72] shows however some interesting spatial segregation of bubbles depending on their sizes. Bubbles sizes are primarily bounded from above by bubble breakup, and owing to surface tension, small bubbles are expected to be more resilient against large acoustic pressures. Defining a main bubble size dependent on the local acoustic pressure in our model is under investigation.
Finally, we would like to stress another flaw in our model, which is related to the location we choose for those inertial bubbles that dissipate energy, as defined by Eqs (2). The latter impose the presence of bubbles in any zone where the acoustic pressure lies beyond the Blake threshold. However, it is well-known that in a standing wave, the primary Bjerknes force can become repulsive above some threshold, owing to the increasingly long expansion phase of inertial bubbles with the driving pressure [106,107,62,108,86,60]. This threshold is close to 1.75bar in water, and this feature is known to yield the so-called "Jellyfish structures" [109][110][111]62], where a large pressure antinodes in a standing wave zone repels the bubbles toward the threshold location. Eqs. (2) do not catch this subtlety and blindly puts some bubbles in such regions, which prohibits the prediction of such bubble structures in our simulations. Our model thus predicts streamers that are not present in experiments, and erroneously compute some dissipation there. In an ongoing work extending the present one, we performed comparisons between luminol maps and the model's predictions, that would tend to support this interpretation. The shift in the predicted locations of the extrema on Figs. 9d and 10d might also stem from the same issue. Whereas the other above-mentioned imperfections of our model lied in the potential inaccuracy in the evaluation of the bubble dissipated power, the present one is due to the incorrect prediction of the spatial location of bubbles. Correcting this problem rigorously would thus require modelling the translational motion of bubbles as in Ref. [72]. However, in order to keep reasonable computation costs, we plan to use an in-between solution in order to patch the condition used in Eqs. (2), taking advantage of the ability of our model to compute the primary Bjerknes force field [60].

Conclusions
Louisnard's simplified model of nonlinear sound propagation in cavitating liquids [59] has been included in a global sonoreactor model, accounting for the vessel walls vibrations, the transducer, and the automatic frequency control strategy of the generator. For a given geometrical configuration, the only required parameter of the simulation is the amplitude of the input current, which is the quantity imposed by the power button of the generator. Therefore, our simulations exactly mimic experiments, without resorting to debatable arbitrary boundary conditions. Our model provides estimates of experimentally measurable quantities such as working frequency, electrical complex impedance, and consumed electrical power. To the best of our knowledge, this is the first reported model allowing such fine and comprehensive simulations of acoustic cavitation experiments, allowing large parametric studies, and directly testable against experimental data.
Experiments were also performed, varying either the input current or the transducer immersion depth, for three different beaker geometries. Our global model was put to test against the electrical data provided by the ultrasound generator in the latter experiments. Despite the correct order of magnitudes for frequency, impedance and electrical power are caught by our simulations, there emerges a net tendency of our model to overestimate the consumed electrical power, by a factor of about two. Parametric variations of the bubble ambient radius and bubble density in our model did not fix the discrepancy observed.
Geometrical effects can be clearly observed both in experiments and simulation, and materialize as extrema in the consumed electrical power when the transducer immersion depth is varied, but some discrepancy is observed on the location of these extrema. This suggests that strong variations of a sonochemical yield may be observed when changing the transducer immersion depth at constant input current and constant liquid volume. Understanding and correctly predicting these effects to optimize the shape of sonoreactors is a future challenge for extensions of

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Olivier LOUISNARD reports financial support was provided by Sinap-Tec. Laurie BARTHE reports financial support was provided by Sinap-Tec. Igor GARCIA-VARGAS reports a relationship with SinapTec that includes: employment and funding grants. the cavitation research group (GDR CAVITATION) that stimulated fruitful discussions with all our colleagues involved in the field. We also thank Benjamin Laulier, Sébastien Barbry from SinapTec for clarifications on the ultrasound generator and Samuel Clere for his help in the design, manufacture and characterization of the transducer.

Appendix A. Bubble dynamics model
The dynamics of the bubbles is computed by a Keller equation, accounting to the liquid compressibility up to first order in Mach number Ṙ /c l [112,113]: with R(t) is the bubble radius, c l , ρ l and μ l are the sound velocity, density and viscosity of the liquid, respectively, p g is the pressure in the bubble, σ the surface tension, and is the sinusoidal driving pressure at frequency f, around the ambient pressure p 0 . We take this opportunity to report an error in our original paper [59]: we reported there that we simulated the Rayleigh equation in order to obtain the dissipated power Π v and Π th . In fact the Keller equation was used and this is true in all our subsequent work, including the present one.
The pressure in the gas is assumed uniform. The gas content is a mixture of several species (here H 2 O and air), with number of moles n i , represented by an approximate van der Waals equation of state: where R is the universal gas constant. Following [114], the reduced volume is computed as averages of the individual species reduced volumes: To account for the bubble thermal behaviour and heat transfer with the liquid we followed the following simplified approach: the bubble internal temperature is assumed uniform except in a near-wall layer of thickness δ th [115][116][117]. The latter is estimated by: where Pe th = RṘ/α g is a thermal Peclet number, α g is the thermal diffusivity of the gas mixture in the bubble. The cutoff 1/π proposed in Ref. [115] is needed to avoid nonphysical large values of δ th during nearly isothermal phases of the bubble oscillation. It was shown by Preston et al. [117] that the correct value of the cutoff is rather 1/5. Following this approach, the temperature of the bubble core can be obtained by applying the first law of thermodynamics to the whole bubble content: In the latter equation, the left-hand side is the time-derivative of the bubble gas mixture internal energy the first term in the right-hand side of (A.6) is the rate of work done by the liquid on the bubble, the second is the enthalpy flux of water evaporating from (or condensing into) the liquid with a molar flux ṅ H2O , and the last one is the conduction heat transfer between the bubble and the liquid, assumed at constant undisturbed temperature from ambient [115].
To account for water evaporation/ condensation at the bubble wall, an approach similar to heat transfer is used. Water vapour concentration is assumed constant in the bubble except in a near-wall layer of thickness The C v (T) and C p (T) are computed from Ref. [119]. Following [120], the gas mixture conductivity is computed using method of Chung et al. see [121] Sections 10.3 and 9.5. This may allow estimations of λ for high density states, but due to convergence issues, we followed Ref. [115] and estimated λ at liquid temperature. We also followed Storey & Szeri [120] to compute the diffusion coefficient between air and water, with the Brokaw correction to account for polar molecules (here water) see [121] Eq. (11.3.7) and followings. Despite the correction proposed in appendix of Ref. [120] to account for high densities is available in our code, the diffusion coefficient was also estimated at liquid temperature.
The set of differential equations (A.1), (A.6) and (A.9) is solved with a FORTRAN code using the stiff solver DDRIV3 from NIST Core math library (CMLIB), over 200 acoustic periods. The solution yields the unknowns R(t), T(t) and n H2O (t), from which the gas pressure can be obtained from Eq. (A.3). The power dissipated by the bubble over the last cycle can then be estimated with the integrals: The graphs reported in our original paper [59] were obtained with the above model and Π th was found to be much smaller than Π v in the parametric range of interest, which justifies that only the latter is retained in Eq. (5). It should be noted that a different conclusion was found by Jamshidi & Brenner [63,72]. This might be due to the fact that the latter authors cut the computation of the integral at the end of the collapse, rather than over a whole bubble cycle. Contrarily to the integrand in the viscous term in (A.12) which is always positive, there are sign changes in the one of (A.11), as the heat flow changes direction along the oscillations (even if the integral over one cycle must be positive). The latter authors would miss therefore some post-collapse events involving an inward heat flow occurring during part of the rebounds, associated with very low transient temperatures attained in the bubble, consecutive to rapid adiabatic re-expansions [122]. Other reasons for discrepancy may lie in the tricky computation of thermal properties of the gas mixture. The bubble response is computed by parametrically varying the dimensionless driving amplitude P * by small steps between 0 and 3, to compute the viscous contribution Π v , for a given ambient radius R 0 and frequency f. The resulting curve can be perfectly fitted in dimensionless form by: , (A. 13) and is injected in COMSOL as an explicit function to compute I(k 2 ) from (5). We emphasize that the whole procedure must be repeated for any change of R 0 ,f, as well as when changing the liquid or the bubble gas content, or the ambient conditions. A universal analytical expression for Π v avoiding this straightforward but lengthy procedure is currently under study.