Unified framework for laser-induced transient bubble dynamics within microchannels

Figure 1A sketch of the experimental setup with representative images of bubble dynamics from inception to collapse—one oscillation cycle. X represents the bubble size and \(t_\mathrm{osc}\) is the bubble lifetime.Dynamic bubble sizeSince the midpoint of the capillary and the bubble coincides, in our theoretical model we analyze only half the geometric domain owing to the bilateral symmetry, see Fig. 1. Thus, the length of the liquid column under analysis (\(l_\mathrm{L}\)) is half the total length of the channel used, and the half-size of the bubble (X) is determined by the position of the vapor-liquid interface from the center as illustrated in Fig. 1. In the post-processing of the images the parameter X was estimated by dividing the bubble volume with the area of the channel cross-section (see Supplementary Section 1). This approach of determining X ruled out the effects of bubble end curvatures, furthermore allowing us to estimate the liquid velocities as a function of rate of volume displaced. The liquid droplets at either end of the channel act as reservoirs, thus compensating for the liquid displaced by the vapor bubble. In addition, these droplets also ensure that the evaporation of the liquid to the ambient surrounding doesn’t deplete the liquid within the channel.After attaining the capillary diameter (\(d_\mathrm{h}\)), the vapor bubble elongates along the axial (along x) direction with the cross-section as that of the channel. Thus, for an elongated bubble, when \(X \ll l_\mathrm{L}\), the equation of motion of the liquid column within the channel is22$$\begin{aligned} l_\mathrm{L} \rho _\mathrm{L} \frac{{\rm d}^2{X}}{{\rm d}t^2} = p_\mathrm{V} (t) – p_{\infty } – {\mathscr {R}}\, \frac{{\rm d} {X}}{{\rm d} t}, \end{aligned}$$
(1)
where \(\rho _\mathrm{L}\) is the liquid density, \(p_\mathrm{V} (t)\) is the pressure inside bubble, \(p_{\infty }\) is the ambient pressure and \(\mathscr {R}\) the hydraulic resistance of the channel. Due to strong confinement in two directions, the resulting flow is quasi-1D along the longitudinal axis x, which justifies the use of a one-dimensional model. The expressions of the steady state hydraulic resistance derived using laminar flow theory for different cross-sections are38: \(32\mu _\mathrm{L} l_\mathrm{L} / a^2\) (circle); \(28.4\mu _\mathrm{L} l_\mathrm{L} / a^2\) (square) and \(\approx 12\mu _\mathrm{L} l_\mathrm{L} / [(1-0.63b/a)b^2]\) (rectangle), where \(\mu _\mathrm{L}\) is the dynamic viscosity of the liquid and a, b are the cross-section’s edge lengths with \(b<a\). For circular and square cross-sections \(d_\mathrm{h} =a\), and for a rectangular cross-section \(d_\mathrm{h} =2ab/(a+b)\). By non-dimensionalizing Eq. 1, we obtain the characteristic velocity (\(\upsilon\)) and timescale (\(\tau\)) for the channel geometry as \(\upsilon =(p_{\infty }-p_\mathrm{V} (t))/\mathscr {R}\) and \(\tau =(l_\mathrm{L} \rho _\mathrm{L} )/\mathscr {R}\), respectively.Figure 2(A) Representative bubble dynamics for different channel geometries. (B) Universal motion of bubbles within channels with different size, shape and length. The dashed line represents the developed theory, Eq. (2). The marker colors represent the hydraulic diameters (\(d_\mathrm{h}\)), the shapes represent the cross-section and the facecolor represent the lengths (L). The graphical marker symbols and colors established here are followed throughout this article. The black arrow represents the region of deviation(s) from the expected dynamics.Since \(p_\mathrm{V} (t)\) is reported to be significant only during the first \(\approx 10\%\) of the bubble’s lifetime23, we simplify the problem by dropping this term (\(p_\mathrm{V} (t)=0\)). Thus, the liquid initially gains momentum from the impulse offered by the vapor pressure. Subsequently, the bubble continues to grow due to the liquid’s inertia, even though the pressure inside the bubble is much less than the ambient pressure22. At the end of the growth phase, the liquid momentum is fully dissipated by the resistance of the channel walls, following which the ambient pressure collapses the bubble. We describe the bubble growth and collapse by solving the ordinary differential equation (ODE) from Eq. (1) and derive the following dimensionless equation for the dynamic bubble size:$$\begin{aligned} \widehat{X}_\mathrm{max} -\widehat{X} = \widehat{t} + \exp (-\widehat{t}) – 1, \end{aligned}$$
(2)
where \(\widehat{X}=X/(\upsilon \tau )\) and \(\widehat{t}=(t-t_0)/\tau\) are the dimensionless bubble size and time, respectively. The ODE is solved using the boundary condition \(X=X_\mathrm{max}\) at \(t=t_0\), where \(X_\mathrm{max}\) is the maximum bubble size. Fig. 2A shows the representative bubble dynamics from experiments for different channel cross-sections and lengths. The bubble’s maximum size, lifetime and the growth/collapse velocities vary significantly for different microchannel geometries. Figure 2B shows the experiments to follow the derived general expression for dynamic bubble size when non-dimensionalized using the characteristic scales (see Supplementary Section 1 for unscaled results).In most cases, we see a deviation in the bubble size from theory in the first 10% of the bubble’s lifetime, as expected. This can be attributed to the significant vapor pressure inside the bubble due to instantaneous phase change during formation (\(\sim\) 100 bar) which then decreases abruptly due to rapid change in the bubble volume22. Following the phase change, in the time range of \(O(\upmu s )\) the vapor pressure is reported to be still very high \(\sim 8\,\mathrm{bar}\) corresponding to a saturation temperature of \(\sim 170\,^{\circ }\)C22. After which the change in pressure and therefore the corresponding saturation temperature is rapid, 0.15 bar/\(\upmu\)s or 2 \(^{\circ }\)C/\(\upmu\)s approximately. Moreover, the force exerted by the channel wall during the radial (r, see Fig. 1) growth of the bubble appearing during this first 10% of the bubble’s lifetime will also influence its dynamics19. Thus a more sophisticated theoretical model incorporating the rapid change in pressure/temperature, change in the bubble geometry from spherical to elongated, and wall forces in radial direction is necessary to accurately match the observed bubble dynamics in experiments. We therefore will overlook the accuracy within this time regime as its effect over the calculation of the characteristic parameters such as bubble lifetime, Womersley and Reynolds numbers are not large enough, discussed in detail in the following (sub)sections.In Fig. 2A, for \(d_\mathrm{h} =300\,\upmu\)m, we see a rapid acceleration during the start of the bubble collapse (pointed using a black arrow). Similarly, in Fig. 2B, we also see a deviation from theory during the start of bubble collapse becoming more significant with increasing hydraulic diameter (\(d_\mathrm{h} \ge 200\,\upmu\)m). We attribute these to instabilities, discussed in detail under the section “Emergence of instabilities” later in this article.Bubble lifetimeWhile the dynamic size of the vapor bubble quantifies the liquid velocity or flow rate, the lifetime of the oscillating bubble governs the oscillatory time period of flow. In this work, we focus on the primary expansion and collapse of the bubble and ignore the bubble rebound caused by high pressures and temperatures39. Thus, the lifetime of the bubble is the time taken for one oscillation, \(t_\mathrm{osc}\). In Eq. (2), by substituting \(\widehat{X}=0\), we obtain the bubble’s time of formation and collapse. The analytical prediction of the bubble’s lifetime is$$\begin{aligned} \widehat{t}_\mathrm{osc} = t_\mathrm{osc} / \tau = W _0 (-e^{-\xi }) – W _{-1} (-e^{-\xi }), \end{aligned}$$
(3)
where W\(_0\) is the principal branch of the Lambert W function and W\(_{-1}\) its only other real branch, and \(\xi =1+\widehat{X}_\mathrm{max}\). The dimensionless times \(W _{-1} (-e^{-\xi }) + \xi\) and \(W _{0} (-e^{-\xi }) + \xi\) correspond to the time span of bubble expansion and collapse, respectively. We note that the dimensionless time in Figure 2B is negative due to the choice of time zero at the maximum bubble size.Figure 3(A) Experimentally determined dimensionless lifetime (\(\widehat{t}_\mathrm{osc}\)) against maximum size (\(\widehat{X}_\mathrm{max}\)) of the bubble (markers) and comparison with the theoretical prediction from Eq. (3) (dashed line). (B) General linear relation between the non-dimensionalized experimentally determined maximum bubble size (\(X_\mathrm{max}\)) and laser energy absorbed by the liquid (\(E_\mathrm{abs}\)). The size and energy are non-dimensionalized using the channel hydraulic diameter (\(d_\mathrm{h}\)) and unitary kinetic energy (\(E_\mathrm{ke}\))-calculated using the channel geometry and liquid properties. The inset illustrates the intercept in the plot which corresponds to the threshold laser energy for bubble formation (\(E_\mathrm{th}\)). The coefficient of determination for the fit line is \(R^2=0.928\).Figure 3A illustrates the relation between the dimensionless bubble lifetime and maximum bubble size. We report a good agreement between the experiments and the theoretical prediction from Eq. (3). See Supplementary Section 1 for unscaled data presented in Fig. 3A. There is a decrease in the bubble lifetime with an increase in the hydraulic diameter for a fixed channel length and \(X_\mathrm{max}\). This is caused by the decrease in hydraulic resistance with increasing hydraulic diameter, resulting in faster bubble dynamics. While the exact solution to the Lambert W function in Eq. (3) accurately captures the experiments, an approximation of this equation is discussed in Supplementary Section 1 to help solve the equation with more commonly used mathematical functions. The approximation however does decrease the accuracy of the results.The above sections discuss the theory for dynamic bubble size (X(t), Eq. (2)) and lifetime (\(t_\mathrm{osc}\), Eq. (3)). The proposed theory demands \(X_\mathrm{max}\) as the only necessary parameter to compare with experiments. Thus it is necessary to have an estimate of \(X_\mathrm{max}\) as a function of the laser energy supplied—a control parameter in experiments.Energy balanceFigure 3B represents a general scaling law using empirical correlation for the maximum bubble size (\(X_\mathrm{max}\)) against the absorbed laser energy, \(E_\mathrm{abs}\). The parameters \(X_\mathrm{max}\) and \(E_\mathrm{abs}\) are non-dimensionalized with respect to the hydraulic diameter and kinetic energy of the liquid with unitary velocity, respectively. The liquid kinetic energy with unitary velocity (\(E_\mathrm{ke}\)) is calculated as$$\begin{aligned} E_\mathrm{ke} = A l_\mathrm{L} \rho _\mathrm{L} , \end{aligned}$$
(4)
where A is the cross-sectional area of the channel. This approach to determine \(X_\mathrm{max}\) via empirical correlation avoids the need for a more sophisticated numerical model with energy balances and phase transition39.For a certain cross-sectional shape, the energy required to achieve a certain \(X_\mathrm{max}\) increases with the hydraulic diameter and length of the channel (see Supplementary Section 1 for unscaled results of Fig. 3B). In Fig. 3B we see the \(E_\mathrm{abs}\) to be at least two orders of magnitude higher than the \(E_\mathrm{ke}\). This observation is in agreement with the numerical simulations by Sun et al.22, where most of the absorbed laser energy heats up the liquid directly surrounding the bubble, rather than transforming into the kinetic energy of the liquid.In the experiments, the amount of laser energy absorbed by the liquid depends on the liquid properties and channel geometry. Based on the Beer-Lambert law40, the absorption coefficient of the liquid and the distance the light travels through the liquid is used to estimate the amount of light absorbed, \(E_\mathrm{abs}\). Furthermore, the channel’s wall thickness, material, cross-section shape and dimension can influence the laser energy available for the liquid to absorb. Thus the energy absorbed by the liquid in experiments was calibrated by measuring the difference of energy transmitted with the microchannel containing water and water-dye mixture (working fluid). Since water has very low absorption coefficient (\(<0.001\,\mathrm{cm} ^{-1}\)) for the laser wavelength used (532 nm)41, we attribute the measured energy difference to the laser energy absorbed by the dye.For rectangular channels, the longer edge length was used for supporting the channel over the stage (see Fig. 1), while the laser light was passed along the shorter edge. The details on the energy absorption measurement technique can be found in Supplementary Section 2 . The estimated absorption coefficient of liquid, \(\alpha\), is \(86.54\,\mathrm{cm} ^{-1}\). This estimated value is different from measurements from the spectrometer (DR6000, Hach), 173 cm\(^{-1}\), which employs an unfocused and continuous light source. We attribute this difference in absorption coefficient to the non-linear absorbance (saturable absorption) of the liquid due to high energy intensities (\(\sim \mathrm{GW\, cm ^{-2}}\)) as the laser is focused and has short laser pulse duration (4 ns).The fit in Fig. 3B has an intercept for \(E_\mathrm{abs}\), implying there exists a threshold absorbed energy (\(E_\mathrm{th}\)) only above which a vapor bubble forms. For a bubble to form, theoretically we use the spinodal temperature of water as the necessary condition42—the temperature at which water explosively turns into vapor. Thus, a sensible heat corresponding to \(\Delta T = 279\,\mathrm{K}\) (at 1 atm and with respect to room temperature of 298 K) and a latent heat proportional to enthalpy of vaporization (\(H_\mathrm{L}\)) is minimum required at the laser spot for bubble formation. Supplementary Section 2 provides the derivation for analytical expression of \(E_\mathrm{th}\) based on the energy balance for liquid. The threshold absorbed laser energy is,$$\begin{aligned} E_\mathrm{th} = (1-10^{-\alpha d_\mathrm{L} })\frac{\pi w_0^2 \rho _\mathrm{L} (c_p\Delta T + H_\mathrm{L} )}{\alpha \log (10) 10^{-\alpha (d_\mathrm{L} /2)}}. \end{aligned}$$
(5)
Where \(w_0\) is the laser spot radius, \(c_p\) is the specific heat and \(d_\mathrm{L}\) is the distance the light will travel through the liquid. \(d_\mathrm{L}\) is equal to the hydraulic diameter (\(d_\mathrm{h}\)) for circular and square channels. For rectangles it is the shorter edge length of the cross-section. The parameter values for the calculation of Eq. (5) are presented in Table 2.Table 2 Parameter values for theoretical calculation of the threshold energy for bubble formation using Eq. (5).Figure 4The threshold laser energy absorbed for bubble formation estimated from experiments (\(E_\mathrm{th,exp}\)) against theory (\(E_\mathrm{th,theory}\)) presented in Eq. (5).In Eq. (5), the other necessary parameter for calculation is the laser spot diameter, \(2w_0\). Theoretically, the laser spot diameter is calculated using the expression \(2w_0 = 4M^2\lambda f / (\pi D_\mathrm{L} ) = 1.69\,\upmu\)m. Where \(\lambda =532\,\mathrm{nm}\) is the laser wavelength, \(f=10\,\mathrm{mm}\) is the objective lens focal length, \(D_\mathrm{L} =4\,\mathrm{mm}\) is the laser beam diameter and the beam quality parameter \(M^2=1\) (assuming a perfect Gaussian profile). The depth of field (DOF) of the focused laser beam is \(2\pi w_0^2 / (M^2 \lambda ) = 8.47\,\upmu\)m. However, these theoretical calculation for spot size and DOF can be different in experiments due to the optical aberrations caused by the channel walls43. From Fig. 4, we see a good agreement between experiments (\(E_\mathrm{th,exp}\)) and theory (\(E_\mathrm{th,theory}\)) presented in Eq. (5) for the threshold laser energy absorbed by the liquid for bubble formation. The laser spot diameter used in the analytical expression is \(9\,\upmu\)m, which is estimated based on the best fit of the theory to the experiments. This value of the laser spot diameter is same order of magnitude as the theoretically calculated value, thus emphasizing its validity. Furthermore, in Fig. 4 we see a deviation in the threshold laser energy absorbed. We attribute the deviation to the laser aberrations due to channel wall curvature and thickness that can influence the laser spot diameter and therefore the threshold laser energy. In addition, the r position of the laser spot can also be affected by the confinement wall curvature and thickness, resulting in a larger \(E_\mathrm{th}\) as the laser spot moves towards \(r=d_\mathrm{L} /2\). The channel specifications used in this work are summarized in Table 1.The above analysis thus provides us with an expression for bubble formation as a function of the energy absorbed by the liquid – a parameter that can be measured in experiments. The estimation of the laser energy threshold is one of the key design parameters necessary to choose the range of the laser energy required based on the channel geometry, liquid properties and objective lens specification.Emergence of instabilitiesIn Fig. 5A, for \(d_\mathrm{h} \ge 200\,\upmu \mathrm{m}\) we see the bubble dynamics deviating from an expected bell curve like profile. Interestingly, the deviation occurs before \(X_\mathrm{max}\) for \(d_\mathrm{h} =200\,\upmu \mathrm{m}\) and after \(X_\mathrm{max}\) for \(d_\mathrm{h} =300\,\upmu \mathrm{m}\). Furthermore, Fig. 5B shows that for \(d_\mathrm{h} =200\,\upmu \mathrm{m}\) the deviation disappears for larger \(X_\mathrm{max}\). We hypothesise the observed deviations to the boundary layer instabilities near channel walls. These instabilities emerge due to the oscillating nature of the flow and unlike unidirectional flows it cannot be characterized only using the Reynolds number. Thus, the flow instability in oscillating channel flow is characterized using a peak oscillatory Reynolds number \(Re_\mathrm{osc} = \rho _\mathrm{L} U_\mathrm{max} d_\mathrm{h} /\mu _\mathrm{L}\) and a Womersley number \(\textit{Wo} = d_\mathrm{h} /2\sqrt{\omega \rho _\mathrm{L} / \mu _\mathrm{L} }\)44,45. \(U_\mathrm{max}\) is the maximum mean flow velocity during the bubble growth, which on average occurs in the middle of the duration of growth, and \(\omega = 2\pi /t_\mathrm{osc}\) is the angular frequency of oscillation. Figure 5C shows these dimensionless numbers for the experiments performed in this work.Figure 5(A,B) Representative dynamic bubble size curves illustrating the emergence of instabilities. The zones of the instabilities are highlighted using a shaded rectangular area. The arrows represent if the instabilities occur before or after \(X_\mathrm{max}\). (A) Illustrates the experimental data for different \(d_\mathrm{h}\) with similar oscillation time. The instabilities emerge with increasing \(d_\mathrm{h}\). (B) Illustrates the data for \(d_\mathrm{h} =200\,\upmu\)m with increasing laser energies. The instabilities disappear with increasing \(E_\mathrm{abs}\). (C) Flow stability diagram with the transition line at \(\textit{Wo}=7\)34. The markers represent the experiments and the lines represent the analytical estimate. The numbers correspond to the channel hydraulic diameters (in \(\upmu\)m) with the dashed and solid lines representing the channel lengths L = 25 and 50 mm, respectively. (D) The dimensionless convective timescale against the \(L/d_\mathrm{h}\) aspect ratio. The partition line is a linear relation between the x and y axes with \(45\times 10^{-6}\) as the slope and the origin as the intercept.A simple approach to estimate the flow velocity, U, using theory can be formulated by using the analytical expression for \(t_\mathrm{osc}\) (Eq. 3) and differentiating Eq. (2) with respect to time. The resulting dimensionless flow velocity is$$\begin{aligned} U/v = -1 + \exp (-\widehat{t}). \end{aligned}$$
(6)
By substituting half the bubble growth time \([W _{-1} (-e^{-\xi }) + \xi ]/2\) for \(\widehat{t}\), we calculate a corresponding maximum mean flow velocity \(U_\mathrm{max}\) and \(Re_\mathrm{osc}\).The lines in Fig. 5C represent the analytical prediction of the dimensionless flow parameters for varying \(X_\mathrm{max}\) in all channel geometries. We report a good agreement between the experiments and the analytical estimate. The horizontal line in the figure, \(\textit{Wo}=7\), represents the empirically observed laminar to distorted laminar transition value from literature for a perfectly oscillating flow with sinusoidal pressure gradient34. To predict the transition, the onset and growth of instability is determined using the parameters: (i) the position of the point of inflection, and (ii) the convective timescale, respectively. Das et al.46, using pulsatile flow theory with a sinusoidal mean velocity profile8, showed the dimensionless position of the point of inflection to be independent of the Wo above a critical Stokes parameter, \(d_\mathrm{h} /(2\delta _\mathrm{st} ) = \textit{Wo}/\sqrt{2} \approx 3.6\), where \(\delta _\mathrm{st} = \sqrt{2\mu _\mathrm{L} /(\rho _\mathrm{L} \omega )}\) is the Stokes layer thickness. This theoretical critical value of \(\textit{Wo}\approx 5\) is below the empirically observed critical transition value of \(\textit{Wo}=7\). Thus, the aforementioned analysis of instability using \(\delta _\mathrm{st}\) demonstrates the critical \(\textit{Wo}\) in flow transition to be independent of the position of the point of inflection. However, the source of instability can still be attributed to the inflection of the velocity profiles near channel walls occurring due to the flow reversal (boundary layer separation) during the deceleration phase of the liquid. To further emphasize the observed deviation in Fig. 5A as a cause of instability and not due to assumptions underlying the theory under consideration, we compare it to laminar flow due to an oscillating pressure gradient. Just before the bubble attains its maximum size, detachment of the boundary layer at the confining walls has been reported for flows between parallel plates based on numerical analysis47. The flow detachment exists due to the mismatch between the direction of the pressure gradient (\(\partial p/\partial x > 0\)) and fluid flow (\(U>0\)). In accordance, our estimated flow profiles using theory with oscillating pressure gradient for microchannels48 also predict a flow reversal occurring close to the channel walls (see Supplementary Section 3), similar to others49. However, the flow reversal near walls should decelerate the flow as the bubble begins to collapse – opposite to what is observed in experiments (in Fig. 5A, \(d_\mathrm{h} =300\,\upmu\)m). This is a consequence of the instabilities causing flow distortion.While the inflection in velocity profile near channel walls can lead to the onset of instabilities, the growth and therefore the time of emergence of instabilities is governed by the convective timescale. Fig 5D shows the dimensionless convective timescale from experiments, \(\delta _\mathrm{st} /(U_\mathrm{max} \tau )\), for different \(L/d_\mathrm{h}\) ratios. Small characteristic timescales of the channels (\(\tau\)) compared to convective timescale (\(\delta _\mathrm{st} /U_\mathrm{max}\)) would correspond to large resistance to flow offered by the channel walls as \(\tau \propto 1/\mathscr {R}\). Large resistance (\(\mathscr {R}\)) in other words would mean a viscosity-dominated flow that effectively would dampen any perturbations/instabilities. However, as \(\tau\) increases the instability would develop during the deceleration phase and prevail in the acceleration phase46. This argument on instability evolution explains the observed orientation of flow distortion relative to \(X_\mathrm{max}\) in Fig 5A, when analyzed using the partition line depicted in Fig 5D. The partition line is adapted based on the experiments from this work which emphasizes the time of emergence of instabilities. As a data point approaches the partition line, the instabilities die out due to the motion-less state (zero velocities) at \(X_\mathrm{max}\), while below the partition line the instabilities sustain and prevail due to larger kinetic energies (\(\propto U_\mathrm{max}\)).In summary, the observed deviations in the dynamic bubble size from a bell curve like profile can be attributed to the disruption of the momentum boundary layer near channel walls36 that would alter the channel hydraulic resistance and hence the mean velocities. However, these distorted laminar states are transient and therefore revert (decay) to laminar flow over time36,44. Consequently, there exists a momentary deviation in dynamic bubble size as seen from the representative images in Fig. 5A.