Shape Oscillation of a Single Microbubble in an Ultrasound Field

,e shape oscillation of a single two-dimensional nitrogen microbubble in an ultrasound field is numerically investigated. ,e Navier–Stokes equations are solved by using the finite-volume method combined with the volume-of-fluid model. ,e 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. ,e vortexes produced due to velocity variations lead to the deformation of the bubble interface.


Introduction
Two-phase flow occurs in a wide range of natural phenomena and engineering applications [1][2][3][4]. Among them, microbubbles have received much attention in the past decades, either for the energy that bubbles can release by cavitation [5][6][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 ultrasoundtargeted 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]. erefore, 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. e 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. ey 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 twodimensional (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. ese phenomena are interesting. However, it is difficult to obtain the detailed flow behaviors inside and around the oscillating bubble experimentally. erefore, 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.

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. e Navier-Stokes equations [21] are used, and the volume-offluid (VOF) model [22] is introduced to track the interface between gas and liquid. e following assumptions are adopted to simplify the model: (1) e liquid is considered to be incompressible, while the gas inside the bubble is compressible, conforming to the ideal gas law. (2) e gas phase is immiscible with the liquid phase, and the mass transfer between two phases is neglected. (3) e effect of gravity is ignored [23].
In the present study, a 2D simulation is carried out. e 2D microbubble is corresponding to the pancake bubble confined between two walls in the microchannel, as described in [20]. e schematic of the simulation domain is shown in Figure 1. In the simulation, the side length of the square domain is 500 μm. e 2D bubble is positioned in the middle of the domain, and the bubble radius is variable. e orthogonal quadrilateral grid is adopted, and mesh refinement is applied near the bubble. e 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. e inlet pressures follow the ultrasound pressure equation: where f is the driving frequency and P d is the acoustic pressure amplitude. e initial pressure inside the bubble follows the equation as where p g0 is the initial gas pressure inside the bubble; p l0 is the initial liquid pressure; c is the surface tension coefficient, and c � 0.0735 N/m is used in this study; and R 0 is the initial bubble radius. To solve the time-dependent Navier-Stokes equations, the finite-volume method [24] is applied. e volume fraction equation is solved through explicit time discretization. e PISO algorithm [25] is used for pressurevelocity 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 . e time step is chosen depending on the frequency of ultrasound.

Model Validation.
e 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. e dynamic behaviors of a microbubble at specified ultrasound frequency and acoustic pressure amplitude are first investigated. e ultrasound frequency is set to be 130 kHz, and the acoustic pressure amplitude is 42 kPa. e 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. ese bubble shapes agree well with those experimental results of Mekki-Berrada et al. [20], as shown in Figure 2. e experimental figures presented in Figure 2 are carried out in a channel confined between two walls 25 μm apart. It is a pancake bubble. e 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. e natural frequency of two-dimensional bubble shape oscillation was studied by Mekki-Berrada et al. [20]. ey 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: e 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 P d � 42 kPa, the oscillation mode is 4 when the bubble initial radius is in the range of 26 to 30 μm.
e 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].

Numerical Simulation of Shape Oscillation
e 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/T. Here, T is the period of the driving ultrasound wave, T � 1/f. e bubble shapes in an oscillation cycle are shown in Figure 6 at t * � 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 < t * < 0.5T). en, the bubble returns to the circular shape in the next half driving ultrasound period (0.5T < t * < T). Later, it changes to a foursided fork shape in the third half driving ultrasound period  (T < t * < 1.5T) and returns to the circular shape when 1.5T < t * < 2T. e bubble's shape variation period is exactly 2 times of the driving ultrasound pressure period. 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. e external liquid pressure is probed near the boundary, while the internal gas pressure is probed at the bubble center. e variations of liquid and gas pressures from t * � iT to t * � (i + 2)T are shown in Figure 7. e amplitude of the liquid pressure is approximate 42 kPa, which is the same as the driving ultrasound pressure. e 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. e 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, ρ(θ, t * ) is the distance from the bubble center to interface at time t * . Distinguishing bubble shape deformation is mainly based on the angle corresponding to the maximum value of ρ. e fork-shaped bubble appears when ρ max corresponds to θ 1 � ((2i − 1)/4)π, and when ρ max corresponds to θ 2 � (1/2)π, the cross-shaped bubble appears. e instantaneous velocity vectors and bubble shapes in a single shape oscillation period are shown in Figure 9. e shape oscillation period is exactly twice of the driving ultrasound period. Let's first take a look at the time t * � (i + 1.75)T. e bubble shape is approximately a square, as shown in the figure. e gage liquid pressure is at its maximum point, and the inner gas pressure is near the minimum, as can be found in Figure 7.
ere 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.
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. is velocity reaches its maximum at t * � iT, as shown in Figure 9. At this moment, the internal pressure equals roughly to the external pressure. e 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. e velocity in the θ 1 direction begins to decrease. e interface in the θ 1 direction keeps moving inward, and the interface in the θ 2 direction keeps moving outward. e cross-shaped bubble begins to appear. e liquid pressure reaches its minimum, and the inner gas pressure is near the maximum at t * � (i + 0.25)T. e bubble begins to expand because of its high pressure. e velocity magnitude in the θ 2 direction is larger than that in the θ 1 direction. e bubble reaches its extreme cross-shaped state at t * � (i + 0.5)T. At this moment, the inward velocity in the θ 1 direction vanishes and these eight vortexes near the gasliquid interface disappear. e gas and liquid pressures are nearly the same. e bubble reaches its maximum volume. After this moment, the outside liquid pressure is larger than the inner gas pressure. e bubble enters the compression stage. e interface in the θ 2 direction begins to flow inward. As a result, it returns to the square shape at t * � (i + 0.75)T but with an angle of π/4 different from the bubble shape at t * � (i + 1.75)T.
Up to now, an ultrasound vibration period has finished. e 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. e 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. e instability of the bubble surface and the oscillation mode are combined effects of liquid viscosity, bubble diameter, and driving ultrasound parameters.

Conclusions
In the present paper, the response of a single 2D microbubble in an ultrasound field was investigated numerically. e 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. e shape oscillation mode increases with the increase of the bubble's initial radius and driving ultrasound frequency. e 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. e 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
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.