The collapse of microbubbles near a fiber is an example often encountered in water treatment situations and cavitation fibrillation processes. However, due to the broken symmetry conditions, this process has not been studied in detail experimentally or numerically, making it difficult to precisely measure or simulate the rapid bubble evolution during collapse. In this work, we present a novel experimental method, allowing for precisely repeatable cavitation events observation, combined with numerical simulations offering insight into pressure and velocity fields distribution developments in time. Both experimental and numerical works focused on small distances between the bubble and the fiber, where the physical interaction between subjects is the strongest. Four different bubble offsets were considered within the scope of this work, and very good agreement of numerical simulations with experiments was found in all cases. Two modes of bubble collapse were identified, leading to mushroom-shaped bubbles at positions closest to the fiber and a pear-shaped bubble at the farthest position. It is noteworthy that in all four cases, a planar jet formation toward the fiber was observed. The formed jet initially assumes an elongated shape, whereas its stability depends on the mode of bubble collapse. Numerical analysis of the planar jet as the defining feature of the collapse defined lower bounds for the actual values of peak jet velocities, ranging between 250 and 330 m/s, and the resulting impact pressures, which range from 100to 500 MPa.

## I. INTRODUCTION

Cavitation as a physical phenomenon describes appearance of small cavities, depictable as vaporous and gaseous bubbles, within the liquid due to a local pressure drop from acoustic sound waves (acoustic cavitation, also referred to as sonication) or a local pressure decrease (hydrodynamic cavitation). The cavitating vaporous structures are unstable, and they often collapse violently. The cavities collapse rapidly once the local static pressure increases above the vapor pressure, giving rise to localized large hydrodynamic forces in the form of pressure waves of several 100 MPa, high velocities >100 m/s, and high shear flow and momentary temperature peaks of several 1000 K, causing water dissociation and the formation of hydroxy radicals (•OH).^{1} While some of the physical mechanisms, which are known to accompany cavitation, are well described, little is known about which, or a combination of which, mechanisms are important for a specific consequence of cavitation. Hence, engineers need to deal with nuisances, such as the fall in turbomachinery efficiency, production of noise, vibrations, and wall erosion.^{2} On the other hand, there is a great potential to utilize cavitation in various important applications in the fields of biology,^{1} chemistry,^{3} medicine,^{4} and environmental protection.^{5}

For more systematic studies, isolated single bubbles are required. The liquid breakdown is usually induced by a focused laser light pulse; however, other approaches also exist.^{6} Yet, even on the simplest level, a single cavitation bubble can give a wide ensemble of phenomena, such as splitting,^{7} shock wave emission,^{8} splashing,^{9} luminescence,^{10} and jetting.^{11} The latter—jetting—being one of the more interesting ones for both fundamental research and applications.^{12} If the bubble is exposed to an asymmetrical pressure field, a jet will normally occur as it collapses. This can occur in the case of various anisotropy drivers, such as nearby boundaries, other bubbles, stationary potential flow, liquid–liquid interface^{13} and even the presence of a gravitational field.^{11}

Single bubble collapse in the vicinity of various boundaries and their underlying mechanisms have been a subject of extensive research in the past due to a wide range of practical applications in a range of disciplines.^{14} The simplest manifestation of the anisotropy is the well-researched “infinite” flat surface, but even there new phenomena, such as very fast jets, are being discovered in the last years.^{15–18} On the other hand, researchers are focusing on other geometries, such as flat surfaces with defined size (comparable to bubble size) and spheres. In all cases, struggling to preserve a symmetrical environment, which eases experimentation and above all simulations.

Almost concurrently, two independent studies were published discussing bubble dynamics on the top of a cylindrical rod^{19,20} based on the experimental and numerical approach. The results showed that the bubble dynamics collapsing near the cylindrical rod differed significantly from the single bubble collapse near solid flat plates. When the bubble was generated close to the flat top of the cylinder and its maximum radius exceeded one of the flat top surfaces, it collapsed in the form of a “mushroom” with a footing on the cylinder, a long stem, and a hat-like cap typical for a mushroom head. In cases of proximate bubble collapses also, fast thin jets (>1000 m/s) could be observed in simulations.

Another geometrical possibility is the collapse of a bubble near a spherical particle. Previous studies report that bubbles will assume a pear-like or mushroom-like shape during collapse, depending on the distance from the sphere and their size ratio.^{21,22} A recent work by the present authors^{23} reports an elongated shape during collapse which causes the bubble to split, resulting in a fast and thin annular jet (peak velocity approximately 2100 m/s) toward the sphere. They show that the jet toward the sphere forms only when a sphere is larger than the bubble.

A special case of anisotropy is a flat (infinite), but elastic surface, which was studied by Brujan *et al.*^{24} Similar mushroom-shaped collapses are observed, but in addition, significant tendency of bubbles to split and develop into two opposing jets (toward and away of the boundary) could be determined. Similar results are also obtained by Zevnik and Dular^{25,26} for the case of elastic spherical surfaces. One could determine that when the bubble is not in contact with the structure at its maximal size, it will retain almost spherical shape also during the collapse phase. On the contrary, for cases when the bubble is similar or slightly larger than the sphere and initially almost in contact with it, a thin jet away from the cell will form and the bubble will also split.

Cavitation bubbles of millimeter and sub-millimeter size have been investigated for more than 100 years. Yet, recent research points to the fact that the dynamics of microscopic cavitation bubbles can significantly differ from the dynamics on larger scales.^{23,27} However, almost no experimental investigations of microscopic bubbles in a precisely controlled environment were performed to the date—due to difficulties in both time domain and spatial domain resolution of visualization techniques.^{28}

In the present work, we go past the comfort of a symmetrical environment and generate a microscopic cavitation bubble on top of a thin fiber (cylinder) of a comparable size. We show that breaking the symmetry causes formation of a never before encountered planar jet. From an applied perspective, understanding the specifics of bubble dynamics near elongated structures at a microscale is important to progress the work on exploitation of cavitation for water treatment. There bubbles are used to eradicate bacteria, typical bacteria (*Escherichia coli* or *Legionella pneumophila*) being a cylinder several *μ*m in diameter and several 10 *μ*m long. Although the methodology provides good results, the bubble–bacteria interaction mechanisms are still only partially understood, slowing down the optimization and the development of the technology. The present work is also significant for the interpretation of recent experimental studies of cavitation fibrillation of the cellulose fiber,^{29} which could, if carefully liberated from biomass, become a promising building block for a bio-based society. There we postulate that cavitation acts on the cellulose by microstreaming, which can stretch or bend the cellulose fiber and fibril.

First, we present Sec. II, divided into the experimental setup and imaging sections, followed by a separate section for imaging, introducing the novel technique applied in this work. The complementary numerical methods are presented in Sec. III and further described in detail in Appendixes A–C. Section IV consists of two separate parts, divided by the observed characteristic bubble shape during the collapse—mushroom-shaped bubble in Sec. IV A and pear-shaped bubble in Sec. IV B. The wrap-up of the findings is provided in a short conclusion.

## II. EXPERIMENTAL METHODS

### A. Experimental setup

The experimental setup for the generation and observation of laser-induced cavitation bubble dynamics near comparably sized cylindrical structures coupled the excitation fiber laser and the short-pulsed illumination system inside the experimental chamber. The excitation fiber laser emitted 60 ps laser pulses at a wavelength of 515 nm and was focused through a dipping microscope objective (Nikon CFI Apo NIR 40× W) with 40× magnification, a numerical aperture (*NA*) of 0.8, and a working distance of 3.5 mm. Plasma was formed, accompanied by a rapid localized rise in temperature, followed by shock wave emission and cavitation, converting some of the laser energy into mechanical energy, most relevantly for this work into bubble energy. An acousto-optic shutter was used to lower the laser pulse repetition frequency to 1 Hz, in order to allow for a steady state condition to be reached in water after each subsequent laser induced bubble (LIB) event and to give residual bubbles enough time to migrate out of the volume of interest. The water environment was controlled by using distilled water with a specified conductivity of 0.8 *μ*S/cm.

Geometrically, the dipping objective was inserted into the 3D-printed experimental chamber with a footprint of 20 × 20 × 20 mm^{3} from above [Fig. 1(A), $x\u2192$ direction], while the imaging cover-glass ports were inserted into the chamber walls in the $y\u2192$ direction. The cylindrical structure was created from a stripped optical fiber with 80 *μ*m diameter attached to transparent support structure, the latter positioned well out of the volume of interest: >10 γ. The dimensionless distance parameter γ is defined as the ratio of the bubble center distance from a solid surface *h* and bubble maximum radius in an infinite liquid $R0$: $\gamma =h/R0$. The original bubble radius in an infinite liquid is taken as a reference in accordance with recent literature,^{16,19} due to the changes in bubble morphology exhibited in the close vicinity of solid structures.

According to the coordinate system definition [Fig. 1(B)], the cylinder was viewed in either the $y\u2192$-axis direction [further referred to as the “frontal view,” Fig. 1(C)] or the $z\u2192$-axis direction [further referred to as the “side view,” Fig. 1(D)], in both cases resulting in a magnification of 1 *μ*m/pixel. A reference, free bubble was generated in an experimental approximation of an infinite liquid, with the cylindrical structure moved away from the bubble, to ensure uninterrupted dynamics observation. The laser pulse energy was controlled to generate a stable cavitation bubble with the maximal diameter $R0\u224893\u2009\xb1\u20092\u2009\mu m$. The whole bubble evolution was measured in time from the initial LIB to the second collapse and compared to two different models. The obtained data, including both model functions, is presented in Fig. 2. For the reference bubble, the parameter was set to $\gamma \u226b5$, while the effects of the nearby solid cylinder were measured at four different parameter values: $\gamma ={0.10,\u20090.32,\u20090.53,\u20090.74}$. The distance $h$ was measured optically, from the LIB center to the nearest point of the solid cylinder.

### B. Imaging

The imaging setup was based on a short-pulsed illumination system,^{30} focused into the area of interest by a 20 mm lens, with the resulting image captured by a combination of an externally triggered camera and a Mitutoyo M Plan Apo NIR 20× imaging objective with 20× magnification, *NA* of 0.40, and a working distance of 20.0 mm. The camera used was Photron Fastcam Nova S9, set to an externally triggered frame rate. The full camera resolution of 1024 × 1024 pixels at low frame rates was used, with the pixel size of 20 × 20 *μ*m.

The illumination system is based on a commercially available single-emitter laser diode module coupled into a delivery fiber and emits light at a nominal wavelength of 793 nm. The system has already been proven for monitoring the shockwave propagation near a concave surface^{31} and shockwave expanded nanobubbles.^{32} In the present experiments, the illumination pulse duration was set to 3 ns, and the illumination system had matched electronic clocks with the excitation laser, allowing for synchronization jitter below 1 ns. For each LIB event, two illumination pulses were emitted, and captured within a single, synchronized camera exposure. The first illumination pulse delay was varied with respect to the excitation pulse, because it was meant to capture the bubble evolution in time. The second illumination pulse was kept at a fixed delay near the moment of initial bubble collapse, ensuring that the observed LIB event was maximally similar to the rest of events [Figs. 1(B) and 1(C)]. This approach effectively superimposes the bubble shape at two separate times into a single frame. Using the technique, we have mostly eliminated the effects of excitation laser pulse energy fluctuations and energy conversion differences between the laser pulse and cavitation bubble. The resulting timing stability of the recorded moment of collapse was estimated to be ±20 ns. Events that did not correspond with the criteria given were not recorded. The bubble evolution in time was then stitched together from multiple separate LIB-initiated cavitation events, one event per each delay at each viewing direction and each $\gamma $. A complete bubble evolution sequence at a given γ, as shown in Fig. 2, is composed from bubble evolutions of 80 separate LIB events. The resulting effective frame rate varies with the speed of the bubble wall, from 1 MHz around the maximal bubble diameter, up to 20 MHz near the point of first collapse.

A complete initial bubble oscillation sequence imaged from just after the LIB event (timestamp 0.00 *μ*s, actual delay after the LIB is 5–10 ns, constituting an unknown, but constant delay), up to 18.20 *μ*s is shown in Fig. 2. The moment of collapse is at 18.00 *μ*s, showing the smallest experimentally obtainable feature during the initial bubble oscillation. Each subpanel shows the bubble imaged from two views at the same time delay.

## III. NUMERICAL METHODS

Simulations were used to provide insight beyond the resolution of the experiment. Bubble dynamics can be mathematically described through the laws of mass, momentum, and energy conservation. Presently, we consider the compressible two-phase flow with the following assumptions:

Two phases are considered—the gas bubble and the ambient liquid. Their interface is assumed to remain sharp through the bubble lifetime.

The gas phase is considered as air and ambient liquid as water. The presence of water vapor inside the bubble is neglected.

Mass transfer between both phases is neglected. The effects of condensation on bubble collapse and rebound intensity are indirectly considered through gradual reduction of equilibrium bubble radius during its lifetime.

Heat transfer between both phases is neglected.

Heat conduction within each phase and viscous dissipation are neglected.

Both phases are considered as Newtonian fluids.

Both phases are described through barotropic equations of state. Material compression and expansion is assumed as an isentropic process.

Material parameters, i.e., dynamic viscosity and surface tension, are considered as constant.

Surface tension is modeled as a body force acting at the phase interface. Other body forces are not considered.

Fiber is considered as a rigid and fixed cylindrical wall with the no-slip condition.

Conservation of mass and momentum is given in the following equations, whereas energy conservation is imposed through the previously given assumptions:

Here, $\rho $ and $p$ denote the density and pressure field, $U$ is the velocity vector field, $\tau $ is the viscous stress tensor, and $s$ is the body force term due to surface tension. $\u2207$ and $\u2207\u22c5$ correspond to the differential operators of the gradient and the divergence, whereas $\u2297$ denotes the tensor product. Viscous stress tensor $\tau $ can be for Newtonian fluids written as

where $\mu $ denotes the dynamic viscosity, and $I$ is the unit tensor. Surface tension term $s$ is considered according to the Continuum Surface Force model^{33} and can be expressed as

Here, $\sigma $ is surface tension, and $\kappa g$ the bubble surface curvature calculated as

Volume fraction field is denoted by $\alpha $, and a subscript is used to denote a phase specific quantity ($g$ denotes the gas phase, and $l$ denotes the liquid phase), e.g., $\rho g$ is the density of the gas phase. The considered material parameters of both phases are given in Appendix A.

The liquid phase is modeled as water through a modified version of the Tait's equation of state

Here, $\rho ref$ is the reference density at reference pressure $pref$. The terms $n$ and $B$ denote the density exponent and the Tait's pressure, respectively.

The gas phase is modeled as air through the Noble–Abel equation of state for an isentropic process in a calorically perfect gas, which states

Here, $k$ and $b$ denote the isentropic exponent and the co-volume of gas molecules. The constant value at the right hand-side of the equation is considered as $p\u221e(1\u2212k)Rg*T\u221ek$, which follows from the original form of the Noble–Abel equation of state. $Rg*$ denotes the specific gas constant, $T\u221e$ is the ambient temperature, while the co-volume is expressed as a fraction of the gas specific volume at equilibrium bubble conditions $b=\beta /\rho \u221e$.

For both fluids, the speed of sound is calculated as

Additionally, a correction term is included in the employed equation of state (see Appendix B), through which we can indirectly mimic the effects of condensation on increased bubble collapse intensity and decreased magnitude of bubble rebound.^{19}

The presently utilized fluid dynamics solver Ansys Fluent^{34} is based on the Finite Volume Method. A two-phase problem is numerically addressed through a single-fluid approach, namely, the Volume of Fluid methodology. There both phases are represented by a single fluid with effective properties and a shared pressure and velocity field. The bubble–liquid interface is captured by solving the continuity equation for the volume fraction field $\alpha l$ of the ambient liquid

The volume-averaged fluid properties $\varphi $ throughout the computational domain can be obtained as

Presently, this holds for density $\rho $ and dynamic viscosity $\mu $. Following that, a single momentum equation is solved, which yield the shared velocity field.

A pressure-based variant of the solver is employed, where the pressure-implicit with splitting of operators pressure-velocity coupling algorithm^{35} is used along with first order implicit temporal discretization. Both density and momentum are discretized by the second order upwind scheme, whereas the Pressure Staggering Option scheme^{36} is used for pressure discretization. Volume fraction field is discretized using an algebraic compressive scheme—a second order reconstruction scheme based on the slope limiter.^{34} Boundary conditions at the end of the computational domain are considered as a wave transmissive pressure outlet.

The geometrical configuration of the considered phenomenon imposes a 3D modeling approach, which greatly increases the computational burden in comparison to numerous other cases where dimensionality of the problem can be reduced, e.g., bubble dynamics near a flat solid boundary can be reduced by taking advantage of axial symmetry. The computational domain is represented by a rectangular block with a cut out cylinder of 80 *μ*m in diameter. Two symmetry planes are considered (x-y and z-x planes, also referred to as a side and frontal plane, see Fig. 1), which reduces the finally considered geometry to a quarter of the whole 3D domain. A structured computational mesh is used throughout the domain. Mesh spacing varies between 0.25 and 2 *μ*m in the bubble-fiber region, with mean cell spacing of approximately 1.5 *μ*m. Mesh resolution gradually increases to approximately 5 R_{max} at the outer boundaries of the computational domain, which are placed reasonably far away from the bubble-fiber region (approximately 50 R_{max}). The final computational cell count in is on the order of 5 × 10^{6}, which further justifies the consideration of two symmetry planes, as the full 3D model would consist of roughly 20 × 10^{6} cells. We do acknowledge that the considered mesh resolution is relatively coarse; however, the presently available computational resources pose a major constraint for consideration of finer grids. Variable time stepping is employed according to the Courant–Friedrichs–Lewy (CFL) condition of the maximum cell CFL number of 0.2.

Pressure at the outer domain boundaries is set to $p\u221e$ of one athmosphere. Symmetry boundary conditions are enforced at both planes of symmetry, whereas the no-slip condition is enforced at the cylinder surface. Bubble expansion is realized by an initial over-pressure of the bubble in relation to the ambient liquid. The internal bubble pressure is determined directly from the initial equilibrium radius R_{eq,0} and the initial bubble radius R_{0}, which can be understood as initial parameters that determine the magnitude of bubble expansion (see Appendix C). In addition, one needs to also determine the parameters that primarily govern the magnitude the bubble's rebound. In the present case, this is achieved by gradually reducing the equilibrium bubble radius after the bubble has reached its maximum size (see Appendixes A and B). A similar approach is common with the use of ordinary differential equation (ODE) based models but was also employed and validated with more advanced simulations in the past.^{16,19} Through this, one can closely match bubble collapse and rebound dynamics to the experiments in a relatively simple manner and without employing more complex modeling techniques that include consideration of heat and mass transfer mechanisms. It is true that the latter type of models might be generally preferred due to physical consistency and accuracy. However, in certain applications, simpler and less computationally demanding models might be preferable.

The employed computational fluid dynamics (CFD) model is first validated for the case of a spherical bubble, which also serves as means to determine the appropriate parameters ( Appendix A) and initial conditions ( Appendix C) for further calculations. The obtained results shown in Fig. 2 are compared to the experimental results (mean error of ±1.9 *μ*m) and the results of the Gilmore's model,^{37} which show a good level of agreement until the second bubble collapse, especially when one considers the estimated measurement uncertainty of ±2 *μ*m. The latter stems from a combination of imaging magnification (1 *μ*m/pixel) and deviations of the bubble shape from the spherical shape.

## IV. RESULTS AND DISCUSSION

A complete bubble evolution sequence at $\gamma =0.32$ is shown in Fig. 3. The early bubble evolution, up to around 200 ns post LIB, is nearly spherically symmetrical, and thus equivalent to bubble evolution in an infinite liquid. Upon closing the distance between the bubble wall and the solid cylinder surface down to a few micrometers, the bubble shape evolution is significantly defined by the solid cylinder shape. The bubble starts to engulf the cylinder and stretch along it. The maximal bubble volume is reached around 9.00 *μ*s after the LIB event.

Afterward, the bubble collapse starts, with the bubble wall speed increasing non-homogeneously, with the greatest movement detected around the cylinder structure, leading to an early formation of sharp bubble corners first visible between 11.00 and 12.00 *μ*s after the LIB event. This non-homogeneous evolution leads to a characteristic mushroom-shaped bubble development (16.00 *μ*s after the LIB event). In the final stages of the collapse, the bubble wall speed differences lead to a planar jet formation within approximately the last 200 ns of the initial bubble oscillation. The jet starts to form between the 17.70 and 17.80 *μ*s timestamps and is visible as an indentation in the top bubble wall on the next frames taken in the frontal view. The jet can hardly be detected in the side view due to its geometry, as clearly visualized in the simulations.

Bubble evolution sequences for other experimentally realized $\gamma $-values are shown in Figs. 4–6. We include shortened sequences, showing the bubble at 12 different timestamps, chosen to show the initial LIB positioning with respect to the solid cylinder, typical shape evolution, and the smallest volume at the time of collapse. The bubble evolution at the two lowest γ-values (i.e., 0.10 and 0.32) exhibit many similarities in shape, including the characteristic mushroom-shaped collapsing bubble (frontal view) and the elongated bubble shape along the $y\u2192$*-*axis visible in the side view perspective. The latter features steep side walls exhibited during the last 5 *μ*s of the initial oscillation (from around the 14.00 *μ*s timestamp onward). Similarities in bubble shape evolution were also observed between the bubble evolution at the two highest $\gamma $-values (i.e., 0.53 and 0.74). The characteristic mushroom-shaped collapsing bubble (frontal view) is either barely pronounced (at 0.53) or not present at all, the bubble reaching a pear-like shape (e.g., Fig. 5, timestamp 16.00 *μ*s), and the side view reveals a significant deviation from the previously elongated bubble shape in the $y\u2192$*-*axis, which now shows an elongation in the $x\u2192$*-*axis. In all cases, a jet toward the cylinder surface is formed in the last 200–300 ns of the initial bubble oscillation, most clearly indicated by the indentation in the top bubble wall visible in the frontal view and only marginally indicated in the same way in sideway imaging of experiments at the two lowest $\gamma $-values.

The experimental results point toward highly non-homogeneous flow formations during most of the bubble collapse phases, identified by sharp corner formation in the bubble shape, as well as by bubble wall speed differences around the bubble sides. The optical imaging of the bubble evolution does not allow for flow measurement. Therefore, we have employed numerical simulations, using the experimentally determined bubble shape evolution in time as a benchmark for evaluating the simulation performance and ability to match the situation within the experimental chamber with a reasonable accuracy.

### A. Mushroom-shaped bubble: Bubble in extreme vicinity of a cylinder—$\gamma \u2009=\u20090.10$ and $0.32$

When a bubble is generated in the extreme vicinity of the fiber ($\gamma =0.10$ and 0.32), it almost instantly begins to deviate from the initial spherical shape. The bubble expansion is initially much more pronounced in the direction away from the fiber (Fig. 7 at 1.00 *μ*s), as one might intuitively expect. Images of the frontal plane view hint that the upper half of the expanding bubble retains a hemispherical shape during expansion (until 10.00 *μ*s), which would be consistent with the previous results of the bubble–sphere interaction studies.^{21,22} This, however, is not exactly the case, as the side plane view reveals bubble flattening and elongation along the fiber axis. Once the bubble reaches its maximal volume ($t\u22489\u2009\mu s$), a significant portion of the fiber is engulfed by the bubble, which reaches below the center point of the cylinder (Fig. 7 at 4.00 and 10.00 *μ*s). This feature later leads to the formation of a previously documented mushroom-shaped bubble.^{19–22,24} As the bubble begins to contract, the flow along the curved wall of the cylinder accelerates in the vertical direction [Fig. 7 at 10.00 *μ*s and Fig. 8(a)]. Two distinct shape features begin to form in the frontal plane—bubble stem and cap, resembling a mushroom-shape [Fig. 7, frontal view at 14.00, 16.50 *μ*s and Fig. 8(b)]. In Figs. 8 and 9, one can notice that the obtained numerical solution implies that when the bubble is in the extreme vicinity of the fiber ($\gamma =0.10$), the liquid film coating the cylindrical surface breaks and the gaseous bubble contents come into direct contact with the wall, i.e., the bubble sticks to the wall. This is also the reason for the occurrence of circular features inside the region where the bubble is located. However, we believe that this is an artifact of the employed numerical resolution, which would be overcome by using a finer mesh in the region where the bubble is present. Although a finer mesh spacing was used closer to the wall (down to 0.25 *μ*m), this problem persists for $\gamma =0.10$. Thus, a finer mesh would have to be employed across the entire bubble region to retain high gradients of volume fraction fields at the bubble–liquid interface.

During bubble contraction, the upward facing vertical flow gradually changes the direction toward the vertical axis, which results in a convergent radial flow (Fig. 7, frontal view at 16.50 *μ*s), and the bubble further shrinks in a way that resembles closing of a folding hand fan. This can be further seen at 17.00 and 17.50 *μ*s [Figs. 8(c) and 8(d)], where the bubble cap develops an indentation, i.e., a “kink,” which is typical for a collapsing bubble in the close vicinity ($\gamma <0.7$) of a flat surface.^{16,18} However, the kink is not axisymmetric, as in the side plane, the horizontal contraction is predominant in the mid part of the bubble. The radial bubble contraction eventually leads to the impact of the opposite bubble walls at the vertical axis [Fig. 7, frontal plane at 17.80 *μ*s, and Figs. 9(b) and 9(c)]. As a result, the convergent radial flow is focused into a jet toward the cylinder [Fig. 7, frontal plane at 17.85 *μ*s and later, Figs. 9(c)–9(e)].

On the other hand, different bubble shape features can be seen in the side plane view, orthogonal to the frontal plane. One can notice that the bubble wall above the cylinder becomes practically vertical during the collapse (Fig. 7, side view at 16.50 *μ*s) and that the bubble contraction is slower than in the frontal plane [Figs. 8(b) and 8(c)]. Also here, a so-called kink is developed [Fig. 7, side view at 17.80 *μ*s and Figs. 9(a)–9(c)], although it is less apparent and does not lead to the impact of the opposite bubble walls, as the jet toward the wall has already been formed. Judging from the experimental sequences, one cannot confirm the presence of a jet from the side view. Nevertheless, the results of the numerical simulation do reveal it just before the bubble collapse at 18.05 *μ*s (Fig. 7, side view).

Looking at the velocity contours and 3D bubble shape from the numerical simulation, one can notice that the jet develops at 17.80 *μ*s, after bubble wall collision along the side plane [Fig. 9(c)]. This mechanism is very similar to the onset of a recently documented fast, thin jet,^{15–18} which develops when a bubble collapses in the extreme vicinity of a flat wall and can reach supersonic velocities. An important difference, however, is the nature of radial flow focusing, which is dependent on the symmetry conditions. In the case of a fast thin jet, the bubble assumes an axisymmetric shape, and the flow is focused toward and along the vertical axis of symmetry. In the present case, the bubble shape is not axially symmetric, and the flow is focused toward and then along the side plane (x-y symmetry plane), because the bubble contraction is more pronounced in the orthogonal direction (along the frontal plane). This leads to the formation of a planar jet with two vastly different horizontal dimensions—width of 5 *μ*m and length of 45 *μ*m [Fig. 9(d)]. At this point, the numerical simulations predict the peak jet velocity on the order of 200 m/s. As the bubble continues to collapse, the jet further accelerates toward the wall reaching 250 m/s and beyond 300 m/s after the impact into the fiber. Along with jet acceleration, one can also observe its gradual change in shape from the initial elongated shape (thin in the frontal plane and wide in the side plane) to a more even cross section [diameter of 20 *μ*m, Fig. 9(f)]. This implies that the initial elongated shape is unstable, which can be attributed to two factors. First, this can be attributed to the effects of the remaining radial bubble contraction along the side plane, which meets the downward flowing jet and causes the rearrangement of its cross section. The shape is expanded along the frontal plane and contracted orthogonally to it (along the side plane), since water is only weakly compressible. Interestingly, the saddle-like bubble shape also reverses its orientation, which can be clearly seen upon comparison of bubble shapes in Figs. 9(e) and 9(g). Second, as the jet propagates along the cylinder, this eventually leads to the impact and outflow along the cylinder walls. It comes as no surprise that the outflow is more pronounced along the curved walls of the cylinder [frontal plane, Fig. 9(h)], as the streamlines must undergo a lesser deviation from the initial downward direction and the flow is thus less restrained by the curved cylinder wall.

The peak jet velocities are on the order of 300 m/s, which upon its impact into the fiber wall results in highly transient (∼10 ns) water hammer pressure with the peak of ∼500 MPa. This value was evaluated as a maximum pressure at the wall during the water jet impact. Although the obtained numerical bubble shapes correspond very well to the experimental bubbles, we cannot directly validate the obtained velocity values due to the obvious constraints of experiments at such small spatiotemporal scales. Additionally, the employed resolution of the numerical simulations is still relatively coarse, and thus, the reported values of peak jet velocities and impact pressures should only be interpreted as the lower bounds of the actual values. Further and more detailed experimental and numerical studies are, thus, needed to provide more insight into the validity of the obtained values and the details of bubble shape progression in the final stages of the collapse.

### B. Pear-shaped bubble: Bubble in moderate vicinity of a cylinder—$\gamma =0.74$

Bubbles generated in the moderate vicinity of the fiber ($\gamma =0.74$) initially exhibit more uniform expansion and begin to significantly deviate from the spherical shape only after the bottom bubble wall reaches the fiber and begins to engulf it at 4.00 *μ*s (Fig. 10, frontal view). Furthermore, judging from the side-plane view, the shape remains practically spherical even at the time of maximum bubble volume at ∼9 *μ*s (Fig. 10, side view at 10.00 *μ*s). However, this does not remain true during the bubble collapse phase (after 10 *μ*s). While most of the bubble wall is already in the contraction phase, one can notice that a large vortex in the side plane still drives the widening of the bubble base just above the cylinder wall [Fig. 11(a)]. At this instance, a smaller vortex of opposite orientation can be seen in the frontal plane along with accelerated inflow along the curved cylinder wall, which eventually leads to the formation of an upside-down pear-like bubble shape (Fig. 10, frontal view at 17.00 *μ*s). Although subtle changes in local bubble wall curvature can be observed [Fig. 10, side and frontal plane at 16.00–18.00 *μ*s, and Figs. 11(b) and 11(c)], they do not fully develop into a more pronounced kink as with smaller values of gamma.

During further bubble contraction, its bottom part remains attached to the cylinder, which leads to the uneven velocity distribution at the bubble wall [Fig. 11(b)]. The top part of the bubble further accelerates toward the cylinder [Fig. 11(c)] and the downward flow is eventually focused into a vertical jet toward the cylinder [Fig. 11(d)]. This closely resembles a well-known phenomenon of near-wall bubble collapse, where an axial jet toward the boundary is developed.^{38} Here, the obtained results are also consistent with the bubble–sphere interaction studies, where bubbles in the moderate vicinity developed a pear-like shape along with uniaxial jets toward the sphere.^{21,22} However, some differences can be observed in the present case, which is in part also due to different symmetry conditions. Here, the flow in the side plane lags behind, and thus, the developed jet initially assumes a planar shape, with the main dimension being along the axis of the cylinder [Figs. 12(a)–12(c)]. The jet velocity reaches 300 m/s; however, its velocity does steadily decline with propagation toward the wall and the peak resulting water hammer pressure is approximately 100 MPa. Overall, the formed jet characteristics are similar to the previously described case with $\gamma =0.10$. Nevertheless, the jet shape in the present case is more stable and is not significantly affected by the remaining horizontal bubble contraction along the side plane—preserving its planar shape until the impact [Figs. 12(a)–12(c)]. After the impact, the jet radially outflows past the cylinder, with a more pronounced flow along the curved part of the cylinder [Fig. 12(d)], as already observed at $\gamma =0.10$.

## V. CONCLUSIONS

In this work, we presented the interaction of microbubbles with nearby solid fibers of comparable sizes, with bubbles reaching up to 100 *μ*m radius and the fiber being 80 *μ*m in diameter. Both experimental and numerical work was performed, focused on small distances between the bubble and the fiber, where the physical interaction between the subjects is the strongest. For all four different γ-values considered within the scope of this work, we found very good agreement of numerical simulations with experiments.

Two modes of bubble collapse were identified, resulting in mushroom-shaped bubbles at positions closest to the fiber ($\gamma =0.10$ and 0.32) and a pear-shaped bubble at the farthest position ($\gamma =0.74$). At $\gamma =0.53$, features of both modes have been observed and can thus be defined as the transition between modes. Similarities were found with cases of bubble collapsing near-wall and near-sphere. However, some differences do exist. This can be attributed primarily to different symmetry conditions, as the present case does not entail axial symmetry. Notably, all four cases resulted in a planar jet formation toward the fiber. The formed jet initially assumes an elongated shape, whereas its stability depends on the mode of bubble collapse.

Further analysis of the planar jet, as the defining feature of the collapse, showed the peak jet velocities between 250 and 330 m/s and the resulting impact pressures between 100 and 500 MPa. The obtained values are in the lower range of the actual values, as the employed numerical resolution is relatively low due to high demand on computational resources for a fully 3D simulation. Further numerical and experimental studies are needed to gain more insight into the details surrounding the final stages of bubble collapse. Another open question is the scalability of jet formation to smaller and/or larger dimensions.

## ACKNOWLEDGMENTS

The authors acknowledge the financial support from the European Research Council (ERC) under the European Union's Framework Program for Research and Innovation, Horizon 2020 (Grant Agreement No. 771567-CABUM) and the Slovenian Research Agency ARRS (Core Funding Nos. P2-0422, P2-0270, and Project No. J2-3057).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Jaka Mur:** Data curation (equal); Investigation (equal); Writing – original draft (equal). **Vid Agrež:** Conceptualization (equal); Methodology (equal); Writing – original draft (equal). **Jure Zevnik:** Data curation (equal); Investigation (equal); Writing – original draft (equal). **Rok Petkovsek:** Funding acquisition (equal); Methodology (equal); Project administration (equal); Resources (equal). **Matevz Dular:** Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Writing – original draft (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: CONSIDERED MATERIAL PARAMETERS

The considered material parameters of both phases (Table I).

Phase . | Quantity (Unit) . | Value . |
---|---|---|

Liquid–water | Reference pressure $pref$ (Pa) | 101 325 |

Reference density $\rho ref$ (kg/m^{3}) | 998.2 | |

Tait's pressure $B$ (Pa) | 3.0 × 10^{8} | |

Density exponent $n$ | 7.15 | |

Dynamic viscosity $\mu l$ (Pa s) | 0.001 | |

Gas–air | Specific gas constant $Rg*$ (J/kg K) | 287.06 |

Ambient temperature $T\u221e$ (K) | 293.15 | |

Ambient pressure $p\u221e$ (Pa) | 101 325 | |

Adiabatic exponent $k$ | 1.4 | |

Dynamic viscosity $\mu g$ (Pa s) | 1.8 × 10^{−5} | |

Normalized co-volume $\beta $ | 0.0015 | |

Surface tension $\sigma $ (N/m) | 0.0728 | |

Initial equilibrium radius $Req0$ (m) | 4.25 × 10^{−5} | |

Final equilibrium radius $Reqfin$ (m) | 1.46 × 10^{−5} | |

Time parameter $ta$ (s) | 1.0 × 10^{−5} | |

Time parameter $tb$ (s) | 1.2 × 10^{−5} |

Phase . | Quantity (Unit) . | Value . |
---|---|---|

Liquid–water | Reference pressure $pref$ (Pa) | 101 325 |

Reference density $\rho ref$ (kg/m^{3}) | 998.2 | |

Tait's pressure $B$ (Pa) | 3.0 × 10^{8} | |

Density exponent $n$ | 7.15 | |

Dynamic viscosity $\mu l$ (Pa s) | 0.001 | |

Gas–air | Specific gas constant $Rg*$ (J/kg K) | 287.06 |

Ambient temperature $T\u221e$ (K) | 293.15 | |

Ambient pressure $p\u221e$ (Pa) | 101 325 | |

Adiabatic exponent $k$ | 1.4 | |

Dynamic viscosity $\mu g$ (Pa s) | 1.8 × 10^{−5} | |

Normalized co-volume $\beta $ | 0.0015 | |

Surface tension $\sigma $ (N/m) | 0.0728 | |

Initial equilibrium radius $Req0$ (m) | 4.25 × 10^{−5} | |

Final equilibrium radius $Reqfin$ (m) | 1.46 × 10^{−5} | |

Time parameter $ta$ (s) | 1.0 × 10^{−5} | |

Time parameter $tb$ (s) | 1.2 × 10^{−5} |

### APPENDIX B: MODIFICATION OF EQUATION OF STATE FOR THE GAS PHASE

Laser-induced bubbles generally consist of both water vapor and incondensable gasses. However, the exact amount of each constituent is generally not known, and their adequate consideration in numerical models remains one of the challenges up to this day. By considering bubble contents as air, we neglect the bubble's vapor content and the underlying mass transfer mechanisms, namely, evaporation and condensation. For this reason, an additional term $mod$ is included in the employed equation of state for the gas phase, through which we can indirectly mimic the effects of condensation on increased bubble collapse intensity and decreased magnitude of bubble rebound.^{16,19} The term $mod$ is represented by a ratio of initial bubble mass $mg0$ and instantaneous bubble mass $mg$ and can be written as

The latter is achieved through a gradual reduction of equilibrium bubble radius $Req$ during its lifetime in the following manner:

Here, $Req0$ and $Reqfin$ denote the initial and final equilibrium bubble radius. A smooth transition between both values is imposed according to a smoothstep function between times $ta$ and $tb$. Density and pressure at equilibrium bubble conditions can be written as

The finally employed equation of state for the gas phase can be written in the following form:

### APPENDIX C: DETERMINATION OF INITIAL CONDITIONS

The initial equilibrium radius $Req0$ along with the initial bubble radius $R0$ can be understood as initial conditions that determine the magnitude of bubble expansion. The latter is realized due to the bubble being over-pressurized in relation to the ambient liquid. The internal bubble pressure is determined directly from parameters $Req0$ and $R0$ as

where an initially spherical bubble in ambient liquid at rest is considered.

In the present case, the initial bubble radius $R0$ was chosen at 10 *μ*m based on the available mesh resolution and an additional constraint of being able to simulate bubble dynamics for $\gamma \u22650.1$. It is worth mentioning that smaller values of γ would be possible to consider, although a new set of initial conditions ($R0$, $Req0$) would have to be determined along with a reconsideration of the employed mesh resolution.

The initial equilibrium radius in then determined through a trial-and-error approach for the case of a spherical bubble. This can be performed until the desired level of error is obtained with respect to the experimental data. Presently, only a few simulation cycles were required to achieve a good match with the experiments, as a good initial guess can be obtained by solving simpler models, e.g., Gilmore's model.