Dynamics study of a laser-induced bubble on a finite metallic surface in water

Pulsed laser ablation in liquid (PLAL) has been widely studied because of its diversified applications, such as underwater laser machining,¹ treatment for kidney stones,² laser propulsion in water,³ and producing nanoparticle colloids with ligand-free surfaces.⁴ Cavitation bubbles, as a common attendant phenomenon of PLAL, have garnered widespread interest among researchers. Understanding bubble dynamics is important for avoiding or decreasing its destructive consequences,⁵ such as unwanted cavitation noise, influence on hydrodynamics, and destructive erosion. When a high-intensity laser pulse is focused on a material in water, some material at the laser’s focal area melts and evaporates after absorbing the laser energy, and plasma is formed at the laser’s focal area due to optical breakdown.^6,7 The plasma expansion is always accompanied by the emission of a shock wave and the generation of a bubble.^8,9 At first, the bubble continues to expand outward until all its initial kinetic energy is converted into potential energy. Then, the bubble collapses under the external pressure of the liquid. As a result, both the temperature and the pressure inside the bubble increase again, leading to its rebound, along with a shock wave.¹⁰ In general, the bubble oscillates several times before its final collapse,¹¹ and a liquid jet is generated after the bubble’s final collapse. Up until now, the theoretical and experimental investigations on bubble dynamics are mainly focused on bubbles generated in an infinite liquid,^11–13 near various boundaries^14,15 or on a rigid boundary.^16,17 Research in the dynamics of bubbles induced on a finite rigid boundary, especially when bubble size closely matches the size of the boundary, is currently insufficient.

In order to investigate the dynamics of a bubble induced on a finite rigid boundary, the time dependence of the bubble’s radius needs to be experimentally measured. The commonly used methods for detecting bubble dynamics include high-speed photography,^9,18,19 shadow photography,^20,21 schlieren photography,^22,23 holography-based photography,^24,25 measurements based on laser beam deflection probe (BDP),^26–28 measurements based on laser beam transmission probe (LBTP),^29–32 and some other methods.^33,34 All the mentioned methods, excluding LBTP technology, require substantial repeatability of the measuring process for obtaining the time dependence of the bubble’s radius, which requires a large workload and statistical measuring of the results. Although high-speed photography can obtain results with only one measurement, it requires an expensive high-speed camera to perform. LBTP technology is able to detect the whole bubble’s dynamics from a single shot, but it cannot detect its 2D spatial distribution. Another drawback of LBTP is that the probe beam does not have a top hat profile, meaning non-uniform probe-intensity distribution.³² In the case that the bubble’s position and geometry are known, it is possible to obtain the time dependence of the bubble’s radius by employing LBTP and shadow photography simultaneously,^29,32 or by combining LBTP and the corresponding algorithm.³⁰ However, the set-up of the former method is complex and the accuracy of the latter method significantly depends on the quality of the probe-intensity’s Gaussian TEM₀₀ distribution. Thus, an appropriate experimental method needs to be developed.

In this paper, we present a simple measuring method, based on the principle of LBTP technology, to detect the dynamics of a bubble induced on a finite metallic surface in water. The time dependence of the bubble’s radius can be obtained from a single shot, and a high-speed photographic method was employed as a comparative method to verify the validity of this method. A correction factor based on Rayleigh collapse time formula was proposed for describing the collapse time of the bubble induced on a finite rigid boundary, and its expression was obtained from the simulation and experimental results. In addition, the energy transfer during the first three oscillations of the bubble induced on the metallic surface was roughly analyzed.

When an isolated spherical bubble collapses in an infinitely large, incompressible, and inviscid liquid with the assumption that the liquid’s pressure is constant during its collapse, the variation of the bubble’s radius with time can be roughly described by the Rayleigh model:¹² 

Here, R is the bubble radius, P_R is the pressure in the bubble, P_∞ is the surrounding liquid pressure, and ρ is the liquid density.

Based on Eq. (1), the bubble energy E_B and the collapse time of the bubble T_C are deduced by Rayleigh:

Here, R_max is the maximum bubble radius in the Rayleigh model, and P_v is the saturated vapor pressure. Eq. (3) describes the relationship between the collapse time and the maximum radius. The Rayleigh model is a relatively simple model for bubble dynamics, but it can still describe the first oscillation of the bubble induced on an infinite rigid boundary. And this can be experimentally verified by the data in Refs. 16 and 17.

However, if the bubble is induced on a finite rigid boundary, especially when the bubble’s size approaches the size of the rigid boundary, Eq. (3) is not suitable and should be modified. To simplify the problem, we assume that the bubble is induced on the center of the finite rigid boundary (Fig. 1). Then, a dimensionless parameter γ, γ = R_max / L, is introduced to represent the spatial relationship between the bubble and the finite rigid boundary. The collapse time T_CF of the bubble induced on a finite rigid boundary can be described as:

Here, κ is the correction factor which is a function of γ, and it will be obtained from the experimental and simulation results in Sec. V. In this work, ρ is 1 × 10³ kg/m³, P_∞ is 1.01 × 10⁵ Pa, P_v is 2.34 × 10³ Pa, and L is 2 mm.

Since the shape of the bubble induced on a rigid boundary is quasi-semispherical, its bubble energy E_BF can be described as 0.5 E_B approximately, that is:

In order to obtain the relationship between the two parameters, κ and γ in Eq. (4), the collapse of a bubble induced on an infinite and finite rigid boundary was simulated with the VOF method. For simplifying the models of numerical simulation, we begin establishing the models when the bubble expands to the maximum radius first time, which means the simulation results may only be valuable for the bubble’s first collapse. The VOF method is based on the volume fraction equation Eq. (6) and the conservation laws: mass conservation Eq. (7), momentum conservation Eqs. (8) and (9), and energy conservation Eq. (10):

Here, α is the volume fraction, $u→$ is the velocity field, ρ is the liquid density, t is the time, x and r is axes, u and v are velocities along axes x and r, respectively, p is the pressure, and T is the temperature.

The models of numerical simulation are shown in Fig. 2: (a) the collapse of a bubble induced on an infinite rigid boundary, (b) the collapse of a bubble induced on a finite rigid boundary, which are called model a and model b in following sections, respectively. The unit of the spatial coordinates is mm, and x is the rotational symmetric axis of the whole model. Region D is the computational domain. Point (0,0) is the bubble center (the laser focus). A, B and C are the pressure boundaries, which are set to 1.01 × 10⁵ Pa. The solid object in Fig. 2(a) is 25 mm high (in r axis), which is approximated as an infinite rigid boundary because it is far larger than maximum bubble radius. The solid object in Fig. 2(b) is 2 mm high (in r axis) and serves as a finite rigid boundary. Both the two solid objects are 3 mm wide (in x axis).

There are two phases existing in the computation domain D: One is water and the other is vapor in the bubble. When the bubble reaches maximum size, the pressure of its content is equal to the vapor pressure P_v of the water,¹⁸ so the initial pressure of the bubble is set to saturated vapor pressure at room temperature 2.34 × 10³ Pa. The environmental pressure is set to 1.01 × 10⁵ Pa.

We developed an experimental system based on LBTP technology for measuring the time dependence of the bubble’s radius on a finite metal surface in water, shown in Fig. 3. It is mainly composed of two parts: the bubble generation system and the optical detection system. In the bubble generation system, a Q-switch Nd: YAG laser (Beamtech Dawa-200, wavelength 1064 nm, pulse width 7 ns) was used to induce the laser plasma and subsequent bubble. A cylindrical metal object (Φ 4 mm × 5 mm, Titanium) with a plane surface was fixed in a glass cuvette filled with distilled water. In order to prevent multiple optical breakdowns, the pulsed laser beam was first expanded to achieve a large focusing angle in water. In addition, the metal object was placed in front of the laser’s focal area to prevent water breakdown.

The optical detection system (Fig. 3) is the most important part of this system. It was fixed on a two-dimensional movable platform with 10 μm spatial resolution, which allows X and Y direction movement of the detection beam. All optical components were aligned on one horizontal line. In this system, the detection beam was emitted from a He-Ne laser (Melles Griot, 25-LHP-925-230, 632.8 nm, 17 mW). Its intensity profile has a Gaussian (TEM₀₀) distribution, and the beam diameter is 0.96 mm. In the detection area, this system offered a ‘—’ shaped detection beam with an almost even intensity distribution in the X direction, and the plasma center (i.e., the center of the bubble) was on the ‘—’ shaped detection beam during the measuring process (shown in the Inset of Fig. 3). At the receiving end, the detection beam was focused into a photomultiplier which was connected to an oscilloscope (Tektronix THS730A). So the optical signal was transformed into electrical signal, fed into the oscilloscope, and stored in a computer. The interference filter filtered light with wavelength outside of 632.8 nm. The circle aperture was used to block the detection beam deflected by the bubble and improve measurement accuracy.

After a bubble was generated on the surface of the metal object, it locally changed the density of the water and consequently its refractive index. It resulted in the deflection of the part of the detection beam going into the bubble, which changed the amount of the detection beam at the receiving end and led to the change of the signal shown by the oscilloscope. Thus, the bubble’s radius can be reflected in real time by the oscilloscope signal.

The details on the experimental method will be presented in the  Appendix, including the verification of the intensity distribution uniformity of the “—” shaped detection beam, comparison with high-speed photographic method, and error analysis.

Fig. 4 shows the relationship between the maximum bubble radius R_max and the collapse time T_CF for the bubble’s first three oscillations on a finite metallic surface under seven different values of incident laser energy E_L. Interestingly, it seems R_max and T_CF present a same linear relationship for the bubble’s first three oscillations, which is shown by the fit curve (dash line) in Fig. 4. In addition, the fit curve is obviously different with the solid line derived from Eq. (3), which means that the correction factor κ in Eq. (4) is not equal to 1. Thus, In order to obtain the expression of κ, we analyzed the fit curve. Through the linear fit, T_CF can be described as:

Since γ = R_max / L (L is 2 mm in the experiment), κ can be obtained by combining Eq. (4) and Eq. (11):

However, Eq. (12) is questionable when γ is very small or larger than 1. In this work, the surface tension was neglected since the maximum bubble radius is ∼mm and the Weber number is considerably larger than 1.¹⁶ But in the case that the bubble size is very small, some factors, such as surface tension and viscosity, cannot be neglected, and Eqs. (3), (4) and (12) are no longer suitable. When γ is larger than 1, the maximum bubble size is larger than the size of the finite metallic boundary. In this case, the bubble divides mainly into two parts during its collapse, resulting in the decrease of the collapse time.³⁵ 

As for the bubble’s first oscillation, we found that the first three experimental data (from left to right) agreed well with Eq. (3) (shown in Fig. 5), so it is possible that the experimental results of the bubble’s first oscillation with smaller bubble size will still agree with Eq. (3), which can be validated by the simulation results shown in Fig. 5. In other words, for smaller values of γ during the bubble’s first oscillation, κ is equal to 1, which can be experimentally validated by the data in Refs. 16 and 17. Thus, Eq. (12) needs to be improved, and κ(γ) should be different for bubble’s first oscillation and subsequent two oscillations.

The simulation results in Fig. 6 show the whole evolution of the correction factor κ with the dimensionless parameter γ (γ < 1) for bubble’s first oscillation. The experimental results agree with the simulation results. From Fig. 6, κ is approximately equal to 1 when γ is smaller than 0.6, and begins to decrease when γ is larger than 0.6. In other words, for the bubble’s first collapse, the relationship between the collapse time and maximum bubble’s radius can be described by the Rayleigh collapse time formula when γ is smaller than 0.6. In the γ > 0.6 case, the boundary effect of the finite metallic surface starts to have an impact on the bubble’s oscillation, and Eqs. (4) and (12) are more appropriate in describing this scenario.

In summary, for the first collapse time of the bubble induced on a finite metallic surface in water, the correction factor κ₁ based on the Rayleigh collapse time formula can be roughly described as:

Fig. 7 shows the the correction factor κ as a function of the dimensionless parameter γ for the bubble’s second and third oscillations. It shows that the correction factor κ for the bubble’s second and third oscillations both agree with Eq. (12). However, it is obviously impossible that the collapse time is infinite when γ is close to 0, so Eq. (12) is not suitable for the case that γ is close to 0. In this work, we focused on the dimensionless parameter γ in the range of 0.1∼1. Regarding the second and third collapse time of the bubble induced on a finite metallic surface in water, the correction factor κ_2,3 based on the Rayleigh collapse time formula can be roughly described as:

Since the time dependence of the bubble’s radius can be obtained by our measuring method, it is possible to make a preliminary energy analysis for the evolution of the bubble induced on a metallic surface based on Eq. (5). The energy analysis would help explain the effect of bubble dynamics on the metallic surface. In this work, we calculated the bubble energies of the first three bubble oscillations, which are represented by E_B1, E_B2 and E_B3 respectively. The ratios of the energy left for the first three bubble oscillations, which are η₁ = E_B1 / E_L, η₂ = E_B2 / E_B1 and η₃ = E_B3 / E_B2, under seven different values of E_L are shown in Fig. 8. We observed that both η₁ and η₂ increase with E_L (shown clearly in the inset of Fig. 8), but the variation trend of η₃ is not clear (Fig. 8). Additionally, the ratio of the energy left for subsequent bubble oscillation increases with the number of bubble oscillation, i.e., η₁ < η₂ < η₃. η₁ and η₂ are also not high values. In summary, after a high-intensity laser pulse is focused on a metallic surface in water, most of the incident laser energy E_L is converted into the heat energy for melting and evaporating the metal, and only a small portion of E_L (< 6% for E_L < 18.87 mJ in our experiment) is converted into the bubble energy. Then, the bubble begins to oscillate. When the bubble rebounds the first time, a shock wave is radiated, and only a small portion of the bubble energy (< 14% for E_L < 18.87 mJ in our experiment) produces the second bubble oscillation. When the bubble rebounds the second time, a shock wave is radiated again, and a larger portion of the bubble energy of the second bubble oscillation (< 65 % for 4.21 mJ < E_L < 18.87 mJ in our experiment) produces the third bubble oscillation.

In order to investigate the dynamics of a bubble induced on a finite rigid boundary, we developed a simple experimental method based on the principle of laser beam transmission probe (LBTP) technology to measure the time dependence of the bubble’s radius on a finite metallic surface under different incident laser energies E_L. High-speed photographic method was employed as a comparative method to verify the validity of this method. A correction factor κ based on the Raleigh collapse time formula was proposed to describe the collapse time of the bubble induced on a finite rigid boundary. κ is a function of the dimensionless parameter γ, which is the ratio of the maximum bubble’s radius to half of the finite rigid boundary when the bubble is induced on the center of the rigid boundary. The experimental and simulation results show that κ is slightly different for the bubble’s first collapse and subsequent two collapses. Its detailed expression was obtained from the experimental and simulation results. Furthermore, the experimental results show that during the process of laser-metal interaction in water, the conversion efficiency of E_L into the bubble energy increases with E_L, and the ratio of the energy left for subsequent bubble oscillation increases with the number of bubble oscillations. The results found in this paper are relevant for laser machining and laser propulsion in water environment.

This work is supported by the National Natural Science Foundation of China for Young Scholars (No.11402120, No.11502116), the Jiangsu Natural Science Foundation for Young Scholars (No. BK20140796), and the Fundamental Research Funds for the Central Universities (No.30915015104, No.30916014112-015).

As shown in Fig. 3, in order to obtain the ‘—’ shaped detection beam, the detection beam was first expanded in the X direction by a combination of two cylindrical concave lenses. It subsequently went through a rectangular aperture (3.9 mm × 0.2 mm), and a rectangular detection beam was formed. In order to reduce the possibility of light diffraction, we use a 0.2 mm-wide rectangular aperture. Then another cylindrical lens was used to focus the rectangular detection beam in the Y direction. It does not change the length of the detection beam. We made sure the focal line of the detection beam was perpendicular to the surface of the metal and the detection beam is parallel to the surface of the metal. Consequently, a ‘—’ shaped detection beam was formed in the detection area.

The relevant schematic diagram is shown in Fig. 9. It can be seen from Fig. 9 that the detection beam’s intensity distribution in the X direction is almost uniform theoretically after these measures.

Before we took measurement, it was necessary to verify the intensity distribution uniformity of the detection beam in the X direction, and do calibration, i.e., build the relationship between the length of the blocked detection beam and the strength of the electrical signal shown in the oscilloscope. First, we moved the detection beam to the left side of the metal object, and the metal object does not block the detection beam. Second, we moved the detection beam towards the metal object until the detection beam contacts the metal object. This was accomplished by observing the amplitude of the signal displayed on the oscilloscope during the movement. The amplitude of the signal will begin to decrease once the detection beam contacts the metal object. Then we recorded the position X₀ of the detection beam in the X axis and the amplitude V₀ of the electrical signal displayed on the oscilloscope. In this case, the length S of the detection beam blocked by the metal object was 0. Third, we continued to move the detection beam towards the metal object by 50 μm, and recorded the amplitude V₁ of the signal displayed on the oscilloscope. In this case, the length S of the detection beam blocked by the metal object was 0.05 mm. We repeated this step until the detection beam was entirely blocked by the metal target. In this condition, the amplitude of the signal displayed on the oscilloscope will not decrease any more. Fig. 10 shows the relationship between the length S of the blocked detection beam and the strength of the electrical signal displayed on the oscilloscope. It can be seen that the intensity distribution of the detection beam has a good uniformity in the X direction. In addition, the slope K of the fit line is -1.5014, and the maximum absolute error is about 27μm.

After calibration, we began the measurement of the laser-induced bubble on the metallic surface. First, the detection beam was moved to the left of the metal object once again. Then, we moved the detection beam towards the metal object until a small part of the detection beam was blocked by the metal object, and recorded the amplitude V_max of the electrical signal displayed on the oscilloscope. When a bubble is induced on the surface of the metal object after a high-intensity laser pulse, the electrical signal V(t) displayed on the oscilloscope can be obtained, and it can be transformed into the bubble’s radius R(t) by:

Fig. 11 shows a typical R(t) curve obtained by our method. It can be seen that the maximum bubble radius R_max and oscillation time T_OS for the bubble’s first three oscillations can be obtained by only one measurement. Furthermore, in our experiment, since the duration of the laser pulse is much shorter than the bubble’s oscillation time and the viscosity of the liquid is negligible, the expansion and the collapse of the bubble are almost symmetrical processes. Therefore, the bubble collapse time T_C can be expressed as T_C = 1/2 T_OS.²⁸ The laser-induced plasma shock wave can also be detected by our method (please see the dashed box in Fig. 11). To verify the validity of our method, the high-speed photographic method (SVSI GigaView)³⁶ was employed as a comparative method to measure the bubble’s radius over time when the incident laser energy E_L is 13.98mJ. The comparison between the measurement results obtained by our method and high-speed photographic method is shown in Fig. 11. There is a clear agreement between the measuring results obtained by the two methods, supporting the validity of our method.

The accuracy of our method is mainly determined by the intensity distribution uniformity of the detection beam in the X direction and the intensity stability of the detection beam. For the former, the error can be estimated from the calibration process. The maximum absolute error was about 27 μm. For the latter, the intensity stability of the detection beam was good enough and the relevant error was smaller than 20 μm since the intensity of the detection beam was very stable (<0.5%) during a 3 hours test. It is possible that a fraction of the detection beam, which goes through the center of the bubble, was transmitted through the bubble. As a result, the measuring result might be somewhat smaller that the real result, so we lengthened the distance between the metal object and the photomultiplier and applied a circle aperture between them. In this case, the error resulted from this reason can be neglected. In our experiment, the maximum radius of the bubble was in the millimeter range, so this measuring method can be used.

Since the time dependence of the bubble’s radius can be obtained by only one measurement, our system does not require the repeatability of the measuring process, which reduces the workload significantly. Furthermore, the time resolution of this system is determined by the sampling frequency of the oscilloscope, and it has the potential to be applied to studying faster phenomena. Compared to LBTP technology, our method has the advantage of having a detection beam with uniform intensity distribution, and boasts a simpler setup for measuring the time dependence of a bubble’s radius.

Disadvantages of our system include its weak accuracy for measuring smaller bubbles (∼100 μm), high requirement for optical alignment, and inability to detect a bubble’s 2D spatial distribution.