Design Exploration of Initial Bubbles for Higher-Speed Laser-Induced Microjets

Laser-induced microjets are advantageous for yielding high-speed and convergent microjets shape. However, owing to the need for high-power laser devices, reducing the required laser power is a prerequisite for convenience in practical use. This study aimed to optimize initial multiple bubbles design and find the dominant parameters for efficient microjet generation. For this purpose, the bubble design was optimized using an evolutionary algorithm—the covariance matrix adaptation evolutionary strategy. In addition, the dataset of solutions obtained through the optimization process was analyzed via principal component analysis (PCA). The optimization and analysis revealed that a higher microjet velocity could be obtained when a single bubble was generated near the tube wall, and an optimal total bubble volume was observed corresponding to this optimization. In addition, this approximate optimal solution we obtained was confirmed to reduce 88% of total initial energy compared to that of the typical solution. The result shows possibility to improve the efficiency of microjet generation and reducing the laser power.


I. INTRODUCTION
As proposed by Tagawa et al., the laser-induced microjet can achieve supersonic speeds [1]. Moreover, traditional microjet-generation methods impact the drop on a liquid interface [2], rupturing bubbles at liquid interfaces [3], and generate an impulsive force at the bottom of a liquid-filled tube [4]. Compared to those formed via the traditional methods, laser-induced bubbles can generate faster microjets, and this feature renders them desirable for engineering applications, such as inkjet printing [5], nanodroplet generation [6], and needle-free drug injection [7], [8], [9]. Fig. 1 presents an illustration of the laser-induced microjet generation. The generation process is as follows: a pulsed laser is focused inside a liquid-filled microtubule, thereby producing high-pressure bubbles. Subsequently, pressure waves from these high-pressure bubbles propagate through the liquid to the air-liquid interface ( Fig. 1 (a)), and the interaction with the interface generates a convergent jet at the concave air-liquid interface (refer to Fig. 1 (b)).
A high-power laser system is required for the generation of supersonic microjets (e.g., a previous study [7] employed The associate editor coordinating the review of this manuscript and approving it for publication was Lei Wang. a 100-mJ Nd:YAG laser system). In particular, improving the efficiency of jet generation and reducing the required laser power are prerequisites for enhancing the feasibility of such systems. In this regard, several studies have focused on laser-induced microjets. Peter et al. characterized the effect of the contact angle between the air and liquid interface and wall on the microjet velocity [10]. Moreover, Kyriazis et al. investigated the effect of the air-liquid interface's shape on the microjet velocity [11]. In particular, they reported that initial bubble parameters, such as bubble volume, number of bubbles, and position, strongly impacted the jet velocity even with the same laser power. Furthermore, Tagawa et al. experimentally investigated the effect of the initial bubble volume by changing the laser-focus position, demonstrating that the jet velocity increases with the initial bubble volume [1]. In addition, they performed experiments by changing the magnification of the laser focusing lens. Based on an experimental investigation of the pressure impulse emanating from the bubbles, they revealed that the bubble position and number significantly affect the jet velocity [12]. Reportedly, the pressure impulse is strongly correlated with the jet velocity [10], [13]. From this perspective, the design exploration of bubble generation is a significant factor for maximizing jet velocity.
In this study, we aimed to optimize initial multiple bubble design and investigate the dominant parameters of efficient microjet generation numerically. The bubble design is optimized assuming that the total energy given is constant in all cases. However, the jet behavior exhibits a strong nonlinearity against the design parameters because it is determined by multiple phenomena such as the pressure wave generation and its propagation and interaction with the channel wall; therefore, design exploration is a challenging task. To address this challenge, we adopted the covariance matrix adaptation evolutionary strategy (CMA-ES) [14] for the design optimization of the initial bubbles. Furthermore, to extract the dominant parameters, the dataset of solutions obtained in the optimization process was analyzed using principal component analysis (PCA). In the optimization, a compressible multiphase-flow solver based on a five-equation model [15] was adopted to evaluate the solution. To reduce the simulation cost, the simulation was conducted in a two-dimensional space. In addition, because microjet tracking is computationally expensive, the impulse of the pressure wave was used as the maximization target. The essential physical phenomena are expected to be included in the two-dimensional simulation, and the results reveal that the pressure wave impulse is strongly correlated with the microjet velocity reported in previous numerical and experimental studies, as mentioned in this section [10], [13]. Therefore, the research approach adopted in this study was suitable for our research purposes. Owing to the experimental difficulty in precisely observing and controlling initial design of bubbles, it is challenging to experimentally confirm our numerical results. However, our results can provide some insights for efficiently generating laser-induced microjets.

A. OPTIMIZATION PROBLEM
We define the optimization problem to be solved as follows: where I (b) denotes the impulse which is related to the microjet velocity [10], [13], [16], b represents the bubble design variable, y i ∈ [−7R, 7R] depicts the position of the four circles, r i ∈ [0.4R, R] indicates the circle radius, and R = 35µm symbolizes the reference length. The initial bubbles are set based on variable b in the optimization process via evolutionary computation. The energy E is uniformly distributed in all the initial bubbles as the internal energy. The total energy E of initial bubbles was constant in all solutions and set as 36.1 mJ. To obtain the impulse, I (b), a multiphase-flow simulation is required. The simulation model and method are described in Section II-C.

B. OPTIMIZATION ALGOLITHM
CMA-ES, an evolutionary algorithm (EA), is one of the most effective and successfully used black-box optimization algorithms [14], [17]. EA is an optimization algorithm based on biological evolution, namely iterative generating solutions (recombination and mutation) and selection. Among these algorithms, CMA-ES features solutions that are generated in the normal distribution defined by the covariance matrix, which is learned self-adaptively from the generated solutions. Owing to self-adaptive learning, the optimization performance of the CMA-ES exhibits a low level of dependence on the properties of the individual optimization problem. In addition, this algorithm functions effectively for multimodality, nonseparability, and discontinuity of optimization problems [17], [18]. Owing to these advantages, CMA-ES has been leveraged in various engineering optimization problems [19], [20], [21], [22], [23], [24]. Considering its advantages, the CMA-ES is expected to be an effective strategy for solving our targeted problem. In particular, considering that the objective function may be multimodal owing to the strong nonlinearity of the fluid dynamics, this algorithm's effectiveness for multimodal objective functions is of great significance.
In this study, we implemented the CMA-ES through pymoo (version 0.4.2.2), which is a single-objective and multi-objective optimization framework [25]. The parameters of CMA-ES implemented were in accordance with the recommended parameters of pymoo.

C. MULTIPHASE FLOW SIMULATION
The details of the simulation utilized for the evaluation of impulse I (b) are described in this section. Fig. 3 illustrates the simulation settings. We evaluated the impulse on the measurement surface where p is the pressure, and S is the area of the measurement surface. The pressure p at each time is obtained based on the simulation result. We considered compressible non-viscous flow without surface tension and phase transition. The governing equation is a five-equation model [15]; it a compressible-multiphase fluid equation, which is derived from the non-equilibrium Baer-Nunziato model [26], assuming velocity and mechanical equilibrium. Particularly, it is a diffuse-interface method and has advantages over the sharp interface method because it does not require geometric reconstruction, remeshing or special thermodynamic treatments [27]. Therefore, this model has been leveraged for solving several practical problems [28], [29], [30], [31], [32]. The five-equation model can be written as follows: ∂ρu ∂t where α b ∈ [0, 1] denotes the volume fraction, ρ depicts the density, u = (u, v, w) represents the velocity vector (Cartesian coordinates), and e indicates the total specific energy. The subscript i ∈ b, l indicates the quantity of the bubbles or liquids. In this study, the fluid flow was expected to be dominated by the strong pressure waves and rapid bubble interface motion coming from the high-pressure laser induced bubbles. Therefore, we disregarded the viscosity, phase transition, and surface tension. Zhao et al., in developing their computational method for simulating thermal cavitation dynamics induced by long-pulsed laser, also neglected viscosity and surface tension [33]. Experimental studies of bubble collapse near walls by Zhai et al. also did not consider the effects of viscosity or surface tension [34]. Peter et al. simulated laserinduced microjets, neglecting viscosity and phase transition, and obtained the result that the microjet velocity was consistent with experimental results [10]. Reese et al. simulated the dynamics of cavitation bubbles, neglecting the effects of phase change as in our study [35]. To close the system, a stiffened gas equation of state (SG EoS) was adopted.
where γ indicates the specific heat ratio, p ref represents the reference pressure, and ϵ i denotes the internal energy. By assuming mechanical equilibrium, (7) can be rephrased as: The parameters of the SG EoS are established as follows: All physical quantities are nondimensionalized as follows: The five-equation model was discretized using the finite volume method. Numerical fluxes are computed using the Harten-Lax-van Leer-contact approximate Riemann solver [28], which is frequently utilized in five-equation modeling [27], [29], [30], [31], [32]. For the reconstruction, we leveraged the third-order monotonic upstream-centered scheme for conservation laws [36] featuring a minmod limiter [37] (MUSCL3) and ρ-THINC [31]. ρ-THINC is a variation of the tangent of the hyperbola interface-capturing method (THINC) [38], [39] designed to sharpen the volume fraction and density at the interfaces. This coupling reconstruction method (ρ-THINC-MUSCL3) enables the capture of sharp interfaces and stabilizes the computation. The second-stage strong stability preserving Runge-Kutta is used for time marching.
The simulation domain is set as [0, 30] × [−7, 7], with the symmetric condition on the left side, outflow condition on the right side, and slip-wall conditions elsewhere. The outflow condition is implemented via zero-order extrapolation and is set sufficiently far from the measurement surface. We utilized isotropic cells for a domain x ∈ [0, 10], and the cell size in x ∈ [10,30] is stretched toward the outflow boundary with a constant stretching ratio of 1.2. Moreover, the cell width is set to 0.04 in x ∈ [0, 10]. The initial conditions of the liquid are established as follows: (ρ l , u, p) = (10 3 , 0, 0.7143), (17) and the initial conditions of the bubbles are set as where V denotes the bubble volume, k symbolizes the bubble number, and the bubble pressure, p b , is computed using (7), (18) and (19). Generally, the volume fraction of the bubble interface is smoothed via a hyperbolic tangent function for computational stability [27], [31], [40]. However, the hyperbolic-tangent-based smoothing function for a spherical bubble cannot be applied because nonspherical bubbles are generated during optimization. To account for this phenomenon, we leveraged a Gaussian kernel for smoothing, as follows: where n denotes the neighborhood cell number, and the variable h is set to be 1.2 times the cell spacing.
The time interval is set as t = 1.0 × 10 −4 based on the Courant-Friedrichs-Lewy condition, and the Courant number corresponds to 0.2.
Validation of our solver and the grid dependence were confirmed and investigated in previous studies [16].

D. VISUALIZATION
A dimensionality reduction algorithm, PCA, was employed for visualizing the datasets of the solutions obtained through the optimization. The first principal component, X 1 , is defined to maximize the data variance. Furthermore, principal component X n is defined to maximize the variance under the constraint of orthogonal to X n−1 , and the dimension of the datasets can be reduced by projecting the solution to a lower-dimensional space defined by our selected principal components.
Two datasets were visualized to investigate the relationship between the bubble design variables and the impulse. The first dataset simply comprises bubble design variables and the impulseb = (y 1 , r 1 , . . . c, y 4 , r 4 , I ). The second dataset contains the initial pressure in each cell and impulsep = (p 1 , p 2 , . . . c, p m , I ), where m denotes the cell number in the region [0, 1.0] × [−7, 7] which includes the initial bubbles. Considering that the pressure distribution is show the bubble position,p can be assumed to represent an adequate solution expression.  Fig. 4 plots the optimization history. Observably, the baseline solution is a bubble with unit radius placed at y =0 in the initial condition, which is the typical solution observed in the experiment. The CMA-ES successfully explores a solution that is superior to the baseline solution, as indicated in Fig. 4. Fig. 5 depicts the transition of the solution during the optimization process. As the optimization progresses, the bubbles are clearly placed near the wall. We obtained the approximate optimal solution by placing a hemispherical bubble on the wall side, as illustrated in Fig. 5.      from the initial pressure. Moreover, a significant weight is distributed near the wall in X2, and this result supports the assumption that the initial bubbles near the wall strongly influence the impulse. Furthermore, as the total volume and initial pressure of the bubbles are related through (6) and (18), X2 is expected to be related to the total bubble volume. As indicated in Fig. 8, the bubble volume evidently varies approximately along the X2 axis. This result conclusively indicates the existence of an optimal bubble volume for maximizing the impulse.

A. OPTIMIZATION RESULT
Based on the results of previous investigations, we conclude that the initial bubble parameters, especially the positions, number of bubbles, and total bubble volume, are dominant factors for the impulse. Finally, for further elaboration, a scatter plot of the solutions against these parameters is displayed in Fig. 10; the vertical axis represents the impulse. The plot results clearly indicate that a single bubble with an optimum volume near the wall is effective for generating a higher impulse.
Although it is difficult to experimentally verify our approximate optimal solution, these results imply that adjusting the laser-focusing position and spot size are effective strategies for laser-induced microjet generation. Inducing bubbles on the wall surface with gold nanoparticles may also be an effective strategy.

C. ENERGY REDUCTION RATE
Finally, we show the energy E R required for the approximate optimal solution to generate the impulse equivalent to that by the baseline solution and the energy reduction rate achieved in the approximate optimal solution is clearly demonstrated. E R is determined using the bisection method via the following equation: where I baseline is the impulse of the baseline solution and I * (E) is an impulse of the approximate optimal solution at the energy E. We perform the computations until the condition |f (E)| < ϵ is achieved, where the threshold ϵ is set to 10 −7 , and obtain the result E R = 4.332 mJ. This result shows 88% reduction of required laser power compared to the baseline solution. Although this reduction rate was achieved under specific conditions and settings in our numerical simulation, it implies that placement of an initial bubble near the wall with appropriate volume has potential to improve the microjet generation efficiency.

IV. CONCLUSION
In this study, we optimized bubble design under the condition that the total energy given is constant in all cases. In addition, the solutions obtained through the optimization process were visualized through PCA. The following conclusions were drawn: placing the bubbles on the wall, reducing the bubbles, and setting optimal bubble volume were important parameter for microjet generation. In addition, we obtained the approximate optimal solution, which is a single bubble placed near the wall, and this phenomenon was accompanied by an optimal total bubble volume. This approximate optimal solution was confirmed to reduce 88% of total initial energy compared to that of the typical solution.