Analysis of dynamic acoustic resonance effects in a sonicated gas–liquid flow microreactor

Highlights • Dynamic acoustic resonance in gas–liquid microreactors is modelled.• Gas-liquid interface deformation and atomization are accurately predicted.• Faraday instability at resonance peaks leads to transient atomization bursts.


Introduction
The integration of ultrasound actuation with micro-scale flow reactors is often chosen as a route to further intensify processes [1][2][3][4], by e. g. improving mixing [5][6][7], and in case of two-phase flow interfacial mass transfer [8,9]. In general, these ultrasonic microreactors are either operated in the low frequency (below 200 kHz) or high frequency (above 1 MHz) range. Low frequency ultrasound is associated with the formation of cavitation bubbles in the microchannels, which are then exploited to intensify mixing and interfacial mass transfer, but also for e. g. synergistic effects in particle synthesis [10][11][12][13][14]. Contrary to that, high frequency ultrasound is not dominated by cavitation effects, as the applied power levels do not exceed the cavitation threshold. However, the wavelength of the applied ultrasound in the fluid matches the channel size, making it possible to form a standing wave within the channel and utilize the associated effects, such as the acoustic radiation force and streaming [15,16]. These effects are commonly exploited to manipulate particles [17][18][19][20][21][22], and to enhance mixing [23][24][25][26].
As highlighted by the examples above, the application of low and high frequency ultrasound in microreactors is quite advanced and has found various applications in biological, pharmaceutical and chemical processes. However, the transition frequency range (between 200 kHz and 1 MHz) is largely unexplored in microreactors. Recently, Mc Carogher et al. [27] reported the observation that gas-liquid segmented flow can act as a set of acoustic resonators in a frequency range from 100 kHz to 500 kHz. At resonance, a large acoustic amplitude is reached, leading to gas-liquid interface deformation, atomization of micrometer sized droplets, and cavitation. This resonance phenomenon in the transition frequency range could be exploited to e.g. increase interfacial area between the gas and liquid phase, and intensifying reactions that are limited either by mass transfer and/or phase distribution. However, so far no model to describe the acoustic phenomena and the acoustic behavior of the two phases in the microchannel is available. The development of such a model is paramount to increase the fundamental understanding of the resonance phenomenon, and to further design and optimize microreactors based on it.
In this work, a numerical modeling framework is presented which predicts the acoustic field in the microchannel, the deformation of the gas-liquid interface, and the atomization threshold. This approach focuses on modeling dynamic acoustic resonance, which is the main mechanism that induces the observed phenomenon. The next section briefly reviews the experiments used for validation, followed by the presentation of the model and its underlying assumptions. Then, the comparison between numerical results and experiments for different relevant configurations is presented, followed by an in-depth discussion and concluding remarks.

Description of the experimental setup
Experiments were performed in two microreactor setups: (i) a silicone etched microreactor with 1.2 × 0.6 mm cross section square channel and (ii) a single glass channel with a 1.2 × 1.2 mm or 1.0 × 1.0 mm square cross section. Both microreactor setups are integrated with ultrasound by gluing a piezoelectric plate to the bottom surface, and the selected actuators were: (i) Pz26, 80 × 40 × 4.06 mm, (Ferroperm TM , Denmark), (ii) Pz26, 20 × 11 × 1.67 mm and 80 × 20 × 1.67 mm (Ferroperm TM , Denmark). A detailed description of the setups is given in Mc Carogher et al. [27]. Various gas-liquid flow rate ratios were delivered to the reactors to achieve gas-liquid segmented flow, while being sonicated at powers from 5 W to 20 W, and specific frequencies associated with actuator resonance from 40 to 600 kHz. The flow and the  resulting acoustic resonance phenomena were recorded with a microscope (SMZ25, Nikon) equipped with a high-speed camera (Fastcam mini Ux100, Photron limited). To validate the model, the results for an actuation frequency of around 500 kHz were used, as the most significant phenomena were observed at this frequency. Long air bubbles (with a length of 2 to 6 times the channel height, typically used in practice) were considered. The liquid and gaseous stream are generated by equipment chosen to deliver a steady continuous volumetric stream without disruptions (Fusion 4000x Syringe pump, Chemyx Inc. TM , Stafford, Texas and EL-FLOW mass flow controller, Bronkhorst TM , Vorden NL). These steady streams are combined in a microfluidic junction built inside the chip. The steady conditions and microfluidic combination led to a repeatable and predictable process, so that each bubble was generated with the same volume after a defined time span. The relevant fluid properties and geometric parameters selected for this study are listed in Table 1.

Classification of phenomena
From the experimental observations in Mc Carogher et al. [27], it is possible to classify the different dynamic behaviors driven by acoustic resonance as follows (see also Fig. 1 for illustration): • Steady flattening of the gas bubble • Strong transient shock on the gas bubble • Nonlinear oscillations of the gas bubble • Atomization at the sonicated surface of the gas bubble In general, interface flattening occurs prior to atomization. As illustrated in Figs. 1 (c.1) and (c.2), flattening can be observed with or without atomization at the flattened interface. All of these observations highlight resonance order in the channel axis. The cutoff frequency for the existence of high order resonance in the channel cross section is f c = c 0 /(2h) = 617 kHz, for a 1.2 × 1.2 mm channel, and 740 kHz for a 1.0 × 1.0 mm channel. The selected frequency was thus below this first cutoff frequency at which a pressure node would appear at the mid-channel height.
It is possible to define a mode using the notation (n,m), where n ∈ {1, 2, …} stands for the order along the channel axis and m ∈ {0, 1, 2, …} stands for the order along the two axes of the square cross section. Since the gas-liquid interface may be approximated as a free surface, with a zero acoustic pressure boundary condition, the liquid acts as an open acoustic resonator. Examples of mode shapes of interest are given in Appendix A in Figs. A.26 and A.27.

Nonlinear acoustics modeling
In linear acoustics, the average motion of a particle oscillating in an acoustic field is null. The non-zero average motion is a second-order nonlinear effect and is described as the effect of the acoustic radiation force 1 . According to Baudoin and Thomas in their review on particle and fluid micromanipulation, a definition of the acoustic radiation force would be: a net force applied at the interface between two media with different acoustic properties, this force enables manipulations of particles but also deformation of fluid interfaces [28]. Since a standing acoustic wave is determined by boundaries, which act as reflectors, moving boundaries result in a dynamic acoustic field. As a consequence, the considered configuration involves an additional type of nonlinearity, which is the deformation of the domain where the acoustic field is set. The acoustic radiation force is directly determined by the acoustic field via its amplitude and pattern (nodes and antinodes location, focusing). A large acoustic amplitude occurs when resonance conditions are met. Consequently, when the system changes over time, resonance conditions vary and can be met at specific instants. This can be summarized under the concept of dynamic acoustic resonance.
Dynamic acoustic resonance effects have been addressed by Simon et al. for the acoustic fountain configuration [29,30]. Other authors have studied similar configurations [31]. The focusing of a standing ultrasonic field, generated by a focused ultrasound transducer, on a free liquid surface induces a rise of the liquid. When the height equals a specific fraction of the wavelength, a sudden local resonance produces droplet formation. The sequence of resonances can give rise to a droplet chain growing vertically. We noticed that the principle involved in Simon's configuration was similar to the one encountered in our microreactor. Contrary to Simon's configuration, in the present study, the high acoustic amplitude is not obtained by focusing but by pure resonance. In such a configuration, the acoustic radiation force effects are much more significant than acoustic streaming effects.

Derivation of the acoustic radiation force acting on a free surface
The explicit formula for the acoustic radiation force acting on a gas bubble boundary can be derived from the perturbation expansion method. This method is commonly used for acoustic streaming and particle manipulation prediction, see e.g. [32][33][34][35][36]. Let us consider a liquid domain Ω, whose boundary is composed of walls Γ walls , and free surface Γ free . Writing the expansions of pressure, velocity and density, for v 0 = 0 (fluid at rest relative to the reference frame) and p 0 = 0 (relative to the static pressure, comprising equilibrium surface tension  pressure and fluid pressure). The speed of sound in the liquid is denoted by c 0 . The boundary layer is neglected, and a slip condition is set at Γ walls .
Since one can show that the capillary wave and acoustic wave solutions are uncoupled, assuming a pure acoustic wave in Ω, the boundary condition at Γ free can be approximated simply by p = 0. The time-averaged second order equation assuming an incompressible and inviscid fluid and the pure harmonic solution of the first order equations defined in Ω are where 〈⋅〉 T denotes time averaging over a period T. Since a change of the free boundary geometry from the equilibrium geometry is expected, a resulting capillary force emerges. The boundary condition at Γ free for the second-order system is p = − γΔ s η 2 , where η 2 represents the normal displacement at the surface (pointing outwards of Ω), γ is the surface tension, and Δ s denotes the surface Laplacian. Finally, assuming no additional forces other than surface tension, one obtains the following synthetic equation at Γ free : The acoustic radiation force acting on a free surface can be also derived from Zakharov and Filonenko's modeling of capillary wave turbulence at the surface of a semi-infinite fluid [37,38]. The modeling involves the normal displacement of the surface η, the scalar velocity potential (irrotational) ψ, defined at the surface and in the fluid. Considering a two-dimensional domain (thus a one-dimensional surface) and applying the perturbation method the linear capillary wave model is: where x stands for the space variable along the one-dimensional surface. ∂ x denotes the one-dimensional 1/2 fractional Laplacian which can be defined, for an infinite surface, by the Fourier transform: with F being the Fourier transform 2 of the function f. Eqs. (7) and (8) describe linear capillary waves whose dispersion relation is ρ 0 ω 2 = γk 3 .
Applying the expansion method while assuming a purely harmonic solution of the first order system, the second order time averaged equations of Zakharov and Filonenko's model are: Eqs. (10) and (11) can be combined to obtain the following synthetic equation: Eqs. (6) and (12) provide two different frameworks to desribe the acoustic deformation of surfaces, which consist of a balance between inertia, surface tension and the acoustic radiation force. Both modeling approaches give the same result for the case of the acoustic deformation of an infinitely long liquid film of a given initial depth, as shown in Appendix B. This case has been studied by Scortesse et al. [39], where a similar expression for the acoustic radiation force in a liquid film was derived.

Numerical scheme involving an updated standing acoustic field
The selected numerical scheme for dealing with dynamic acoustic resonance in two-phase media is summarized in Fig. 2. The main challenges of the definition of the numerical method are: (i) a relatively strong nonlinear behaviour is expected, (ii) dealing with a large acoustic amplitude range, (iii) update of the computational domain (i.e. remeshing). The main idea of our approach was to combine a frequencydomain acoustic solver with a time-domain surface Lagrangian solver for the deformation of the gas bubble. This approach is valid if the establishment time of the standing acoustic wave is much smaller than the characteristic deformation time of the surface. Since the typical maximum velocity of the observed moving interfaces is in the range of 2 m/s, with the speed of sound in water being around 1480 m/s, one can propose that this is a reasonable approximation. The solvers are described in detail in the following paragraphs.

The 2D approximation
A two-dimensional approximation was chosen to model the gas bubble, which represents a slice of the actual three-dimensional gas bubble, as illustrated in Fig. 3. The gas bubble was thus modeled as a closed curve in the x, y plane. This choice of a 2D model has two motivations. Firstly, the relative simplicity of the acoustic field does not imply significant 3D motion, in the sense that the motion can be reduced in the channel axis and in one of the cross section plane dimensions of Typical configuration for the two-dimensional finite element acoustic solver. In this example, excitation is set at the bottom.
the channel due to symmetry. Secondly, the 2D model allowed to perform a large set of simulations with fine mesh both in time and space, which aided the complete analysis.

The frequency-domain acoustic solver
The linear acoustic equations are reduced to the frequency-domain pressure formulation: where ω is the pulsation of the monochromatic acoustic wave (in rad/s).
The boundary conditions are shown in Fig. 4, which were: p = 0 at Γ free , ∇p⋅n = 0 at Γ fwall and Γ symmetry , and ∇p⋅n = ρ 0 ω 2 u a at Γ excitation , where u a is the oscillation amplitude of the given excitation. Due to the significant acoustic mismatch between air and water, the acoustic field transmitted inside the gas bubble can be neglected. A common finite element method (FEM) was implemented. For a given gas bubble geometry the domain Ω was meshed with triangular elements. The size of the regular mesh elements was 0.6e − 4 m, and the size of the fine mesh elements at the bubble interface was 1.5e − 5 m (the mesh is illustrated in Fig. 4). Linear elements with three integration points where used. The computed acoustic field is then assumed to be in the form p 1 (x, y, t) = p(x, y)sin(ωt), where p(x, y) is the solution of the frequency domain problem. Time averaging, when computing the acoustic radiation force, introduces a 1/2 factor from 〈sin 2 (ωt)〉 T or 〈cos 2 (ωt)〉 T . The particular oscillation velocity v 1 and amplitude u 1 were computed according to: The resonance consists of a singularity in the sense that the acoustic field is not defined at the resonance frequency. In practice, the acoustic field at resonance has a finite value because of physical saturation. To tackle the singularity, a loss factor of 5e − 4, set empirically, was chosen for artificial dampening in the following results.

Lagrangian surface evolver approach
A Lagrangian "Surface Evolver" was implemented in MATLAB (version R2021a) to capture the gas bubble deformation, which is a similar method as the Surface Evolver by Kenneth Brakke [40,41]. The method consists of considering the motion of the surface driven by surface tension, external forces and constraints, without computing the flow field in the surrounding fluids. The prediction of the motion of the liquid in the bulk domain was not of essential interest compared to the motion of the interface. One of the benefits of such an approach is that it makes the computation straightforward, since the number of nodes is considerably reduced, which allowed us to perform a large number of simulations.
The closed curve associated with the gas bubble geometry was meshed by evenly spaced nodes. The cell size was 1.5e − 5 m. Volume conservation was ensured by adding a constraint to the numerical area enclosed by the curve at each time step, which then also approximates the incompressibility condition. As a consequence of the Surface Evolver approach, Zakharov's formulation in Eq. (12) was used, which is valid for closed surfaces. The fractional Laplacian in Eq. (9) is then transposed to its periodic domain equivalent involving a Fourier series form. A central difference scheme and an explicit strategy was used for the acoustic and capillary forces, so that the latter are computed according to the current time step only. For each node position x = (x(t), y(t)) of the surface the Lagrangian integration had the following form: where the subscripts denote the time step indexes, Δt the time step, and n t represents the outward liquid surface normals at the current time step. The capillary force was computed by The curvature at the current time step κ t is given explicitly by

Fig. 5.
Illustration of atomization generation stages at a gas bubble interface exposed to a resonant standing wave in a microreactor. Artificial numerical contrast was applied to the images. The images are taken from a 32,000 fps video highlighting a small part of the interface, see Video 1 in the ESI.
where s denotes the local coordinate of the curve. The computation of the fractional Laplacian in Eq. (17) was done using the Fast Fourier Transform (fft) equivalent of Eq. (9). For a given current solution of the frequency-domain acoustic solver, the acoustic force f acoustic,t was computed according to the third term in Eq. (12).
The channel walls were modeled as a force preventing penetration of the surface through them. The setting of the force constant of the walls was done empirically according to the equilibrium shape problem of a gas bubble in a channel (one result of the equilibrium shape problem is provided in the Supplementary Material, see reference in Appendix D).
When dealing with strong deformation mesh distortion arises, hence a remeshing strategy is necessary to ensure mesh regularity. The main issue with remeshing of a surface during motion is the error and possible aberrations coming from the interpolation of the speed from the current mesh to the next one. This error can be reduced by reducing the mesh size. We noticed that the interpolation error is quite dependent on the remeshing partitioning strategy. One way to define partitioning is to define Lagrangian trackers that are not affected by the remeshing and to define the partition of the surface according to them. While computing the motion, one needs to update the acoustic force according to the new boundary. In this study, updating time steps from 5e − 6 s (fine) to 1e − 4 s (coarse) were chosen.
To validate the selected solver, comparisons with OpenFOAM volume of fluid (VOF) simulations were performed for the following initial shape problems: (i) merging of a droplet on a flat surface, (ii) bump on a flat surface. Another comparison was also done for the case of the rising droplet based on Hua and Lou's benchmark [42]. Videos of these preliminary simulations are provided in the Supplementary Material, see Appendix D. In general, the behavior of the interface is well reproduced and the surface solver could deal with strongly concave shapes in a stable way.

Methodology for atomization prediction
Atomization consists of the breakup of bulk liquid in multiple ejected droplets at a free surface. As observed by Mc Carogher et al. [27], different mechanisms of atomization are involved in a gas-liquid microreactor. These can be divided into three types: (i) capillary, (ii) cavitation, and (iii) coupled capillary-cavitation. The focus of this work is on the Faraday instability, which belongs to the capillary type, as one of the mechanisms resulting in atomization. In the two-phase microreactor, the (0, n)-mode standing wave results in an oscillation of the bubble interface in the channel axis direction at the ultrasonic excitation frequency. This oscillation, when reaching a large amplitude, leads to the Faraday instability and, at some point, to atomization. Based on Puthenveettil and Hopfinger's analysis of the oscillating tank experiment [43], one can divide the capillary atomization generation mechanism into 4 stages associated with different excitation amplitude ranges: • Linear regime • Standing subharmonic nonlinear Faraday wave • Chaotic motion below the atomization threshold • Atomization The same stages are observed in the ultrasonic microreactor, as illustrated in Fig. 5 and in Video 1 (ESI).
In the linear regime, a standing linear capillary wave oscillating at the driving frequency can be excited at the surface. In this regime, the superposition of waves is valid and simple patterns, such as a plane wave or interference patterns, can be observed. Depending on the configuration, the linear capillary wave can be caused by the presence of a contact point with a vibrating solid surface, as in the vibrating tank experiment. Above a threshold, a standing subharmonic Faraday wave is observed. Faraday crystals can form at the gas bubble surface, where mostly the square pattern is observed, as shown in Fig. 5 (c). The chaotic stage is characterized by an unstructured surface deformation with salient crests. Finally, atomization happens when breakup occurs at the crests, and at the onset of breakup ligament or elongated droplets can be observed. Fig. 6 provides an overview of the different standing wave patterns observed in the experiments.
The theoretical description of the Faraday instability was developed by Kumar and Tukerman [44]. The Kumar and Tuckerman theory (K&T theory) is based on analytical solutions of the Navier-Stokes equations for the interface between two fluids (of finite or infinite depth) oscillating at a given driving frequency. The K&T theory is known to be an accurate method to predict the instability threshold. It is possible to define the Faraday instability threshold with respect to the critical point of the first subharmonic instability region. This corresponds physically to the emergence of a standing wave oscillating at ω/2, where ω is the driving pulsation (in rad/s). The critical point is determined by the dynamic viscosity of the two fluids and their interfacial tension. In case of an air-water interface the effect of the air viscosity is negligible. The instability chart for an air-water interface with a driving oscillation frequency of 500 kHz is depicted in Fig. 7. The oscillation is assumed to be a pure harmonic forced displacement of the interface, of the form u =  Two-dimensional volume of fluid (VOF) simulations were performed to aid the atomization prediction. For this, the standard OpenFOAM  solver interFoam was modified to include a smoothing method to suppress spurious currents (see [45] for solver details). The geometry consisted of a rectangular 2D tank with a driven vibrating bottom (oscillating at 500 kHz for different oscillation amplitudes, which were: 0.4 μm, 0.5 μm, 1.0 μm, 1.5 μm, 2.0 μm, and 2.5 μm) and periodic boundary conditions at the sides. To enable the vibration of the bottom wall, a dynamic mesh was applied using the displacementLaplacian solver, which solves cell motion based on the Laplacian of the diffusivity and the cell displacement, see [46] for details. The schematic of the configuration with dimensions is shown in Fig. 8 (a), and the fluid properties are given in Table 1. This 2D approximation allowed to perform the different simulations associated with different excitation amplitudes within a reasonable computation time. The mesh was generated using blockMesh and then refineMesh was implemented to refine the mesh near the interface (for a thickness of 5e − 5 m), see Fig. 8 (b). The final mesh consisted of 50360 cells with a cell size of 4.7e − 6 m in the regular region, and 1.3e − 6 m in the refined region. The time step size was automatically changed based on a maximum Courant number of 0.2. These VOF simulations enabled the prediction of the air-water interface behavior oscillating at 500 kHz for the different oscillation amplitudes specified above, and their position against theoretical instability tongues are indicated in Fig. 7.

Natural deformation patterns
First, a qualitative view of the natural deformation patterns associated with particular modes is given in Fig. 9, which highlights the dependence of the deformation of the sonicated interface on the resonance mode. The (0, 1) and (0, 2) modes (not displayed) have a similar deformation pattern as the (0, 3) mode. This particular deformation pattern consists of a flattening of the interface. The middle of the   interface is pushed towards the inside of the gas bubble while the sides are pulled towards the walls. For the (1, 1) mode, the middle of the bubble is pulled towards the liquid, leading to a bump instead of flattening. The number of bumps-concavities increases with the order of the mode in the channel cross section. A bump towards the liquid is associated with a pressure node near the interface, while a concavity is associated with a pressure antinode near the interface.

One-side transient flattening on free bubble
First, the simulation results for the case of a gas bubble without flow constraint submitted to a 500 kHz one-side resonance are discussed. The excitation was located on the right boundary to excite the (0, 3) mode, and the simulation parameters for 5 different excitation magnitudes are provided in C in Table C.6. Video 8 illustrates the obtained behavior for the magnitude 4 case, and additional videos in the Supplementary Material show the evolution of gas bubble and acoustic field magnitude over time, see reference in Appendix D. The evolution of the acoustic magnitude, given in terms of maximum acoustic pressure and maximum particular oscillation amplitude, and the mid-channel point (y = 0) velocity, are depicted in Fig. 10. For each magnitude, a resonance peak emerges resulting in interface flattening near the ideal rectangle resonator length for (0,3) mode, which in this case is at about x = 2.4 mm. The velocity of the mid-channel point quantifies the speed of this deformation, and the obtained values for each simulation are provided in Table 2. The higher magnitudes being pushed faster to the resonance position, it can be seen that the associated resonance peak is reached sooner than for lower magnitudes. The flattening shapes are compared to experimental results corresponding to Video 4 in Fig. 11. The measured flattening velocity for this experiment was about 0.12 m/s. The flattening magnitude increases with the acoustic field magnitude, from weak deformation to a strongly flattened shape. An acoustic pressure peak in the range of 10 MPa is needed to produce a visible transient flattening, which requires a velocity in the range of 0.2 m/s. A maximum flattening magnitude is obtained for this configuration, see associated shape in Fig. 11 (b). This limit is reached for a large acoustic magnitude exceeding 80 MPa. For these elevated pressure values a high dynamic range is obtained in a short time of about 1e − 5 s. After the flattening, the gas bubble, in this particular free condition, moves away from the resonance position due to inertia.

One-side transient flattening on constrained bubble
To take into account the presence of the surrounding liquid in the channel, a volume constraint on both sides of the gas bubble was added to the previous configuration. Accordingly, the gas bubble is unable to move away from the resonance position. Seven excitation magnitudes were chosen, see Table C.7. Video 9 illustrates the obtained behavior for the magnitude 4 case, and additional videos in the Supplementary Material show the evolution of the gas bubble and acoustic field magnitude over time, see reference in Appendix D. The plots of the acoustic magnitude are shown in Fig. 12. The introduction of the volume constraint leads to a different behavior compared to the free bubble case. After a resonance peak, a steady state is reached. The acoustic field stabilizes at around |p 1 | = 2 MPa, |u 1 | = 0.55 μm for magnitudes 3 to 7, see Table 3. For the selected simulation time, no peak is observed for magnitudes 1 and 2, as the transient deformation towards the resonance position is too low to reach a peak.

Both-sides resonance
In this section, simulation results for the case of both-sides resonance are presented. Unlike the previous 2 configurations, the excitation is  now applied at the bottom boundary, for which four excitation magnitudes were selected, see Table C.8. Consequently, this configuration does not exhibit symmetry along the mid-channel axis anymore, and acoustic standing waves appear on both sides of the centered gas bubble. Video 10 illustrates the obtained behavior for the magnitude 3 case, and additional videos are provided in the Supplementary Material, see reference in Appendix D. As the initial mesh is not perfectly symmetric, a delay is observed between the left and right resonance peaks, see for example the magnitude 2 curve in Fig. 13. After the two resonance peaks, a steady state is reached. The comparison of steady flattening of the right side for the different magnitudes is shown in Fig. 14, and the acoustic field and deformation values are summarized in Table 4.

One-side resonance on an advected gas bubble
To deal with the case of a gas bubble advected towards a resonance position, the following configuration was considered: A flow velocity of 0.02 m/s was applied to the gas bubble starting from a position about 0.5 mm away from the resonance position, with the excitation located on the right boundary at four different magnitudes, see Table C. 9. The evolution of the gas bubble shape for the case of magnitude 4 is illustrated in Fig. 15 and visualized in Video 11, additional videos are provided in the Supplementary Material, see reference in Appendix D. Figs. 16 (a) and 17 (a) depict the simulation results at the start of acoustic resonance, from 14 ms to 35 ms, and Figs. 16 (b) and 17 (b) plot the simulation results at the end of acoustic resonance, from 49 ms to 57 ms. Initially, an increase of the acoustic amplitude is observed, reaching values which are independent of the excitation magnitude and where flow and the acoustic force are balanced. For magnitude 1 the acoustic force is not strong enough to maintain this balanced regime, and the gas bubble returns to the advected equilibrium shape. For all other magnitudes resonance peaks are observed after the almost steady regime, i.e. one peak is observed for magnitude 2, two for magnitude 3, and three for magnitude 4. Then the bubble returns to the advected state.

Resonance shock on one side
The obtained results highlight that an acoustic resonance peak can lead to a strong transient deformation of the gas-liquid interface when the acoustic force reaches a large magnitude for a short time. However, the value predicted by the acoustic solver may not be exact due to the singular nature of resonance. Therefore, it may be more rigorous to quantify acoustic resonance in terms of energy density of the associated peak, which would represent the product of the acoustic field amplitude by its extent in time. Using this description, a trivial resonance peak in the time-domain can be modeled as a rectangular pulse of period T and amplitude A, so that the energy density of the peak would be T × A. More generally, for a given time period T, it is possible to use the following mathematical definitions for the acoustic force energy density E f : and for the pressure energy density E p : where f acoustic and p 1 are the maximum acoustic force and the acoustic pressure, respectively, computed by the acoustic solver. Simulation result are given for a 2.6 mm long gas bubble subjected to a strong resonance shock from equilibrium shape. The peak energy density values were selected according to the strongest deformations that were   Fig. 15. Evolution of the gas bubble shape for the case of the advected gas bubble exposed to one-side resonance for the case of magnitude 4.
observed in practice. The values are given in Table C.10, and associated videos are available in the Supplementary Material, see reference in Appendix D. The single shock solution is also illustrated in Video 12 corresponding to the magnitude 3 case. The excitation consisted of a pulse of 1e − 4 s duration, which corresponds to 1000 integration time steps and to 1 frame for a 10,000 fps video. Fig. 18 depicts the qualitative comparison of maximum shape deformations for different resonance peak magnitudes with experimental observations (Video 5). The evolution of both position and velocity of the mid-channel point for the five magnitudes is shown in Fig. 19. The shock produces a concave shape, whose extent depends on the magnitude. Traveling capillary waves are also observed along the interface. When the gas bubble recovers its initial shape after the acoustic resonance shock, a new shock can arise, as shown in Video 6. According to this footage, standing Faraday waves emerge for about three frames, which corresponds to about 3e − 4 s. In the following results, the excitation consisted of a sequence of pulses, equal in magnitude, of length 3e − 4 s (sequence of square inputs). The pulse starts to be active when resonance is met, that is to say when a resonance peak is detected numerically. Simulations for different magnitudes of the resonance shock are performed, see Table C.11. Video 13 shows the periodic solution corresponding to the magnitude 3 case, and additional videos are available in the Supplementary Material, see reference in Appendix D. The mid-channel point oscillation is measured and plotted in Fig. 20  (a). After a transient stage of about 4 ms, the deformation of the interface follows a steady cycle. As for the previously discussed single shock results, velocities towards the inside (acoustic force pushing the interface) and the outside (surface tension recovering the shape) of the gas bubble are different. The measurement of the mid-channel point position in the experiment suffers from inaccuracies due to the difficulty to track the point in the footage. Consequently, the values of velocity are approximate. The Fourier spectrum of the oscillation for the different magnitudes is given in Fig. 20 (b). The measured oscillation frequency for Video 6 is about 1000 Hz, giving a period of 1 ms. Numerical results show that the oscillation frequency depends on the magnitude of the peak. Globally, one can observe that a stronger shock produces a more significant deformation that takes longer to recover, extending the oscillation period. It is expected that the recovery period depends on the static radius of curvature of the gas bubble, as the capillary pressure increases with decreasing radii of curvature. To confirm this, similar simulations were performed for two additional radii of curvature R of the gas bubble, varying slightly the channel height, for the case of   Table 5.

Faraday instability VOF results
VOF simulations of the Faraday instability are presented in this section. The phase fraction (extracted at α = 0.5 ) at maximum deformation of the free surface for the six selected magnitudes is shown in Fig. 21. The complete videos of the simulated interface evolution are provided in the Supplementary Material, see reference in Appendix D.
Small amplitude standing waves are observed for u = 0.4 μm and u = 0.5 μm. Typical crest-shaped standing wave starts to appear at u = 1.0 μm, which become higher at 1.5 μm. Chaotic motion and breaking are observed at u = 2.0 μm and u = 2.5 μm.   Fig. 23 (b). At the beginning of the simulation, the free surface oscillates homogeneously at the excitation amplitude. Then, the instability emerges, leading to periodic motion for u = 0.4, 0.5, 1.0 and 1.5 μm, or chaotic motion for u = 2.0 and 2.5 μm. The onset of the instability depends on the excitation amplitude. Fourier analysis of the surface motion is given in Fig. 24. The two synthetic quantities are the space-averaging of the amplitude of time-fft, and the time-averaging of the amplitude of spacefft. The amplitudes are normalized to facilitate a relative comparison.
The synthetic time-fft highlights a peak at 500 kHz, corresponding to the forcing frequency, and a peak at 250 kHz, corresponding to the ω/2 subharmonic Faraday instability. The synthetic space-fft shows a peak at a wavelength of 0.2e − 4 m, corresponding to the wavelength of the standing wave. Chaotic motion is demonstrated through widening of the spectrum for u = 2.0 and 2.5 μm. Peak amplitudes and instability onset delay are plotted in Fig. 25. As expected, the 500 kHz peak follows an almost linear relation with respect to the forcing amplitude. The 250 kHz peak drops at 2.0 μm and 2.5 μm because of peak widening. The onset delay is very dependent on the forcing oscillation, with values from tens of μs to hundreds of μs for the smallest excitation amplitude.

Discussion
The selected simulation method correctly predicts the gas bubble flattening in a consistent manner with the experiments. The flattened shape is only obtained when updating the acoustic force with respect to the moving boundary, but not as a solution of the static problem. The most flattened shape is obtained for the advected case, for which the gas bubble is squeezed between the liquid slug and the acoustic standing wave, see e.g. Video 11 frame No. 2600. The comparison between free, constrained and advected configurations allows to reveal the dependence of the gas bubble behavior on the surrounding medium conditions. Since it is expected that the volume of the liquid slug between bubbles does not change, or only changes relatively slightly, the constrained bubble configuration may be more realistic than the free bubble configuration. However, it has to be noted that the thin liquid film

Table 5
Simulation results for the oscillating gas bubble interface compared to experiments (Video 6).
⃒ ⃒ v 2 | in : maximum speed towards the inside of the gas bubble (acoustic force pushing the interface). between the gas bubble and the channel wall was not modeled, as it cannot be captured with the current modeling approach. Constrained and both-sides resonance configurations highlight a steady deformation state corresponding closely to what is observed in the microreactor. The gas bubble advected towards a resonance position exhibits a more complex behavior which may be encountered in practice. In this case, the resonance is met when the bubble reaches a particular position. Steady flattening is then obtained for a finite period, depending on the acoustic force magnitude. Furthermore, both experiments and simulations demonstrate that there is a coexistence between the flattened shape solution and the shock solution. Both solutions can generate atomization at the deformed interface. We assert that the existence of a shock solution, which leads to transient concave deformation, is due to a short acoustic peak of high amplitude. The acoustic peak depends on the resonance conditions. Relatively good agreement is obtained for the concave deformation of the bubble, as shown in Fig. 18. For the case of strong deformation (e.g. see Fig. 18 (c.1) and (c.2)), it may be expected that the simplified Lagrangian Surface Evolver reaches its limits, since more incompressible flow effects occur near the strongly curved surface. Moreover, the standing wave approximation may reach its limits when dealing with such short acoustic excitation. A time-domain acoustic solver may be more relevant to deal with strong deformation that takes place for hundreds of microseconds. Notwithstanding, the dynamic behavior is well reproduced, but might be overestimated due to the overestimation of the acoustic force by the frequency-domain solver.
The predicted acoustic pressure is in the range of hundreds of kPa to tens of MPa. A large acoustic pressure manifests itself first with gas bubble flattening at speeds up to about 2 m/s, then atomization can emerge. Cavitation is determined directly by the acoustic pressure. Assuming that, at 500 kHz, cavitation arises likely around 100 kPa, the quantitative approach also predicts weak to strong cavitation in the liquid. This is consistent with the experimental observation of Mc Carogher et al. [27], where cavitation bubbles are generated and focused at the pressure nodes of the acoustic mode and at the gas bubble surface. In summary, the comparison of the experiments and the numerical results shows that a stable acoustic pressure value can go up to 2 MPa, and that atomization arises from 10 MPa. We propose that a resonance peak can lead to values up to 80 MPa, which we expect are then responsible for strong atomization and cavitation happening for tens of microseconds, as observed for example in Video 7.
VOF results allow to characterize atomization conditions for this precise configuration. The measured critical wavelength of about 20 μm ( Fig. 24 (b)) is consistent with the K&T theory. This value can also be approximated by Kelvin's formula λ = 2π(4γ/(ρ 0 ω 2 )) (1/3) , where ω is the forcing pulsation. Slightly above the theoretical critical oscillation amplitude of 0.35 μm, for the 0.4 μm simulation, a weak Faraday instability is observed. One can assert that for all the standing acoustic waves whose particular oscillation amplitude is below the K&T theoretical threshold, no standing Faraday wave can be observed. The associated acoustic pressure for the selected 500 kHz configuration is about 1.5 MPa. According to VOF results, it seems that the breakup threshold is between 1.5 μm and 2.0 μm. The associated acoustic pressure range would be from 6 to 10 MPa. Assuming pure (0, n) mode, one can use the formula relating the acoustic pressure and the horizontal oscillation amplitude: p = u x c 0 ρ 0 ω. This gives an acoustic pressure of Goodridge's equation for the selected constants [47,48]: where a d is the critical driving acceleration of the onset of drop ejection. For 500 kHz, this formula gives a critical driving acceleration of 4.977e + 6 m/s 2 , thus an oscillation amplitude of 0.504 μm. Consequently, we found that this formula overestimates the atomization for the selected frequency. However, it has to be noted that this formula was validated in a frequency range from 25 Hz to 100 Hz in Puthenveettil and Hopfinger's work [43], and therefore might not be directly applicable in the ultrasonic range. To obtain an atomization threshold of about 1.75 μm, corresponding to the one predicted in this study, the coefficient in Eq. (21) should be about 0.9 instead of 0.26. Results for the simplified VOF configuration showed that the interface vibration pattern of mode (0, n) is quite homogeneous around the mid-channel point, as depicted in Fig. 9 first row second column. One can expect that more complex atomization patterns are associated with more complex vibration patterns, which may have nodes along the channel height.
Assuming that an acoustic pressure greater than 10 MPa is required to generate a sufficiently observable atomization, the numerical results demonstrate that atomization only happens during strong transient deformations associated with resonance. The simulated steady states  involve steady acoustic pressures below 3 MPa, which is only sufficient to produce at most a standing Faraday wave. This assertion is consistent with the experiments, see e.g. Video 5 frame No. 182 and Video 7 frame 288. Most of atomization appears as what one can call atomization bursts. Thus, we found that a possible steady atomization behavior would consist of a sequence of atomization bursts. This can happen in both the flattened and shock case. The shock solution, for the constrained gas bubble, can be periodic, since the gas bubble can recover its initial shape. For the oscillating solution, the dynamic behavior is well reproduced, as a similar oscillating shape is obtained. The predicted oscillation frequency is from 850 to 920 Hz, compared to an experimental     result of 1000 Hz (from Video 6). The recovery time is determined inherently by the surface tension to density ratio, but also by the magnitude of the acoustic force. Little adjustments of these parameters may give a better fit of the oscillation frequency. But given the simplicity of the assumptions leading to these results, we propose that the oscillation frequency discrepancy is acceptable. Furthermore, in practice, due to the fact that the acoustic resonance is highly sensitive to geometry changes, a perfectly periodic behavior is not observed. The typical recovery period obtained by the simulations can take values from 1 to 5 ms, see Fig. 19 (a), so irregular oscillations with variable period in this range can be expected.

Conclusions and perspectives
A numerical approach combining an acoustic frequency-domain solver and a Lagrangian Surface-Evolver solver was developed to describe the acoustic deformation of gas-liquid interfaces. The method was designed to describe phenomena driven by dynamic acoustic resonance. In such problems, nonlinearity and high dynamic range are expected. We focused on the acoustic resonance happening between gas bubbles in a segmented two-phase flow in a microchannel sonicated at 500 kHz. Globally, the behavior of the bubble submitted to acoustic resonance is well reproduced, with good qualitative agreement with experiment and consistent order of magnitude. The direct relation between acoustic field and interface deformation was established. We concluded that the coexistence of steady flattening and strong transient concave deformation of gas bubbles is consistent with selected hypotheses. We found that the acoustic pressure magnitude can take values from hundreds of kPa to tens of MPa. These values are consistent with the observation of atomization and cavitation in the experiments. VOF simulations allowed to highlight atomization generation conditions and were put into perspective with acoustic resonance. We found that dynamic acoustic resonance gives rise to atomization bursts at the gas bubble surface.
In the presented developments, we focused on a simplified configuration of the microreactor and we considered vibration homogeneous in space along the microreactor plane. The presented approach can be applied to more complex acoustic fields involving more complex channel geometry, the actual vibration pattern of the ultrasound transducer, or more complex two-phase flow. The method can also be extended to 3D to model Faraday crystals and the actual atomization pattern. For the selected configuration, the 3D extension would have led to similar conclusions regarding gas bubble deformation. The selected numerical framework allows not only to predict the phenomena numerically, but can also be used in inverse problems to characterize the design of ultrasonic microreactors.

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.   and where p a is the amplitude of the acoustic pressure, v a is the amplitude of the vertical acoustic oscillation velocity. The relation v a = k y p a /(ρ 0 ω) holds. One can demonstrate by applying both Navier-Stokes and Zakharov approaches, described in Section 2.4, with associated assumptions, that the steady (∂ t ∝0) normal displacement of the upper surface at y = H is: From Eq. (B.5), the characteristic length L in the case of the liquid film is the wavelength in the x direction. The acoustic deformation of fluids can be described synthetically by We acoustic . In the presence of dynamic acoustic resonance, We acoustic can vary greatly because of the large dynamic range of acoustic magnitude.