Abstract
The shape oscillation of a single two-dimensional nitrogen microbubble in an ultrasound field is numerically investigated. The Navier–Stokes equations are solved by using the finite-volume method combined with the volume-of-fluid model. The numerical results are in good accordance with experimental and theoretical results reported. According to the analyses of the shape oscillation process, the bubble deformation period is twice the driving acoustic pressure period and the shape oscillation is mainly caused by the change of interface velocity. The vortexes produced due to velocity variations lead to the deformation of the bubble interface.
1. Introduction
Two-phase flow occurs in a wide range of natural phenomena and engineering applications [1–4]. Among them, microbubbles have received much attention in the past decades, either for the energy that bubbles can release by cavitation [5–7] or for their utility as ultrasound contrast agents in medicine [8]. Recently, shape oscillations of a single microbubble in microfluidic chips in an ultrasound field have been reported by many researchers. Rabaud et al. [9] found that the bubble tended to form periodic crystal-like lattices within a finite interbubble distance when excited in an external acoustic field. From the figures and videos they provided, apparent shape oscillations of a single microbubble can be found, although the authors did not mention this. Liu et al. [10] found that the surface instability of a microbubble can lead to the bubble rupture, thus shortening the microbubble residence time. However, in gene therapy, the ultrasound-targeted microbubble destruction technique (UTMD) needs to induce the instability of the bubbles. A certain intensity of ultrasound is applied on specific microbubbles for carrying the target gene in the vicinity of the cells, causing microbubble cavitation collapse, then formatting holes in the cell membrane, finally promoting target gene into the cell, to achieve the purpose of gene transfection [11, 12]. Therefore, it is important to analyze the stability of a single microbubble in an ultrasound field.
So far, the characteristics of a single microbubble in the ultrasound field are mostly studied through experiments [13]. For the first time, Benjamin et al. [14] used an ultrahigh-speed imaging technique to systematically study aspheric oscillations and surface modes of bubbles. It was pointed out that the nonspherical deformation was caused by the parameter instability driven by radial oscillations. And the relationship between the driving frequency and the resonance frequency was concluded. In the study of the impact of bubble oscillation on the ultrasound cleaning, Kim and Kim [15] found that bubble behaviors in the ultrasound field could be divided into four types: volume oscillation, shape oscillation, collapse, and chaos oscillation. The shape oscillation was the manifestation of the bubble surface instability. In addition, they also gave a bubble behavior classification map, which is very helpful for studying bubble surface instability. Versluis et al. [16] experimentally studied the response characteristics of air bubbles in water. They found that there was a linear relationship between oscillation mode n and the bubble radius, while the mode n had no certain relationship with the amplitude of the driving pressure. McDougald and Leal [17] numerically simulated the nonlinear mode vibration of spherical bubbles under nonviscous conditions by using the boundary integral method. Recently, Liu et al. [10, 18, 19] numerically simulated the shape oscillation of the encapsulated microbubble under viscous conditions. Mekki-Berrada et al. [20] carried out a detailed investigation on the dynamics of two-dimensional (2D) bubbles at the wall interface of a microfluidic channel. By studying the response of the nitrogen bubble in the surfactant-added liquid to ultrasound, it was demonstrated that the wall deformation had no significant effect on the bubble dynamics characteristics. Furthermore, they found that, above a critical pressure threshold, the bubble exhibited a two-dimensional shape oscillation around its periphery with a period doubling characteristic of a parametric instability. These phenomena are interesting. However, it is difficult to obtain the detailed flow behaviors inside and around the oscillating bubble experimentally.
Therefore, the main objective of this work is to study the shape oscillations of a single microbubble in an ultrasound field based on numerical simulation, and the relationship between the bubble shape variation period and the driving pressure period is investigated.
2. Numerical Methods
In this study, a nitrogen microbubble in water in an ultrasound field is investigated numerically to analyze the characteristics of the bubble shape oscillation. The Navier–Stokes equations [21] are used, and the volume-of-fluid (VOF) model [22] is introduced to track the interface between gas and liquid. The following assumptions are adopted to simplify the model:(1)The liquid is considered to be incompressible, while the gas inside the bubble is compressible, conforming to the ideal gas law.(2)The gas phase is immiscible with the liquid phase, and the mass transfer between two phases is neglected.(3)The effect of gravity is ignored [23].
In the present study, a 2D simulation is carried out. The 2D microbubble is corresponding to the pancake bubble confined between two walls in the microchannel, as described in [20]. The schematic of the simulation domain is shown in Figure 1. In the simulation, the side length of the square domain is 500 μm. The 2D bubble is positioned in the middle of the domain, and the bubble radius is variable. The orthogonal quadrilateral grid is adopted, and mesh refinement is applied near the bubble. The grid independence is first carried out, and a total grid number of 150,000 is adopted in the simulation.
In all simulations, the operating pressure is set to be 1 atm and the initial velocity is zero. All boundaries are set as the pressure inlet. The inlet pressures follow the ultrasound pressure equation:where is the driving frequency and is the acoustic pressure amplitude. The initial pressure inside the bubble follows the equation aswhere is the initial gas pressure inside the bubble; is the initial liquid pressure; is the surface tension coefficient, and = 0.0735 N/m is used in this study; and is the initial bubble radius.
To solve the time-dependent Navier–Stokes equations, the finite-volume method [24] is applied. The volume fraction equation is solved through explicit time discretization. The PISO algorithm [25] is used for pressure-velocity coupling, and implicit time discretization is used. Numerical criteria are equally computed in double precision with a segregated solver. For all cases, the convergence criteria for continuity and momentum equations are set to be 10−5 and for energy equation, are set to be 10−7. The time step is chosen depending on the frequency of ultrasound.
3. Results and Discussions
3.1. Model Validation
The simulation results of different shape oscillation modes are first compared with experimental results of references [20] (as shown in Figure 2) and [26] to verify the model.
The dynamic behaviors of a microbubble at specified ultrasound frequency and acoustic pressure amplitude are first investigated. The ultrasound frequency is set to be 130 kHz, and the acoustic pressure amplitude is 42 kPa. The shape oscillation of microbubbles at different initial bubble radii is shown in Figure 3. In these snapshots of phase concentration contours, color blue stands for gas, while color red stands for liquid. These bubble shapes agree well with those experimental results of Mekki-Berrada et al. [20], as shown in Figure 2. The experimental figures presented in Figure 2 are carried out in a channel confined between two walls 25 μm apart. It is a pancake bubble. The boundaries of the bubble in the experiment are not flat. It is different from the 2D bubble in the simulation. So the numerical results seem to give smoother boundaries of bubbles than that of the reference.
(a)
(b)
(c)
(d)
(e)
(f)
The natural frequency of two-dimensional bubble shape oscillation was studied by Mekki-Berrada et al. [20]. They presented a theoretical formula:
In general, the bubble oscillation mode increased monotonously with the radius within a certain frequency. As shown in Figure 4, the relation between bubble modes and radius at 130 kHz is in good agreement with (3) [20]. Lamb [27] investigated 3D bubble oscillations and found the following equation:
The shape oscillation modes of a 2D bubble in simulations are slightly different with that of a 3D bubble. However, the overall trend remains the same. It can be found from the figure that the same oscillation mode occurs over a certain bubble radius range. For example, when f = 130 kHz and Pd = 42 kPa, the oscillation mode is 4 when the bubble initial radius is in the range of 26 to 30 μm.
The simulation results and theoretical results at four frequencies are displayed in Figure 5 to further verify the reliability of the simulation results. In our present simulations, when the bubble’s initial radius ranges from 20 μm to 60 μm, the bubble’s shape oscillation mode increases with an initial radius. It can be found from Figures 4 and 5 that the simulation results are in good agreement with the theoretical formula [20, 26].
4. Numerical Simulation of Shape Oscillation
The dynamic behaviors of bubble oscillation in the ultrasound field at different modes are similar. Without loss of generality, the shape oscillation processes of a single bubble at the fourth mode are discussed in detail.
A dimensionless time is first defined, = t/T. Here, T is the period of the driving ultrasound wave, T = 1/f. The bubble shapes in an oscillation cycle are shown in Figure 6 at = iT, (i + 0.5)T, (i + 1)T, (i + 1.5)T, and (i + 2)T. Here, we choose a circular bubble shape as the start of a shape oscillation period, as shown in the figure. During a whole shape oscillation circle, the microbubble changes from its initial circular shape into the four-sided cross shape in the first half driving ultrasound period (0 < < 0.5T). Then, the bubble returns to the circular shape in the next half driving ultrasound period (0.5T < < T). Later, it changes to a four-sided fork shape in the third half driving ultrasound period (T < < 1.5T) and returns to the circular shape when 1.5T < < 2T. The bubble’s shape variation period is exactly 2 times of the driving ultrasound pressure period.
(a)
(b)
(c)
(d)
(e)
In order to further explore the shape variation during the oscillation process, the internal and external gage pressure variations and the velocity vector near the gas-liquid interface are analyzed. The external liquid pressure is probed near the boundary, while the internal gas pressure is probed at the bubble center. The variations of liquid and gas pressures from = iT to = (i + 2)T are shown in Figure 7. The amplitude of the liquid pressure is approximate 42 kPa, which is the same as the driving ultrasound pressure. The amplitude of the gas pressure is approximately 18 kPa, which is much smaller than that of the liquid pressure. We can also find from Figure 7 that the fluctuating frequency of the gas pressure is exactly the same as that of the liquid pressure. However, there is a hysteresis in the gas pressure. The hysteresis time is about 0.45T.
According to the extreme outline of the bubble shape deformation, we presented the contrast chart of the instantaneous bubble morphology under the extreme state. As shown in Figure 8, with the bubble center as the origin, ρ(θ, ) is the distance from the bubble center to interface at time . Distinguishing bubble shape deformation is mainly based on the angle corresponding to the maximum value of ρ. The fork-shaped bubble appears when ρmax corresponds to , and when ρmax corresponds to , the cross-shaped bubble appears.
The instantaneous velocity vectors and bubble shapes in a single shape oscillation period are shown in Figure 9. The shape oscillation period is exactly twice of the driving ultrasound period. Let’s first take a look at the time = (i + 1.75)T. The bubble shape is approximately a square, as shown in the figure. The gage liquid pressure is at its maximum point, and the inner gas pressure is near the minimum, as can be found in Figure 7. There exists a pressure difference between the liquid and the gas across the interface, and therefore, the interface moves towards the bubble from four corners along the direction of θ1. Eight vortexes appear at the gas-liquid surface. As a result, the interface in the θ1 direction moves inward, and the interface in the θ2 direction moves outward.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Later, the external liquid pressure starts to decrease, and the internal gas pressure starts to increase. But the liquid pressure is still larger than the gas pressure, so the inward flow velocity in the θ1 direction keeps increasing. This velocity reaches its maximum at = iT, as shown in Figure 9. At this moment, the internal pressure equals roughly to the external pressure. The bubble shape is approximately a circle, which is the first stage in Figure 6. After this moment, the liquid pressure keeps decreasing and it is smaller than the gas pressure. The velocity in the θ1 direction begins to decrease. The interface in the θ1 direction keeps moving inward, and the interface in the θ2 direction keeps moving outward. The cross-shaped bubble begins to appear. The liquid pressure reaches its minimum, and the inner gas pressure is near the maximum at = (i + 0.25)T. The bubble begins to expand because of its high pressure. The velocity magnitude in the θ2 direction is larger than that in the θ1 direction.
The bubble reaches its extreme cross-shaped state at = (i + 0.5)T. At this moment, the inward velocity in the θ1 direction vanishes and these eight vortexes near the gas-liquid interface disappear. The gas and liquid pressures are nearly the same. The bubble reaches its maximum volume. After this moment, the outside liquid pressure is larger than the inner gas pressure. The bubble enters the compression stage. The interface in the θ2 direction begins to flow inward. As a result, it returns to the square shape at = (i + 0.75)T but with an angle of π/4 different from the bubble shape at = (i + 1.75)T.
Up to now, an ultrasound vibration period has finished. The bubble undergoes a compression and expansion process. In the next ultrasound vibration, the bubble undergoes another compression and expansion process and the bubble varies from circle to cross-shaped. The shape oscillation period is exactly twice the driving ultrasound period.
It can be concluded from Figure 9 that the shape oscillation of the bubble in an ultrasound field is mainly induced by the interface velocity and vortex, which is an intuitive result of the pressure difference between the outside liquid and inner gas across the gas-liquid interface. The instability of the bubble surface and the oscillation mode are combined effects of liquid viscosity, bubble diameter, and driving ultrasound parameters.
5. Conclusions
In the present paper, the response of a single 2D microbubble in an ultrasound field was investigated numerically. The Navier–Stokes equations are used, and the VOF mode was applied to capture the bubble interface between gas nitrogen and liquid water. Bubble’s shape instability was investigated, and the effects of the initial radius and ultrasound frequency on the unstable mode of the microbubble were discussed. It can be found from the simulations that numerical results are in good accordance with experimental and theoretical results. The shape oscillation mode increases with the increase of the bubble’s initial radius and driving ultrasound frequency. The same oscillation mode occurs over a range of the bubble radius. In order to explore the mechanism of bubble shape oscillation in the ultrasound field, the detailed velocity variation of the bubble in a single shape oscillation period was presented. The bubble oscillation period is exactly twice the driving ultrasound pressure period, and the shape oscillation of the bubble is mainly caused by the change of interface velocity. Variations in the velocity result in vortexes, which will lead to the deformation of the bubble interface.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (Grant nos. 11672284 and 11474259) and the National Key R&D Program of China (Grant nos. 2017YFB0603701 and 2016YFF0203302).