Abstract
This paper presents the development of a numerical model able to track in time the behavior of nonlinear focused ultrasound when interacting with tiny gas bubbles in a liquid. Our goal here is to analyze the frequency components of the waves by developing a model that can easily be adapted to the geometrical restrictions and complexities that come out in several application frameworks (sonochemistry, medicine, and engineering). We thus model the behavior of nonlinear focused ultrasound propagating in a liquid with gas bubbles by means of the finite-element method in an axisymmetric three-dimensional domain and the generalized- method in the time domain. The model solves a differential system derived for the nonlinear interaction of acoustic waves and gas bubble oscillations. The high nonlinearity and dispersion of the bubbly medium hugely affect the behavior of the finite-amplitude waves. These characteristics are used here to generate frequency components of the signals that do not exist at the source through nonlinear mixing (parametric antenna). The ability of the model to work with complex geometries, which is the main advantage of the method, is illustrated through the simulation of nonlinear focused ultrasound in a medium excited from two spherical sources in opposite directions.
1. Introduction
Finite-amplitude ultrasonic waves are of interest in many applications [1]. Specific mathematical models are required to analyze these nonlinear waves, especially when the medium in which they propagate is nonhomogeneous [2]. Bubbly liquids, for which small gas bubbles are introduced into a homogeneous liquid, have been intensively studied in the last few years. However, the nonlinear interaction between ultrasonic waves and bubble vibrations requires more theoretical studies to increase the knowledge about the underlying effects of such waves [3]. A liquid with a small bubble concentration is dispersive. Its sound speed, nonlinear parameter, and attenuation coefficient can be raised over determined frequency bands. This characteristic greatly affects the propagation of ultrasound [4–7]. In this paper, we study ultrasound of finite amplitude that focuses on a bubbly liquid by means of a numerical model. Finite-amplitude focused ultrasound can be useful in medicine (e.g., noninvasive diagnostic and therapeutic techniques) [8, 9] and industry (e.g., sonochemistry) [1, 10–12]. We carry out the study through a finite-element model (FEM). We consider a spherical transducer exciting a three-dimensional volume of liquid with bubbles evenly distributed. The wave equation, in terms of acoustic pressure variable, is used to model the acoustic field. We couple the wave equation to a Rayleigh–Plesset equation, written in terms of bubble volume variation, to model the bubble vibrations [4, 13, 14]. The numerical model developed to solve the differential system considers cylindrical symmetry around the transducer axis. It allows us to track the acoustic pressure in time. The time-dependent signal is then decomposed by the fast Fourier transform (FFT) to obtain its frequency components. The pressure amplitude imposed at the source is modified in the numerical experiments to compare the linear and nonlinear responses of the system.
Few models based on the FEM have been developed for the modeling of nonlinear ultrasonic fields in sonoreactors. Moreover, most of them are based on Helmholtz-type equations. Chu et al., in a recent paper [12], give a summary on this topic in their introduction, from the seminal papers of Dähnke et al. [15, 16], the works by Yasui et al. [17, 18], the work by Tudela et al. [19], the review papers by Tudela et al. [3] and by Lebon et al. [20], up to Sarac et al. [21], before comparing several Helmholtz-based models solved through the FEM.
It must be noted that some seminal studies on nonlinear focused ultrasound in bubbly liquids have been published [22–24]. They are based on the Khokhlov–Zabolotskaya–Kuznetsov (KZK) equation model (including some variants). They show the effect of one single bubble or a bubble population (evenly distributed or not) on the focused ultrasonic field through nonlinearity, diffraction, and absorption and also consider the thermal conduction [24]. None of them assumes a multifrequency excitation from the source, unlike our work presented in Section 3.2. Moreover, their modeling is based on equations valid for focused ultrasound in homogeneous liquid adapted to bubbly liquids through the introduction of equivalent parameters such as sound speed and attenuation. None of them takes the bubbles into account as a series of oscillators that induce nonlinearity, dispersion, and attenuation affecting ultrasound as we do here. Others study the behavior of the focus by taking the dispersion of a medium (not due to gas bubbles) into account through the resolution of KZK-based equations [25, 26]. It is worth noting that [25] shows the focus shift when dispersion is considered, but with a low coefficient of nonlinearity compared to a bubbly liquid, which is several orders of magnitude higher than the corresponding homogeneous liquid [5]. Other works study one-dimensional bubbly cavitating flows in elastic fluids, such as [27], which models this phenomenon through microsized nozzles of different shapes.
The FEM developed here in the time domain allows us to consider the study of the nonlinear effects generated from a dual-frequency source and to easily modify the geometrical configuration of our problem (Section 3).
This paper deals with the development of a numerical model able to track in time the behavior of nonlinear focused ultrasound when interacting with tiny gas bubbles in a liquid. Our goal here is to analyze the frequency components of the waves by developing a model that can easily adapt the geometrical restrictions and complexities in several application frameworks (sonochemistry, medicine, and engineering).
This paper is organized as follows. Section 2 describes the differential system derived for the physical problem of focused ultrasound in bubbly liquids we consider here. Section 3 presents illustrative simulations in a medium with an even density of bubbles by considering a specific physical configuration and parametric excitation for which the pressure amplitude at the source is varied. This section also shows simulations with a more complex geometry carried out easily from the previous ones. The conclusions of this work are described in Section 4.
2. Materials and Methods
We consider a liquid in an open cylindrical domain in which tiny gas bubbles are added. A spherical transducer emits an ultrasonic signal into the bubbly liquid (Figure 1). Bubbles oscillate under the effect of ultrasound. is the initial radius of the bubbles. is their initial volume. is their current volume. The spherical transducer is defined by its diameter and radius of curvature . The bubble density, i.e., number of bubbles per unit volume, , is constant in the entire domain. All the bubbles have the same size and are spherical. The three-dimensional system is assumed to have a symmetry around the axis defined by the main direction of acoustic propagation, . The other space coordinate in the axisymmetric plane is the radial coordinate . is the time-independent variable. and are the lengths of the space domain in the and directions. is the last instant of the study.
(a)
(b)
The coupled differential system that models the nonlinear interaction of the ultrasonic field, wave equation in terms of acoustic pressure , and the vibrations of the bubbles, Rayleigh–Plesset equation in terms of volume variation , is [4, 13]in which is the viscous damping coefficient of the bubble due to the viscosity of the fluid, with the kinematic viscosity of the liquid, is the bubble resonance, with the specific heats ratio of gas, the atmospheric pressure of gas, the equilibrium density of gas, the small-amplitude sound speed of gas, the equilibrium density of liquid, is the second-order nonlinear coefficient due to the adiabatic law used to model the gas pressure inside the bubble, is the nonlinear coefficient associated to the bubble dynamics, is a constant in the source term, and is the small amplitude sound speed of liquid. The system (liquid and bubbles) at rest at the onset of the study, the axial symmetry, and the open-field character of the space domain lead to the initial and boundary conditions that close the differential system:including the excitation at the source defined byin which is the pressure amplitude. Note that when two spherical sources are considered in Section 3.3, the above source equation also applies at . Some approximations are considered in this work: the bubble is spherical (spherically symmetric pulsation); the bubble content is adiabatic gas only, without vapor; the surface tension at the gas-liquid boundary is neglected; the bubble does not radiate sound itself; bubble oscillations are of moderate values; buoyancy, Bjerknes, and viscous drag forces are neglected.
In order to seize the geometrical and graphical aspects given by a commercial FEM solver, we chose to solve equations (1)–(5) with the COMSOL Multiphysics® software [28], through the definition of two “Coefficient Form PDE” coupled in the axisymmetric domain. To optimize the accuracy of the results in relation to the computation cost, this space domain is discretized with quadratic triangular Lagrange elements, and we consider two distinct areas. The first one, called the focal region in which the energy of both the pressure and bubble volume variation fields is concentrated and large gradient values of dependent variables are expected, is an ellipse centered on the focal point according to the radius of curvature. In this focal region, a dense mesh is employed to ensure that the sharp gradients in both fields are approximated correctly. A maximum mesh size of is chosen (where is the shortest wavelength of the pressure field imposed at the curvature radius). The second one that encompasses the rest of the space domain uses 6 elements per wavelength. This number is a reasonable compromise between computation time and accuracy. Figure 2 displays the meshes used to carry out the simulations in the following sections. The number of degrees of freedom (DOFs) ranges from around 36,000 to around 470,000.
(a)
(b)
The generalized- scheme is used to discretize the problem in the time domain [29]. Compared to backward differentiation formula, this implicit method is known to lead to minimal numerical dissipation and to enhance the stability and convergence of the approximate solution. The time discretization step used in the following simulations is , where is the sound speed in the bubbly liquid at the given frequency and is the maximum element size of the mesh, defined previously in the largest part of the domain as the sixth part of the wavelength This leads to the stability of the simulations carried out hereafter. The Courant–Friedrichs–Levy number is set at 0.1 in all the simulations shown in this paper. It is worth noting that simulations carried out with give roughly the same results in the time domain, but the quality of the FFT of the two dependent variables is somehow altered.
The coupling of the wave equation to the bubble equation induces dispersion, which mitigates the formation of shock waves. This important feature of the coupled differential system allows us to observe harmonic components of the pressure wave without strong difficulties. Moreover, the high nonlinear characteristic of the medium allows it to trigger nonlinear behaviors at quite low-pressure amplitudes, which makes it very much easy to model as well. The densification of the FEM mesh, where the highest acoustic-pressure gradients are physically expected definitely ensures the good behavior of the simulation results (see previous paragraphs and Section 3.1).
3. Results and Discussion
The numerical model described above is used here to solve the physical problem given in Section 2 considering a medium made of water and air bubbles. The following parameters are set: m, m, , m, m, s, ( MHz, MHz), /m3, m/s, kg/m3, m/s, kg/m3, , and m2/s.
The dispersive curves of the bubbly medium defined above with bubbles evenly distributed are shown in Figure 3 [4]. They allow us to localize the frequency ranges for which nonlinearity and attenuation are high or low, in order to set the working frequencies according to what we want to prevail in the simulations carried out in Section 3. The vertical red lines indicate the frequency values which are useful in these simulations (from left to right in the diagram): , , , , , and ; dashed, dash-dotted, and solid lines correspond to, respectively, primary, second harmonics of primary, and mixing-created frequencies. belongs to a huge attenuation frequency range. This feature is also valid, in a lower measure, for and . and are located in a frequency range of quite high compressibility, which means that the nonlinearity of the medium is high for these frequencies [4].
This section is structured as follows. We first consider a single-frequency pressure source (Figure 1(a), Section 3.1), for which a very low-pressure amplitude is employed. We then use a high-amplitude pressure to observe the nonlinear effects due to the presence of bubbles. In the second part of this section (Figure 1(a), Section 3.2), a parametric study is done by assuming a dual-frequency pressure source in two configurations: linear and nonlinear. Finally, the potential use of the model, due to its attractive feature concerning the modification of the geometry, is commented in Section 3.3 (Figure 1(b)). Note that the useful model developed in this paper, and tested in this section, will allow us to simulate more complex nonlinear ultrasonic fields easily in the future (more complex geometries, nonevenly distributed bubbles in the liquid).
3.1. Single-Frequency Source
The source is driven at the single frequency kHz and amplitude Pa by means of the continuous signal , where is the frequency of the source (Figure 1(a)). The travelling wave in the liquid impinges the bubbles. The mesh, made of 9,454 elements (36,736 DOF, 2,304 internal DOF), is presented in Figure 2(a).
Figure 4 shows the acoustic pressure vs. time (a) and vs. frequency (b), the bubble volume variation vs. time (c) and vs. frequency (d), all of them at the focus, i.e., the point of maximal pressure amplitude on the axis, and the acoustic pressure amplitude after 10 periods in the axisymmetric plane (e). Note that Diagram (e) represents the whole (z, r) domain, not only its half-calculated part, for a better observation of the focalization effect. The same presentation of the results applies to Figures 5–8. Although s are calculated, the FFT shown in these figures are evaluated from the last 65 periods. Only few periods are displayed in the time representations to clearly observe the behavior of the waveforms.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
The focal zone is situated between m and m from the source. The pressure and bubble volume variation waveforms are symmetric around Pa and m3, respectively. Only one significative frequency component, the source frequency, is seen in both spectra.
Now, we work with the source amplitude kPa in the same configuration as above. Figure 5 shows the corresponding results. Due to the high-pressure amplitude, the nonlinearity of the medium at the working frequency (Figure 3) implies modifications of the ultrasonic field.
Although the space pattern is quite similar to the linear one and the focal zone is still situated between m and m from the source, differences occur at the focus, where acoustic energy is concentrated. Both waveforms are no longer symmetric around their horizontal axis: they are shifted to the positive pressure values and to the negative bubble volume variation values, respectively. The broadening of the negative pressure zone is also observed. Besides the source frequency , the high nonlinearity of the medium at (Figure 3) generates the second , third , and fourth harmonics, seen in the spectra. Their respective pressure amplitude values are about , , and of the fundamental frequency amplitude . This means that the pressure wave is strongly nonlinear.
A mesh refinement procedure has been carried out for the simulations shown in this paper. We present the results only in this particular nonlinear case, but similar results about convergence have been obtained in the configurations presented in the following sections. Figure 6 shows the maximal relative and values obtained vs. the number of DOFs used in the mesh (DOF variation is from 2,000 to 65,000). The reference values and used here are the ones obtained from a 100,000 DOF mesh. Also, it is important to mention that the amplitude values employed in this paper (Sections 3.1, 3.2, and 3.3) are high enough to observe the physical phenomena that we wanted to describe in this work and that the number of DOFs used in these sections is within the stable region of the convergence curve (Figure 6).
3.2. Dual-Frequency Source
We now excite the medium with the following dual-frequency continuous source, , where and are the primary frequencies of the source (Figure 1(a)). In this section, kHz and kHz. We set the pressure amplitude at the source at Pa (Figure 7). Since the working frequency is higher than in Section 3.1, the mesh, presented in Figure 2(a), is now made of 118,418 elements (469,684 DOF, 8,024 internal DOF).
The acoustic and bubble volume variation responses are completely different from the single-frequency case presented in the previous section. Although the focal zone is still situated between m and m from the source, the volume of high acoustic energy is much smaller than in the single-frequency case. Both waveforms are linear. They include the two primary frequencies shown in the spectra.
Now, we aim to analyze the nonlinear mixing of the primary frequencies while they are travelling through the medium. We set the pressure amplitude at the source at kPa. We draw special attention to the difference and sum frequency components (Figure 8). Note that, from equation (5) and the function defined in a previous paragraph, the pressure amplitude at the source can reach 8 kPa.
We adopt the following notation. is the difference-frequency component. is the pressure amplitude of . is the sum-frequency component. is the pressure amplitude of . Since kHz and kHz, kHz, kHz, kHz, and kHz.
The location of the focus in this high-amplitude case is still the same as in the linear case (between m and m from the source). However, the pressure waveform presents a quite strong asymmetry around Pa, showing a shift towards the positive pressure values. Moreover, besides and , the spectra include the components ( of the pressure amplitude at the primary frequency ), ( of the pressure amplitude at the primary frequency ), and ( of the pressure amplitude at the primary frequency ), as well as a slight component. The high nonlinearity of the medium at and (Figure 3) is such that the components and are visible and quite intense (a slight component is even noted).
This model is thus able to track in time the behavior of ultrasound in the bubbly medium in nonlinear configurations. It has a strong advantage over other models [14, 30, 31], since the FEM allows us to modify the geometrical conditions of the physical problem easily. This important benefit will be the topic of future works. In the following section, we use the model described above to illustrate that asset, considering a system with a double spherical source excited at a single frequency.
It must be noted here that to verify and validate the numerical model presented in this paper, it was used with the parameters given in Sections 3.1 and 3.2 and compared with another model based on a different kind of numerical methods, published earlier in the literature [30, 31]. Although the numerical model in the above references did not allow to get the same geometrical shape at the source, this comparison led to a quite good concordance. With the single-frequency source, the maximal pressure amplitude at with the model from [30] was of the fundamental frequency amplitude , whereas it is here. In the case of the dual-frequency source, the maximal pressure amplitude at with the model from [31] was of the pressure amplitude at the primary frequency , whereas it is with the current model presented here.
3.3. Double Spherical Single-Frequency Source
This section is an illustration of the easy applicability of the model presented and tested in the previous sections in relation to the modification of the geometry of the space domain. We consider a double spherical source emitting in the bubbly liquid at a single frequency in opposite directions (Figure 1(b)). With this system, we can check whether the double-source configuration induces nonlinear effects at the focus spot easier than the single-source configuration (Section 3.1) or not, i.e., whether a lower intensity level (lower electrical power) at the source is required to generate a second harmonic or not. The geometry, shown with the mesh used here in Figure 2(b), is easily constructed from one of the previous sections (Figure 2(a)). The mesh is made of 9,428 elements (36,642 DOF, 2,282 internal DOF).
We use the same parameter values as in Section 3.1. The diagrams represent the solution in the same way as in the above sections. Figure 9 shows the linear results, when Pa is set at the double source. Figure 10 shows the results in the nonlinear regime, when kPa is used at the double source. In both cases, the focal zone is situated at m. With Pa, the linear waveforms are sine-like functions and the spectra show only one frequency component at the source frequency . However, when kPa is imposed at the sources, the nonlinear waveforms at the focus present a clear asymmetry and the spectra reveal a quite strong second harmonic component. The quite low-pressure amplitude at the source leads to a nonlinear response of the system, with a significant second pressure harmonic component of amplitude kPa, i.e., of the sum of the pressure amplitude at the sources. The corresponding results obtained with the one-source model described in Section 3.1 using the same amplitude at the source kPa reveals an amplitude of the second harmonic Pa, i.e., of the pressure amplitude at the source. When k Pa is used at the source, kPa, i.e., of the pressure amplitude at the source. These data show that the system with two sources has a lower nonlinear efficiency. The higher efficiency of the single-spherical source system is most likely due to the propagation of a higher amplitude wave from the source which leads to a strong nonlinear process within the high nonlinear bubbly medium. This result could be interesting for the design of nonlinear systems based on focused waves in bubbly liquids.
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
4. Conclusions
The model based on the FEM developed here allows us to simulate the behavior of nonlinear focused ultrasound propagating in a liquid with gas bubbles in an axisymmetric three-dimensional domain. It solves a differential system derived for the nonlinear interaction of acoustic waves and gas bubble oscillations. The simulations carried out here have shown that high-amplitude ultrasound is affected by the high nonlinearity and dispersion of the bubbly medium. The model has been used to generate frequency components of the signals that do not exist at the source through nonlinear mixing. The ability of the model to work with complex geometries has also been shown here through the simulation of nonlinear focused ultrasound from a double-source system.
It is worth noting that the versatility of the FEM proposed here is expected to make further developments easier, specifically those concerning modifications of the geometrical characteristics of the physical problem, as well as changes on bubble-density distribution.
Data Availability
Data presented in the paper are all obtained through numerical simulations carried out with the model described in this article.
Conflicts of Interest
The authors declare that they have no conflicts of interest regarding the publication of this article.
Acknowledgments
This work was supported by the National Agency for Research, Ministry of Science and Innovation of Spain, and the European Regional Development Fund (Grant number DPI2017-84758-P), by the Ministry of Economy and Competitiveness of Spain (Grant nos. DPI2012-34613, BES-2013-064399, and EEBB-I-15-09737), and by the Artois University (Faculty of Applied Sciences) through its Visiting Professors program in 2021 and 2022.