Long-Pulse Laser-Induced Cavitation: A Race Between Advection and Phase Transition

Figure 2 illustrates the problem investigated in this paper, showing a test case that generates a pear-shaped vapor bubble. The computational analysis is designed to start at the time when laser is just activated. At this time, the fluid domain is completely occupied by liquid water. The analysis is expected to predict the localized water vaporization due to laser radiation, the subsequent bubble and fluid dynamics, and the dissipation of laser energy in this two-phase fluid medium. Therefore, we solve the following compressible inviscid Navier-Stokes equations that include radiative heat transfer.

with

Here, $\Omega\subset\mathbb{R}^{3}$ denotes the domain of the fluid flow. Open subsets $\Omega_{0}$ and $\Omega_{1}$ represent the subdomains occupied by the liquid and vapor phases, respectively. They are time-dependent. In most experiments, a sharp boundary between the liquid and vapor phases can be clearly captured by high-speed cameras. Therefore, we assume $\Omega_{0}\cap\Omega_{1}=\emptyset$. $\rho$, $\bm{V}$, $p$, and $T$ denote the fluid’s density, velocity, pressure, and temperature, respectively. $e_{t}$ is the total energy per unit mass, given by

where $e$ denotes the fluid’s internal energy per unit mass. $k$ is the thermal conductivity coefficient, which takes different values in $\Omega_{0}$ and $\Omega_{1}$. $\bm{q_{r}}$ denotes the radiative heat flux induced by laser. Equation (1) needs to be closed by a complete equation of state (EOS) for each phase, including a temperature equation. The computational model and solver utilized in this study supports arbitrary convex EOS (Ma et al., 2023). In this study, we adopt the Noble-Abel stiffened gas equation (Le Métayer & Saurel, 2016) for both phases. Specifically,

in which the subscript $\mathcal{I}\in\{0,1\}$ identifies the liquid ($0$) and vapor ($1$) phases. For each phase, $\gamma$, $p_{c}$, $q$, and $b$ are constant parameters that characterize its thermodynamic properties. For example, a non-zero $b$ allows the model to have a variable Grüneisen parameter that depends on $\rho$. Clearly, (4) is a generalization of perfect gas, stiffened gas, and Noble-Abel equations of state (Le Métayer & Saurel, 2016).

For a given pressure equation like (4), the choice of temperature equation is not unique. We adopt the one proposed in Le Métayer & Saurel (2016), i.e.,

where $c_{v}$ denotes the specific heat capacity at constant volume, assumed to be a constant. It can be shown that the specific heat capacity at constant pressure, $c_{p}$, is also a constant, given by $c_{p}=\gamma c_{v}$. Combining (4) and (5) gives a complete EOS that satisfies the first law of thermodynamics.

Two groups of EOS parameter values are tested in this study, as shown in Table 1. Neither of them was calibrated specifically for laser-induced cavitation (Le Métayer & Saurel, 2016; Zein et al., 2013, see). We will show in Section 4 that the simulation result is indeed influenced by these parameter values.

We model the bubble surface as a sharp interface with zero thickness. It is defined by

On the interface, we assume continuity of normal velocity and pressure, i.e.

where $\bm{n}$ denotes the normal to $\Gamma$.

$\Gamma$ is time-dependent, and must be solved for during the analysis. We represent it implicitly as the zero level set of a signed distance function, $\phi$, defined in the closure of $\Omega$. That is,

where $\overline{\Omega}$ denotes the closure of $\Omega$.

In this way, the aforementioned phase identifier, $\mathcal{I}$, is given by

The evolution of $\Gamma$ in time is driven by both phase transition and advection. The advection of $\Gamma$ by the fluid flow is governed by the level-set equation,

where $\bm{V}$ is the flow velocity.

At the beginning of the analysis, $\Omega=\Omega_{0}$, and $\Gamma=\emptyset$. Therefore, we initialize $\phi$ to be a constant positive value everywhere in the domain, and start solving (10) only after phase transition starts. The detection and handling of laser-induced phase transition will be discussed in Section 2.4.

Let $\Omega_{L}\subset\Omega$ denote the region in the fluid domain that is exposed to laser. The laser generators used in this study have a flat surface with a small diverging angle. Therefore, $\Omega_{L}$ is in the shape of a truncated cone (Fig. 2). Within $\Omega_{L}$, energy conservation implies

where $\mathcal{L}=\mathcal{L}(\bm{x},\hat{\bm{s}},\eta)$ denotes the spectral radiance (dimension: [mass][time]^-3[length]^-1) at position $\bm{x}\in\mathbb{R}^{3}$, along direction $\hat{\bm{s}}$, at wavelength $\eta$. $\mu_{\alpha}$ and $\mu_{s}$ denote the coefficients of absorption and scattering, respectively. They depend on the laser’s wavelength and the medium. $\mathcal{L}_{b}$ denotes the blackbody radiance. $\Phi(\bm{\hat{s}_{i}},\bm{\hat{s}})$ is the scattering phase function, which gives the probability that a ray from one direction $\bm{\hat{s}_{i}}$ would be scattered into another direction $\bm{\hat{s}}$. If $\hat{\bm{s}}$ is independent of $\bm{x}$, the left-hand side of Eq. (11) simplifies to $\nabla\mathcal{L}\cdot\bm{s}$, which gives the well-known radiative transfer equation (RTE) (Modest, 2013; Howell et al., 2020).

The radiative heat flux $\bm{q}_{r}$ in Eq. (1) is obtained by integrating $\mathcal{L}$ over all directions and the interval of relevant wavelengths, $(\eta_{\text{min}},~{}\eta_{\text{max}})$. That is,

The special properties of laser light allows us to simplify (11) and (12). The intensity of the laser light is much higher than the blackbody radiance. Therefore, we assume

Also, $\mathcal{L}$ is non-zero only along the direction of laser propagation ($\bm{s}$) and at the fixed laser wavelength. That is,

where $\delta(\cdot)$ denotes the delta function, and the variable $L(\bm{x})$ on the right-hand side is radiance (dimension: [mass][time]^-3). With these assumptions, a laser radiation equation is obtained by substituting Eq. (14) and Eq. (13) into Eq. (11), i.e.,

Here, $\bm{s}=\bm{s}(\bm{x})$ is a known function, but not a constant unless the laser light has a parallel beam. Similarly, substituting Eq. (14) into Eq. (12) yields

In this study, we measure the time history of laser power, $P_{0}(t)$, in laboratory experiments. We assume the laser radiance on the surface of the laser fiber (i.e., $\Gamma_{LS}$ in Fig. 2) follows a Gaussian distribution, also known as a Gaussian beam (Welch et al., 2011). Therefore, on $\Gamma_{LS}$ we have

where $\bm{x}_{0}$ denotes the center of $\Gamma_{LS}$, and $w_{0}$ is the waist radius. (17) serves as a Dirichlet boundary condition for the laser radiation equation, (15).

We employ the method of latent heat reservoir proposed in Zhao et al. (2023) to predict laser-induced vaporization. As shown in Fig. 3, the fundamental idea of this method is to introduce a new variable, $\Lambda(\bm{x},t)$, to track the intermolecular potential energy in the liquid phase. The vaporization temperature, $T_{\text{vap}}$, and latent heat, $l$, are assumed to be constant and used in the method as coefficients. When the analysis starts, $\Lambda(\bm{x},t)$ is initialized to $0$ everywhere. As the liquid absorbs laser energy, temperature increases gradually. At any point $\bm{x}\in\Omega_{0}$, once $T_{\text{vap}}$ is reached, additional heat — due to radiation, advection, or diffusion — is added to $\Lambda$. When $\Lambda$ reaches $l$, phase transition occurs. In Fig. 3, this time instant is denoted by $t_{4}$. We assume that phase transition completes instantaneously at each point through an isochoric process. The state variables after phase transition are given by

where the superscripts $-$ and $+$ indicate the variable’s value before and after phase transition, e.g.,

In (22), $\Delta x$ denotes the local element size in the computational mesh. The pressure and temperature after phase transition are obtained naturally from the EOS of the vapor phase, i.e.,

The pressure rise, $p^{+}-p^{-}$, drives the expansion of the vapor bubble.

At each time step, if phase transition occurs, the level set function $\phi$ is restored to a signed distance function by solving the reinitialization equation,

to the steady state. Here, $\tilde{t}$ is a fictitious time variable. $\phi_{0}$ is the level set function before reinitialization. $S(\phi_{0})$ is a smoothed sign function, given by

where $\varepsilon$ is a constant coefficient, set to the minimum element size of the mesh. The steady state solution of Eq. (26) is then used as the new initial condition to integrate the level set equation (10) forward in time.

We apply a recently developed laser-fluid computational framework to solve the above model equations (Zhao et al., 2023). The fluid governing equations are semi-discretized using a fixed, non-interface conforming finite volume mesh, denoted by $\Omega^{h}$ (Fig. 4). The laser fiber is modeled as a fixed slip wall, and the associated boundary condition is enforced using an embedded boundary method (Cao et al., 2018, 2021b; Wang, 2017; Cao et al., 2021a). Therefore, $\Omega^{h}$ also covers the region occupied by the laser fiber. Around each node $\bm{n}_{i}\in\Omega^{h}$, a control volume $C_{i}$ is created. Applying the standard finite volume spatial discretization to Eq. (1) yields

where $\bm{W}_{i}$ denotes the average value of $\bm{W}$ in $C_{i}$. $Nei(i)$ denotes the set of nodes connected to node $\bm{n}_{i}$. $\partial C_{ij}=\partial C_{i}\cap\partial C_{j}$. $\bm{n}_{ij}$ is the unit normal to $\partial C_{ij}$. $|C_{i}|$ denotes the volume of $C_{i}$.

The FIVER method is used to calculate the advective flux across each facet $\partial C_{ij}$. Depending on the locations of nodes $\bm{n}_{i}$ and $\bm{n}_{j}$, there are four different cases.

1. [(a)]

2. 1.

   If $\bm{n}_{i}$ and $\bm{n}_{j}$ belong to the same phase subdomain ($\Omega_{0}$ or $\Omega_{1}$), the conventional monotonic upstream-centered scheme for conservation laws (MUSCL) is used to compute the flux across $\partial C_{ij}$.

3. 2.

   If $\bm{n}_{i}$ and $\bm{n}_{j}$ belong to different phase subdomains, a one-dimensional bimaterial Riemann problem is constructed along the edge between $\bm{n}_{i}$ and $\bm{n}_{j}$, i.e.,

   $$\frac{\partial\bm{w}}{\partial\tau}+\frac{\partial\mathcal{F}(\bm{w})}{\partial\xi}=0,\quad\text{with}\quad\bm{w}(\xi,0)=\begin{cases}&\bm{w}_{i},\qquad\text{if}\ \xi\leq 0,\\ &\bm{w}_{j},\qquad\text{if}\ \xi>0,\end{cases}$$

   (29)

   where $\tau$ denotes the time coordinate. $\xi$ denotes the spatial coordinate along the axis aligned with $\bm{n}_{ij}$ and centered at the phase interface between $\bm{n}_{i}$ and $\bm{n}_{j}$. The initial state $\bm{w}_{i}$ and $\bm{w}_{j}$ are the projections of $\bm{W}_{i}$ and $\bm{W}_{j}$ on the $\xi$ axis. This 1D two-phase compressible flow problem can be solved exactly (Ma et al., 2023). The solution is self-similar, and satisfies the interface conditions (7). The states to the left and right of the interface are used in the numerical flux function to compute the advective fluxes across $\partial C_{ij}$.

4. 3.

   If one of $\bm{n}_{i}$ and $\bm{n}_{j}$ belongs to a fluid phase subdomain ($\Omega_{0}$ or $\Omega_{1}$), while the other one is covered by the laser fiber, a 1D half-Riemann problem is constructed and solved exactly (Cao et al., 2021a; Wang, 2017). The solution at the wall boundary is used to compute the advective fluxes across $\partial C_{ij}$.

5. 4.

   If both $\bm{n}_{i}$ and $\bm{n}_{j}$ are covered by the laser fiber, the flux across $\partial C_{ij}$ is set to zero.

The right-hand side of Eq. (28) includes two parts, namely the heat diffusion term ($\displaystyle\int_{C_{i}}\nabla\cdot(k\nabla T)d\bm{x}$) and the radiation term ($\displaystyle\int_{C_{i}}\nabla\cdot\bm{q_{r}}d\bm{x}$). The heat diffusion term is evaluated using a finite volume method, i.e.

where $A_{ij}$ is the area of facet $\partial C_{ij}$. $(k\nabla T)_{ij}\cdot\bm{n}_{ij}$ denotes the heat flow due to diffusion cross the facet $\partial C_{ij}$. In this work, it is computed by

with

Here, $k_{i}$ and $k_{j}$ denote the thermal conductivity at nodes $\bm{n}_{i}$ and $\bm{n}_{j}$, specifically. We assume that the liquid-vapor interface is isothermal, and it can be shown that if $\bm{n}_{i}$ and $\bm{n}_{j}$ belong to different phase subdomains, Eqs. (31) and (32) enforce energy balance at the interface, with interface temperature (Faghri & Zhang, 2006)

We set the thermal conductivity at all the nodes covered by the laser fiber to zero. In this way, Eqs. (31) and (32) enforce the adiabatic boundary condition at the surface of the laser fiber.

The radiation term is evaluated by

At each time step, the divergence of the radiative flux, $\nabla\cdot\bm{q_{r}}$, is computed using the current solution of the laser radiation equation (15). In this work, the laser radiation equation is discretized using the same finite volume mesh ($\Omega^{h}$) created for the fluid governing equations (Fig. 4). Specifically, integrating (15) within an arbitrary cell $C_{i}$ yields

where $L_{i}$ is the cell-average of $L$ within $C_{i}$. $\bm{s}_{ij}$ denotes the laser direction at $\partial C_{ij}$, which is set to the laser direction at the midpoint between nodes $\bm{n}_{i}$ and $\bm{n}_{j}$. $L_{ij}$ is the value of $L$ at facet $\partial C_{ij}$. In this work, it is evaluated by the mean flux method, i.e.

where $\alpha\in[0.5,1]$ is a numerical parameter. Substituting Eq. (36) into Eq. (35) yields a system of linear equations with the cell-averages of laser radiance as unknowns. The Gauss-Seidel method is applied to solve this system to get $L_{i}$. Then, $\nabla\cdot\bm{q_{r}}$ is obtained by

When solving the laser radiation equation, the fluid mesh $\Omega^{h}$ does not contain a subset of nodes, edges, and elements that resolve the boundary of $\Omega_{L}$ or the laser propagation directions $\bm{s}(\bm{x})$. In this work, the embedded boundary method proposed in Zhao et al. (2023) is adopted to resolve this issue. This method involves the population of ghost nodes outside the side boundary of the laser domain using second-order mirroring and interpolation techniques.

The numerical methods described above are implemented in the M2C solver (Wang et al., 2021), which is used to run the simulations reported in this paper. Several verification studies can be found in Islam et al. (2023); Zhao et al. (2023); Cao et al. (2021a), and earlier papers cited therein. An outline of the solution procedure within each time step is provided below. For simplicity, the Forward Euler method is assumed here for time integration. In the actual simulations presented in this paper, a second-order Runge-Kutta method is used.

As mentioned in Sec. 4.1.1, the laboratory experiment reveals a delay of bubble nucleation, that is, the time of bubble nucleation is clearly after the time when the vaporization temperature ($T_{\text{vap}}$) is reached at the fiber tip. The same phenomenon is observed in the simulation.

Figure 9 shows the temperature field at $8$ time instants during the early period of the simulation. At the beginning, the temperature of water in front of the laser fiber increases continuously, as it absorbs laser energy. This can be seen in the solution snapshots taken at $0$ to $7~{}\text{\textmu s}$. At $7~{}\text{\textmu s}$, the temperature of water next to the center of the fiber tip first reaches $T_{\text{vap}}$. This is about $10~{}\text{\textmu s}$ before bubble nucleation, which occurs at $17.4~{}\text{\textmu s}$. Within this time period, the region that reaches $T_{\text{vap}}$ expands continuously, but phase transition does not occur.

This time delay is due to the fact that water has a high latent heat of vaporization. Based on the values of $l$ and $c_{v}$ adopted in the simulation, the latent heat is about $8$ times the energy needed to raise water temperature from $T_{0}$ to $T_{\text{vap}}$. Using the phase transition model described in Sec. 2.4, as soon as $T_{\text{vap}}$ is reached, any additional heat contributes towards increasing the intermolecular potential energy, thereby overcoming the required latent heat of vaporization. To examine this process more closely, we define

which measures the total amount of latent heat in the computational domain. Figure 10 shows the time-history of $\Lambda_{\Omega}$. At $6.9~{}\text{\textmu s}$, $\Lambda_{\Omega}$ becomes non-zero and begins to increase. Up to the time of bubble nucleation ($17.4~{}\text{\textmu s}$), the total energy stored in the latent heat reservoir is around $6.3$ mJ. By integrating the power profile (Fig. 7(d)), we find that this value is approximately $16.84\%$ of the laser energy input up to the same time.

In this simulation, vaporization starts at $17.4~{}\text{\textmu s}$, and continues for a short period of time (less than $1~{}\text{\textmu s}$). Figure 10 shows that after vaporization stops, $\Lambda_{\Omega}$ drops to around $0.8$ mJ. Given that the total laser pulse energy is $0.2$ J, this result implies that only a small fraction of the laser energy input ($(6.3~{}\text{mJ}-0.8~{}\text{mJ})/0.2~{}\text{J}=2.75\%$) is directly used to create the bubble.

After vaporization, the temperature inside the vapor bubble is up to $2000$ K, as shown in the solution snapshots at $17.4~{}\text{\textmu s}$ and $17.6~{}\text{\textmu s}$ in Fig. 9. The temperature near the bubble interface is higher than that in other regions, which is related to the direct energy input to the bubble after vaporization stops as depicted in Fig. 10. Subsequently, the temperature inside the bubble gradually decreases as the bubble expands, owing to the combined effects of advection and diffusion.

In both the experiment and the simulation, a pear-shaped vapor bubble is obtained. To explain the formation of this shape, we look at the pressure and velocity fields obtained from the simulation (Figures 11 and 12).

First, the two pressure snapshots at $0.4~{}\text{\textmu s}$ and $1~{}\text{\textmu s}$ (Fig. 11) capture an outgoing acoustic wave that emanated from the fiber tip. This is a thermal shock caused by the rapid increase of laser power. The vapor bubble nucleates at $17.4~{}\text{\textmu s}$. The pressure and velocity fields at this time are shown in Figs. 11 and 12. At this time, the bubble is already non-spherical. The pressure inside the bubble is high and nonuniform. Specifically, the pressure at the forward end is around $500$ MPa. The pressure at the backward end — that is, near the fiber tip — is much lower, around $100$ MPa. This pressure variation is a result of the brief continuation of vaporization. Regions that have just undergone phase transition have high pressures. Then, the pressure drops as the bubble expands. Therefore, the continuation of vaporization along the beam direction leads to the higher pressure at the forward end of the bubble. In summary, the simulation result at the time of bubble nucleation ($17.4~{}\text{\textmu s}$) already implies that the bubble will likely grow into a non-spherical shape, longer in the laser beam direction.

The high pressure inside the initial bubble also generates a weak shock wave that propagates outward at approximately the speed of sound in water. This wave can be seen in the pressure snapshots taken at $17.6~{}\text{\textmu s}-19~{}\text{\textmu s}$ (Fig. 11). Afterwards, the pressure field becomes quiet, and the bubble continues growing due to inertia. The velocity field at $50~{}\text{\textmu s}-120~{}\text{\textmu s}$ (Fig. 12) shows that the velocity around the bubble is nonuniform. It is higher at the forward end of the bubble compared to other regions. Furthermore, the velocity distribution within the bubble exhibits significant non-uniformity, as evidenced by the snapshots captured at time instances of $50~{}\text{\textmu s}$ and $120~{}\text{\textmu s}$. The velocity is notably higher in the vicinity of the central axis of the fiber, particularly near the forward end. This also explains the formation of a pear-like shape, instead of a sphere.

We show that the result obtained from our laser-fluid coupled computational framework can be influenced by the choice of EOS, and more precisely, the choice of EOS parameter values in this case. Here, we present another simulation with a different group of EOS paramete values, listed as Group $2$ in Table 1. All the other (physical and numerical) parameters remain the same.

Figure 13 compares the solutions obtained from the two simulations at five time instants. In each sub-figure, the left half of the figure is the solution obtained with Group $1$ parameter values, and the right half is the solution obtained with Group $2$ parameter values. Overall, both simulations predict the nucleation of a vapor bubble over a very short period of time. After that, the high pressure inside the bubble drives its expansion. In both simulations, a rounded bubble is obtained at the end.

However, a few differences can be found between the two simulations. First, there is a difference in the time when vaporization occurs. With Group $1$ parameter values, vaporization takes place at $17.4~{}\text{\textmu s}$. With Group $2$ parameter values, it occurs later, at $18.8~{}\text{\textmu s}$. The laser parameters are the same in both cases. Therefore, this discrepancy should be attributed mainly to the different temperature increase rates determined by the EOS, as defined in Eq. (5). Moreover, the speed of bubble growth is found to be lower with Group $2$ parameter values than with Group $1$. In Fig. 13, at both $25~{}\text{\textmu s}$ and $30~{}\text{\textmu s}$, the bubble obtained with Group $1$ parameter values is clearly larger than that obtained with Group $2$. The difference in growth speed can be explained by the pressure field. As shown in Fig. 13 (c), the bubble’s internal pressure reaches a maximum of $73~{}\text{MPa}$ with Group $2$ parameter values. This is much lower than the maximum pressure obtained with Group $1$ parameter values, which is $500~{}\text{MPa}$. Furthermore, the difference in bubble size also influences the delivery of laser energy. The larger bubble obtained with Group $1$ parameter values allows the laser energy to be delivered further, as demonstrated in Fig. 13(a). Lastly, the shapes of the bubbles obtained from the two simulations are slightly different. A pear-shaped bubble is captured with Group $1$ parameter values. With Group $2$, the bubble is more rounded. This discrepancy arises because in the latter case, no significant velocity difference is found at the bubble surface, as shown in Fig. 13(d).

To investigate the mechanism of bubble elongation, we examine the temperature, pressure, and velocity fields obtained from the simulation.

Figure 17 presents the evolution of the temperature field near the fiber tip in the first $5~{}\text{\textmu s}$. The laser radiance field is also shown at three time instants ($1.2$, $2.0$, and $2.4~{}\text{\textmu s}$) to facilitate the discussion. Similar to the previous case (Fig. 9), water temperature increases due to the absorption of laser energy, especially in the region around the central axis of the laser beam. The time it takes to reach the vaporization temperature is significantly shorter, only $0.4~{}\text{\textmu s}$, compared to $7~{}\text{\textmu s}$ in the previous case. The faster temperature increase can be attributed to two factors. First, the smaller beam waist of the laser results in a more concentrated distribution of laser radiance on the source plane. Although the laser power is lower in the current case (compare Figs. 8(d) and 16(b)), the more concentrated distribution leads to a higher laser radiance along the central axis, that is, $110~{}\text{kW/mm}^{2}$ (Fig. 16(a)) versus $80~{}\text{kW/mm}^{2}$ in the previous case (Fig. 8(a)). Second, the absorption coefficient of TFL in water is significantly higher than that of the Ho:YAG laser used previously ($14~{}\text{mm}^{-1}$ versus $2.42~{}\text{mm}^{-1}$). The time delay in bubble nucleation is also observed in this case. After $T_{\text{vap}}$ is reached, it takes another $0.8~{}\text{\textmu s}$ before vaporization occurs at the fiber tip, at $1.2~{}\text{\textmu s}$. Within this time period, the absorbed laser energy is converted to the intermolecular potential energy of liquid water.

A major difference from the previous case is that with the TFL, vaporization continues for a much longer period of time, that is, from $1.2~{}\text{\textmu s}$ until $53.5~{}\text{\textmu s}$. Also, it happens mainly along the central axis of the laser beam, which drives the bubble to grow in the same direction.

Initially, a small, rounded bubble emerges in front of the laser fiber. This is shown in Fig. 17, in the snapshot taken at $1.2~{}\text{\textmu s}$. Because the vapor phase does not absorb laser energy ($\mu_{\alpha}$ set to $0.001~{}\text{mm}^{-1}$), the small bubble extends the laser beam along the axial direction. This effect can be seen in the inset images in Fig. 17. As a result, the liquid water next to the bubble’s forward end — that is, the forward-most point along the axial direction — experiences a sudden increase in laser radiance, which accelerates the accumulation of energy there to overcome the latent heat of vaporization. From the simulation result, we observe the continuation of phase transition in this direction. For example, the snapshot taken at $2.0~{}\text{\textmu s}$ captures a bulge at the bubble’s forward end, which is a newly vaporized region. At the same time, the high pressure inside the bubble also drives it to expand in both axial and radial directions. The combination of these two processes — that is, phase transition and advection — drives the bubble to grow into a long, conical shape.

Figure 18 presents the evolution of the pressure field up to $70~{}\text{\textmu s}$. At the beginning, the sudden increase of temperature due to the absorption of laser leads to a thermal shock at the fiber tip. The snapshot taken at $1~{}\text{\textmu s}$ captures this phenomenon, where the maximum pressure is found to be less than $0.5~{}\text{MPa}$. Next, the snapshot taken at $1.2~{}\text{\textmu s}$ captures the pressure field inside and around the initial bubble. The peak pressure inside the bubble is found to be $94~{}\text{MPa}$ at this time. This high pressure drives the bubble to expand in all directions. It also generates a weak shock wave that propagates outwards, at approximately the speed of sound in liquid water. Because of the small size of the initial bubble, the pressure magnitude of this wave quickly decays to less than $1$ MPa. Compared to the pressure field obtained from the Ho:YAG laser (Fig. 11), the main difference here is that as phase transition continues, acoustic waves keep emanating from the bubble’s forward tip. This phenomenon is captured by all the snapshots taken between $2~{}\text{\textmu s}$ and $50~{}\text{\textmu s}$. In the current simulation, phase transition stops at $53.5~{}\text{\textmu s}$. Afterward, the pressure field becomes quiet. As shown in the snapshot at $70~{}\text{\textmu s}$, the main feature is that the bubble’s internal pressure is higher than the ambient value. Therefore, the bubble continues growing. By this time, it has already formed a long, conical shape.

Figure 19 shows the evolution of the velocity field. The inset image at $1.2~{}\text{\textmu s}$ shows that when the initial bubble has just formed, the high internal pressure leads to high velocity in both axial and radial directions. This phenomenon is also found in the previous case (Fig. 12, $17.4~{}\text{\textmu s}$). In the previous case, the bubble’s velocity quickly starts to decrease. In the current case, however, the velocity inside the bubble remains high. This is again because of the continuation of phase transition, as it keeps adding energy to the existing bubble. In addition, multiple vortexes are observed inside the vapor bubble as shown in the snapshots at $70~{}\text{\textmu s}$, which is related to the propagation of multiple acoustic waves emitted from the bubble’s forward tip. Therefore, the evolution of the velocity field also suggests that phase transition plays a substantial role in the bubble’s dynamics.

Compared to liquid water, the absorption of laser by water vapor is negligible. Therefore, the formation of a vapor bubble along the path of the laser beam allows laser energy to be transmitted over a longer distance. This phenomenon, shown in Fig. 17, is sometimes referred to as Moses effect, after the story of Moses parting the sea (Ventimiglia & Traxer, 2019; van Leeuwen et al., 1993). To investigate this effect more closely, we introduce four sensor points along the central axis of the laser beam, at difference distances from the fiber tip. Their coordinates (in mm) are $\bm{x}_{1}:(-0.3,0,0)$, $\bm{x}_{2}:(0.3,0,0)$, $\bm{x}_{3}:(1.1,0,0)$, and $\bm{x}_{4}:(2.7,0,0)$. The fiber tip is at $\bm{x}_{0}:(-0.5,0,0)$, as shown in Fig. 7(b). Figure 20 shows the time-history of laser radiance at the four sensor locations, as well as the fiber tip.

Before bubble nucleation ($1.2~{}\text{\textmu s}$), most of the laser energy is absorbed by a small volume of water next to the fiber tip. Although Sensor $1$ is only $0.2~{}\text{mm}$ from the fiber tip, the laser radiance at this point has already dropped to $37.6\%$ of the value at the fiber tip. The laser radiance at the other sensor locations is negligible. This means without the vapor bubble, the laser cannot reach Sensors $2$, $3$, and $4$.

After bubble nucleation, the laser radiance at all the sensor points starts to increase. At $2.4~{}\text{\textmu s}$, the bubble reaches $\bm{x}_{1}$. At this time, the laser radiance at Sensor $1$ reaches the maximum value, $98~{}\text{kW}/\text{mm}^{2}$. It is still lower than the laser radiance at the fiber tip, but this is only because the laser beam has a $9.78^{\circ}$ divergence angle.

The time-histories of laser radiance at Sensors $2$ and $3$ follow the same trend. As the vapor bubble’s forward end gets close to the sensor, laser radiance increases. The maximum value is reached when the bubble reaches the sensor. The maximum laser radiance decreases along the axial direction, due to beam divergence. Sensor $4$ is placed at $3.2~{}\text{mm}$ from the fiber tip. At $70~{}\text{\textmu s}$, the bubble’s forward tip is still more than $0.29$ mm from it. As the result, the laser radiance at Sensor $4$ remains nearly zero.

In summary, the simulation result shows that the vapor bubble essentially creates a channel that allows laser to pass through. Compared to rounded bubbles, the long, conical shape obtained in the current case (both the experiment and the simulation) can be advantageous as it provides a longer channel.

Using different laser settings, we have obtained two vapor bubbles in different shapes, namely a rounded, pear-like shape shown in Sec. 4 and a long, conical shape shown in Sec. 5. Depending on the application, one or the other may be preferred. By examining the simulation result, we find that a major difference between the two cases is that in the first case, vaporization only lasts for a short period of time, less than $1~{}\text{\textmu s}$. In the second case, vaporization continues along the laser beam direction for over $50~{}\text{\textmu s}$. It is also clear that in both cases, a newly vaporized region carries a high pressure (from the accumulated latent heat) that drives the existing bubble to expand by means of advection. Therefore, the simulation result suggests that when laser energy input is maintained in time (i.e. long-pulsed laser), the vapor bubble’s shape is influenced by two factors:

1. (1)

   the speed of bubble growth by advection, and

2. (2)

   the speed of bubble growth by phase transition.

Furthermore, the simulation result indicates that the transition between pear-shaped and elongated bubbles may be related to a competition between these two speeds. At least one of the two speeds must have changed from one case to the other. In other words, at least one of them is controllable.

Unfortunately, the fluid dynamics is highly nonlinear and multi-dimensional. It is impossible to separate the two speeds from the governing equations. In order to define and examine the two speeds analytically, we resort to a simplified model problem.

As illustrated in Fig. 21(a), we consider an initial vapor bubble of spherical shape, with radius $R_{0}$. We assume it has an internal pressure $p_{Go}$ that is higher than the ambient pressure $p_{\infty}$, which drives the bubble to expand. For this problem, the bubble dynamics can be modeled by the Rayleigh-Plesset equation (Brennen, 2014). After dropping the viscosity and surface tension terms for simplicity, we get

where $\rho_{0}$ denotes the density of the liquid phase, and $\gamma$ the specific heat ratio of the vapor phase. In this case, the bubble’s dynamics is isotropic, and driven only by advection. Therefore, we define the speed of bubble growth by advection as

where $R(t)$ is the solution of Eq. (41).

From Eq. (41), it is clear that $v_{\text{adv}}$ depends on $p_{Go}$ and $R_{0}$. If the initial bubble is created through vaporization, $p_{Go}$ is determined by the thermodynamics of water, including its latent heat of vaporization (Sec. 2.4). Therefore, it may not always be adjustable. To see the effect of $R_{0}$, we first note that if we rewrite (41) with

$R_{0}$ can be eliminated from the equation. Specifically, we have

and the solution $\widehat{R}(\tau)$ is independent of $R_{0}$. Also, substituting (43) into (42) yields

Therefore, changing the value of $R_{0}$ leads to a linear scaling (i.e. stretching or compressing) of $v_{\text{adv}}$ in time, while the peak values remain the same. If the initial bubble is created by vaporization, $R_{0}$ may be controlled by adjusting the spatial distribution of the heat source. For example, a more uniform distribution of the laser power on the source plane may lead to a larger $R_{0}$.

Next, we assume that at a time $t>0$, a uniform, parallel laser beam is activated, and it creates a bulge on the bubble surface through vaporization (Fig. 21(b)). We assume that over a short period of time, the vaporized region (i.e. the bulge) has a cylindrical shape, with a depth of $\Delta d_{\text{vap}}$. To model the continuation of vaporization, we assume that at time $t$, $T=T_{\text{vap}}$ within this cylindrical region, and $\Lambda=l$ at $x=R(t)$. This assumption is justified by the results of the simulations shown in Secs. 4 and 5.

The energy required to vaporize the forward end of the cylindrical region, i.e. $x=R(t)+\Delta d_{\text{vap}}$, can be estimated by

where $\rho_{0}$ denotes the density of the liquid phase, assumed to be a constant. The time to obtain this amount of energy from the laser beam can be estimated by

Now, we define the speed of bubble growth by phase transition as

Substituting (47) into (48), and noting that $\Lambda(R(t),t)=l$, we get

The derivative $\partial\Lambda/\partial x$ is negative, as the laser energy input decreases along the beam direction. Therefore, $v_{\text{vap}}$ is positive. It is notable that the latent heat of vaporization, $l$, is not involved in (49). Intuitively, the latent heat represents an energy threshold for phase transition to occur. Once this threshold is reached, it may not influence the speed of bubble growth by phase transition. Eq. (49) also indicates that $v_{\text{vap}}$ may be controlled by adjusting the laser’s wavelength and power, which will lead to variations in $\mu_{\alpha}$ and $L$.

We hypothesize that a race between advection and phase transition determines the morphing of the vapor bubble. In the previous subsection, we have defined speeds $v_{\text{adv}}$ and $v_{\text{vap}}$ for a simplified model problem. In real-world applications, these quantities would have to be estimated using available data. We hypothesize that at any time, if $v_{\text{adv}}$ is greater than $v_{\text{vap}}$, the bubble tends to grow spherically. If $v_{\text{adv}}$ is smaller than $v_{\text{vap}}$, the bubble tends to elongate along the laser beam direction.

This hypothesis can be tested using our simulation results. Figure 22 illustrates the method adopted here to estimate $v_{\text{adv}}$ and $v_{\text{vap}}$. Although this figure only shows a solution snapshot obtained from the TFL simulation, the same estimation method is also applied to the other simulation with a Ho:YAG laser.

Because vaporization mainly continues along the laser beam direction, we estimate the advection speed, $v_{\text{adv}}$, by measuring the speed of bubble expansion in the radial direction, outside the beam waist. Specifically, we specify a small time interval, $\Delta t_{\text{adv}}=0.2~{}\text{\textmu s}$. The radial expansion of the bubble over this time interval, $\Delta d_{\text{adv}}$ (Fig. 22), is measured at $255$ time points for the Ho:YAG simulation and $320$ time points for the TFL simulation. At each time point, $v_{\text{adv}}$ is estimated by

To estimate the phase transition speed, $v_{\text{vap}}$, we look at the forward tip of the bubble, denoted by $x_{\text{max}}(t)$ in Fig. 22. We specify a small distance, $\Delta d_{\text{vap}}=0.0015~{}\text{mm}$, along the laser beam direction. Then, we extract the simulation result of $L$ and $\Lambda$ to evaluate (47), which gives us $\Delta t_{\text{vap}}$, that is, an estimate of the time needed to extend the bubble front by $\Delta d_{\text{vap}}$. Then, we estimate $v_{\text{vap}}$ by

which is consistent with its definition in (48). Again, we calculate $\bar{v}_{\text{vap}}$ at $255$ time points for the Ho:YAG simulation and $320$ time points for the TFL simulation.

Figure 23 shows the time-histories of $\bar{v}_{\text{adv}}$ and $\bar{v}_{\text{vap}}$ obtained from the two simulations. For the Ho:YAG simulation that generated a pear-shaped bubble, $\bar{v}_{\text{adv}}$ is found to be roughly two orders of magnitude higher than $\bar{v}_{\text{vap}}$ over the entire time period shown in the figure. This is consistent with the finding that in this case, bubble expansion is mainly driven by advection, while phase transition only lasts for less than $1~{}\text{\textmu s}$. It also supports our hypothesis that if $v_{\text{adv}}$ is greater than $v_{\text{vap}}$, the bubble tends to grow spherically. For the TFL simulation that generated an elongated bubble, $\bar{v}_{\text{vap}}$ is found to be higher than $\bar{v}_{\text{adv}}$. Their difference is more than one order of magnitude in the early period of the simulation. But after $20~{}\text{\textmu s}$, the difference starts to become smaller. This is consistent with the finding that in this case, phase transition continues for a long period of time, until $53.5~{}\text{\textmu s}$. It also supports our hypothesis that if $v_{\text{adv}}$ is smaller than $v_{\text{vap}}$, the bubble tends to elongate along the laser beam direction.

In summary, the simulation results suggest that the transition between pear-shaped and elongated bubbles is determined by a race between advection and phase transition. These two speeds can be characterized by $v_{\text{adv}}$ and $v_{\text{vap}}$, which are mathematically defined for a simplified model problem. $v_{\text{adv}}$ and $v_{\text{vap}}$ depend on the laser setting and the properties of the fluid medium. Therefore, in real-world applications it may be possible to obtained a preferred bubble shape by adjusting the relevant parameters. For example, in our TFL experiment, the laser absorption coefficient and the source laser radiance are both higher than their counterparts in the Ho:YAG experiment. These changes lead to a significant increase of $v_{\text{vap}}$ by about $2$ orders of magnitude. In comparison, the variation of $v_{\text{adv}}$ is much smaller. Therefore, the changes in laser setting make $v_{\text{vap}}$ higher than $v_{\text{adv}}$. As a result, an elongated bubble is obtained.