Ring Vortex Dynamics Following Jet Formation of a Bubble Expanding and Collapsing Close to a Flat Solid Boundary Visualized via Dye Advection in the Framework of OpenFOAM

Abstract

A bubble expanding and collapsing near a solid boundary develops a liquid jet toward the boundary. The jet leaves a torus bubble and induces vortices in the liquid that persist long after the bubble oscillations have ceased. The vortices are studied numerically in axial symmetry and compared to experiments in the literature. The flow field is visualized with different methods: vorticity with superimposed flow-direction arrows for maps at a time instant and colored-liquid-layer flow-field maps (dye advection) for following the complete long-term fluid flow up to a chosen time since bubble generation. Bubbles with equal energy—maximum radius in a free liquid $R_{\max }^{\infty }$= 500 µm—are studied for different distances $D_{\mathrm{init}}$ from the solid boundary. The interval of normalized distances $D^{*}$ = $D_{\mathrm{init}}$/$R_{\max }^{\infty }$ from 0.4 to 1.8 is covered. Two types of vortices were reported in experiments, one moving toward the solid boundary and one moving away from it. This finding is reproduced numerically with higher resolution of the flow field and in more detail. The higher detail reveals that the two types of vortices have different rotation directions and coexist with individually varying vorticity amplitude throughout the interval studied. In a quite narrow part of the interval, the two types change their strength and extent with the result of a reversal of the dominating rotational direction of the fluid flow. Thereby, the experimentally found transition interval could be reproduced and refined. It is interesting to note that in the vortex transition interval, the erosion of a solid surface is strongly augmented.

1. Introduction

The fluid flow around bubbles in liquids develops complicated structures as soon as the bubbles do not stay spherical. A prominent example is a bubble collapsing in the neighborhood of a solid boundary. The most spectacular phenomenon is the liquid jet that develops toward the solid from the distal side of the bubble by involution of the bubble surface (for mainly numerical studies, see, e.g., [1,2,3,4,5,6]). The jet traverses the bubble interior, hits the opposite bubble wall and pushes it ahead. That way, the bubble attains the form of a gaseous torus with a liquid flow through the torus hole (the jet flow). This “standard jet” develops velocities of the order of 100 m/s.

Recently, a different type of jet formation toward a solid surface has been found for bubbles expanding and collapsing very close to the solid: the “thin, fast jet”. It was found first numerically (Lechner et al. [5,7]) and soon after experimentally (Koch [8] and Koch et al. [6]) with later confirmation (Reuter et al. [9]). The fast jet is not produced by involution of the distal side of the bubble but by a ring-shaped involution with self-impact at the axis of symmetry. The fast jet develops velocities of the order of 1000 m/s.

In addition to the jet flow with torus bubble formation, further distinct liquid flows, notably vortices, are generated in the liquid in the vicinity of the bubble. Lauterborn et al. [10] posed the question of whether a collapsing bubble is a short-living vortex. To really decide the question, advanced methods had to be developed in fluid flow visualization. Lauterborn et al. [11] discuss how fluid flows could be measured and visualized. As the methods were all not well suited to the fluid flows that appear in bubble dynamics, Vogel et al. [12,13] developed a time-resolved particle-image-velocimetry (PIV) method to study the flow in the vicinity of a collapsing cavitation bubble. The method consists of a sequence of double-exposure and multi-exposure photographs of tracer particles to give a time sequence of flow vectors in space. It has since then been developed further significantly for all kinds of flows, resolving both space and time (for an overview see, e.g., [14,15,16,17] and references therein, as well as, e.g., [18,19,20,21,22] for a very small selection of recent developments).

Previously, a distinction was made for applications, where the density of the tracer particles of the PIV method is too high, by introducing the particle-tracking-velocimetry (PTV) method with lower particle densities. In case of the PTV, instead of making use of the cross-correlation of the subsequent photographs, the particles are tracked individually. Nowadays, the two methods share an overlap of applicability (see, e.g., [20]) and indicate whether the image-processing algorithm makes use of the cross-correlation between subsequent images (PIV) or performs Lagrangian particle tracking (PTV).

The PTV method has also been applied to acoustically driven bubbles. Bolaños-Jiménez et al. [23] applied the astigmatic particle tracking velocimetry (APTV) [24] to acoustically driven bubbles in small holes in a plane solid to determine the induced flow (acoustic streaming) of a radially oscillating bubble. In this method, spherical tracer particles are imaged as ellipsoids. The elongation of the ellipsoid depends on the depth of the optical path by using the astigmatism induced by a cylindrical lens. This way, depth information is added to normal particle tracking velocimetry. Three-dimensional (3D) information of the flow can then be obtained from a succession of high-speed images. As a result, it was found experimentally that the flow of a fixed, acoustically driven bubble is rotationally symmetric about the axis through the bubble center perpendicular to the surface. Fauconnier et al. [25] applied the PTV method to a single, wall-attached, acoustically driven bubble undergoing asymmetric shape modes, which are excited subharmonically. The authors were mainly interested in the particle trajectories, more than in the velocity field. Hence, for the main body of their study, the method was simplified to backlight illumination, where the particles appear as small dark spots. After recording an image series, the minimum value of the brightness for each pixel over time is plotted. That way, the authors revealed the flower-like micro-streaming patterns around the bubble that oscillates in spherical harmonics.

The focus of the present study is on single, laser-generated bubbles in the vicinity of a plane, solid boundary. The PTV method has been applied for laser-induced, single cavitation bubbles in the literature. A low particle density is favorable in this case, because a higher particle density would obstruct the laser pulse producing the bubble-forming plasma via optical breakdown. Kröninger and Kröninger et al. [26,27] made use of fluorescent tracer particles and double-exposure photography to give more details about the fluid flow forming the jet in the case of a bubble collapse near solid boundaries. In particular, they experimentally determined flow velocities in the jet very near to the bubble and very near the time of complete collapse. Figure 1 gives an example. The photograph of the bubble and the velocity field was taken 10 µs before the collapse. The velocity arrows were gained by PTV with fluorescent tracer particles. The velocity is seen to reach about 30 m/s. This is consistent with numerical calculations in view of the time of 10 µs before collapse. The corresponding numerically found maximum jet speed is about 60 m/s [5,28]. The normalized distance $D^{*}$ in the experiment of Figure 1 is 0.5. $D^{*}$ is defined as [29]

$$\begin{array}{c}\hfill D^{*}=D_{\mathrm{init}}/R_{\max }^{\infty },\end{array}$$

(1)

where $D_{\mathrm{init}}$ is the initial distance of the bubble center (focus of the laser light in the experiment) to the solid boundary and $R_{\max }^{\infty }$ is the maximum radius the bubble would attain by expansion in an unbounded liquid. In the present work, $R_{\max }^{\infty }=500\mu \mathrm{m}$ is used throughout.

The results from a PTV/PIV measurement, i.e., the obtained velocity field, can be used to visualize the temporal evolution of the flow in a post-processing step. Reuter et al. [30] made use of the PTV method to follow the flow field around such a bubble, like the one in Figure 1, also long after the collapse. During post-processing, the authors turned the obtained velocity field into a color-coded map—just as in dye advection, which is a form of the Lagrangian–Eulerian advection (LEA), as described by Jobard et al. [31] and Laramee et al. [32] —mimicking the disturbance of a layer of submerged ink by the bubble. Although the authors computed the texture warping with a forward Eulerian time-step and millions of particle traces instead of using the optimized algorithm described for LEA, the dynamical evolution of the flow is visualized, because each fluid parcel is given a specific color that is advected with it. They found that the jet generates a liquid vortex that, depending on the distance of the bubble to the solid boundary, travels either orthogonally away from the solid boundary or toward the solid boundary. These two types of vortices were called by them a free vortex and wall vortex, respectively.

Saini et al. [33] showed one single experimental case of a collapsing bubble cap (negative value for $D^{*}$ according to Equation (1)) disturbing a dye layer. They found a vortex traveling orthogonally away from the solid boundary. The authors calculated the corresponding vorticity field by means of their in-house solver. As will be shown in the present work, the vortices generated by the bubble in the work of Saini et al. [33] stem from different dynamics than the one induced by the bubble in the work of Reuter et al. [30].

Similar vortex dynamics as the ones in Reuter et al. [30] were observed by Sieber et al. [34] with a bubble close to a sand bed. The sand gravels acted as tracer particles after the bubble collapse and were partially carried away from the bed by vortex shedding after the bubble dynamics have ceased.

Figure 1. Example of an experimental photograph of a jetting bubble with its jet flow close to a solid boundary shortly (10 µs) before collapse (reprinted from [26] under license CC BY-NC 2.0 [35]). The velocity arrows are gained by the particle-tracking velocimetry method. The white cone above the bubble (black) is the plasma of the optical breakdown at $t=0$. $D^{*}$= 0.5.

Figure 1. Example of an experimental photograph of a jetting bubble with its jet flow close to a solid boundary shortly (10 µs) before collapse (reprinted from [26] under license CC BY-NC 2.0 [35]). The velocity arrows are gained by the particle-tracking velocimetry method. The white cone above the bubble (black) is the plasma of the optical breakdown at $t=0$. $D^{*}$= 0.5.

In the present numerical study, the experimental findings of Reuter et al. [30] on vortex formation and motion are compared to numerical work. The simulations are performed using a modified finite-volume CFD-solver from the framework of OpenFOAM [36,37]. In this framework, the concept of a passive scalar exists, meaning a scalar quantity that is advected passively with the flow, not interfering with it. This is completed via solving the advection equation for the passive scalar during solver run-time. The method has been termed color-layer tracer-field advection when applied to bubble dynamics [28]. It is further described in the body of the present work. Here, the simplified and historically more correct term dye advection is used because the flow-field visualization concept of the color layer tracer-field advection method falls into the subcategories dye injection and dye advection of LEA [31,32]. The difference is that no streamlines or cross-correlations need to be computed anymore since the visualization is computed by solving a partial differential equation during the solver run-time.

The article is organized as follows. In Section 2, the bubble model and the numerical solution method are introduced with more information and definitions in Section 2.1 and Section 2.2. In Section 2.3, the reason for the choice of time instants for vortex visualization is given, as the presentation of the full sequence is too space filling. The two types of vortex visualization used here (flow and vorticity maps) are presented in Section 2.4. In Section 3, the results are presented with series of flow fields with their vortices and a comparison with experiments from the literature (Reuter et al. [30]). The results are discussed in Section 4 in view of the erosion problem, among others. Section 5 outlines the conclusions to be drawn, together with an outlook to further research.

2. Bubble Model and Numerical Solution Method

A bubble model for a cold liquid (a liquid far from its boiling point, cf. [38]) with the following properties is used. The bubble contains a small amount of non-condensable gas to comply with experiments [39,40,41]. The vapor pressure is small compared to the ambient pressure and is neglected. The liquid is taken as compressible for the inclusion of pressure waves up to weak shock waves, as a significant amount of energy is radiated in the form of acoustic waves (see, e.g., [12,37,42,43,44,45]). Thermodynamic effects and mass exchange through the bubble wall are neglected. Gravity can be omitted due to the small size of the bubble.

The equations of motion of the two-phase flow are formulated in the “one-fluid” approach, i.e., with one density field $\rho (\mathit{x},t)$, one velocity field $\mathit{U}(\mathit{x},t)$, and one pressure field $p(\mathit{x},t)$, satisfying the Navier–Stokes Equation (2), including a term for the surface tension $\sigma$ and the continuity Equation (3) [46,47,48]:

$$\begin{array}{ccc}\hfill \frac{\partial \left(\rho \mathit{U}\right)}{\partial t}+\nabla ·(\rho \mathit{U}\otimes \mathit{U})& =& -\nabla p+\nabla ·\left[\mu \left(\nabla \mathit{U}+{(\nabla \mathit{U})}^T-\frac{2}{3}(\nabla ·\mathit{U})\mathbb{I}\right)\right]\hfill \\ & & +\sigma {\kappa }_v\nabla {\alpha }_l\mathrm{where}{\kappa }_v=-\nabla ·\left(\nabla {\alpha }_l/\left|\nabla {\alpha }_l\right|\right),\hfill \\ \hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathit{U}\right)& =& 0.\hfill \end{array}$$

(3)

∇ denotes the gradient, $\nabla ·$ is the divergence, and ⊗ is the tensorial product. $\mu (\mathit{x},t)$ denotes the viscosity field and $\mathbb{I}$ denotes the unit tensor. The surface tension coefficient $\sigma$ of water is set to $\sigma =0.0725\mathrm{kg}{\mathrm{m}}^{-2}$. With the one-fluid approach, it is possible to use the Volume of Fluid (VoF) method, where volume fraction fields ${\alpha }_l(\mathit{x},t)$ and ${\alpha }_g(\mathit{x},t)$ are introduced to distinguish between the two fluids, liquid (l) and gas (g), with ${\alpha }_l=1$ in the liquid phase, ${\alpha }_l=0$ in the gas phase, and ${\alpha }_g(\mathit{x},t)=1-{\alpha }_l(\mathit{x},t)$. The position of the interface is then given implicitly by the transition of ${\alpha }_l$ from 1 to 0. The viscosity field $\mu (\mathit{x},t)$ can be written as $\mu (\mathit{x},t)={\alpha }_l(\mathit{x},t){\mu }_l+{\alpha }_g(\mathit{x},t){\mu }_g$ (see, e.g., [49]). The dynamic viscosities ${\mu }_l$ of the liquid and ${\mu }_g$ of the gas are taken to be constant (${\mu }_l=1.002\times {10}^{-3}\mathrm{kg}{\mathrm{s}}^{-1}{\mathrm{m}}^{-1}$, ${\mu }_g=1.86\times {10}^{-5}\mathrm{kg}{\mathrm{s}}^{-1}{\mathrm{m}}^{-1}$). The density field $\rho (\mathit{x},t)$ is given by $\rho (\mathit{x},t)={\alpha }_l(\mathit{x},t){\rho }_l(\mathit{x},t)+{\alpha }_g(\mathit{x},t){\rho }_g(\mathit{x},t)$ with ${\rho }_l$ and ${\rho }_g$ denoting the densities of the liquid and gas, respectively. As there is no mass transfer between the bubble interior (gas) and exterior (liquid), the respective phase-fraction density fields ${\alpha }_l{\rho }_l$ and ${\alpha }_g{\rho }_g$ separately obey the continuity equation:

$$\frac{\partial \left({\alpha }_i{\rho }_i\right)}{\partial t}+\nabla ·\left({\alpha }_i{\rho }_i\mathit{U}\right)=0,i=l,g.$$

(4)

The equations of motion are closed by the equations of state for the gas and the liquid. For the gas in the bubble, the change of state is assumed to be adiabatic,

$$\begin{array}{c}\hfill {\rho }_g\left(p\right)={\rho }_{gn}{\left({\displaystyle p/p_n}\right)}^{1/{\gamma }_g},\end{array}$$

(5)

with $p_n$ and ${\rho }_{gn}$ representing the pressure and the density of the gas in the bubble at normal conditions, respectively, and ${\gamma }_g=1.4$ indicating the ratio of the specific heats of the gas (air). For the liquid, the Tait equation of state for water is used (see, e.g., [39]):

$$\begin{array}{c}\hfill {\rho }_l\left(p\right)={\rho }_{\infty }{\left((p+B)/(p_{\infty }+B)\right)}^{1/n_T},\end{array}$$

(6)

with $p_{\infty }$ representing the atmospheric pressure, ${\rho }_{\infty }$ representing the equilibrium density, the Tait exponent $n_{\mathrm{T}}=7.15$ and the Tait pressure $B=305$ MPa.

The pressure-based two-phase solver compressibleInterFoam of the open source package OpenFOAM, precisely foam-extend, is used for the implementation of the equations. For details, see the description and validation in, e.g., Koch et al. [28,37] and Lechner et al. [5]. Discretization of the above partial differential equations is performed with the finite volume method [48,50].

2.1. Mesh, Initial Conditions, Boundary Conditions and Time Stepping

Simulations are carried out in axial symmetry on a computational domain shaped like a half circle with a radius of $100R_{\max }^{\infty }$. The domain is meshed such that the cells spread radially outward and align to the center core with an edge length of 80 µm, where the cells are oriented in a Cartesian manner. A sketch of the mesh is shown in Appendix A. The idea is to align the cells as much as possible to the bubble interface, i.e., a polar orientation with an apex at the initial bubble center, while at the same time dissipating the shock wave in the outer regions and avoiding a high total cell amount for the whole mesh.

The initial data for the bubble are set in the following way. The compressed bubble is initiated at a distance $D_{\mathrm{init}}$ from the solid boundary. The normalized distance $D^{*}$ is defined as in Equation (1). From the user-given initial bubble radius $R_{\mathrm{init}}\begin{array}{c}\hfill =20\mu \mathrm{m}\end{array}$, the bubble volume is calculated analytically. The cells within that radius are set to ${\alpha }_l=0$ (gas). The discretized volume of the gas is determined numerically. The pressure then is recomputed to yield the same compression energy as given by the original set of parameters $R_{\mathrm{init}}$ and $R_n$, $R_n$ being the radius of the bubble at equilibrium. The resulting adaption of the pressure is usually in the range of 1 %, but it matters in terms of convergence [8]. The bubble pressure at $t=0$ in the cases investigated is $p_{\mathrm{init}}=1.132\mathrm{GPa}$ for a sphere with $R_{\mathrm{init}}=20\mu \mathrm{m}$, $R_n=184.1\mu \mathrm{m}$ and the smallest cell edge length of 2 µm in the mesh center core. The technical details about bubble mass correction methods and bubble mass controls are given in Appendix C and Appendix D. The parameters for the fluids entering the calculation are given in

Table A2. Fluid properties entering the calculation. ${\mathrm{R}}_{\mathrm{spec}}$ is the specific gas constant of air.

Fluid	$\mathit{\mu }$ [$\mathbf{kg}{\mathbf{s}}^{-1}{\mathbf{m}}^{-1}$]	$\mathit{\rho }(\mathit{p}={\mathit{p}}_{\mathit{\infty }})$ [$\mathbf{kg}{\mathbf{m}}^{-3}$]	${\mathit{n}}_{\mathit{T}}$; ${\mathit{\gamma }}_{\mathit{g}}$	B [Pa]; ${\mathbf{R}}_{\mathbf{spec}}$ [$\mathbf{J}{\mathbf{kg}}^{-1}{\mathbf{K}}^{-1}$]
liquid	$1.002\times {10}^{-3}$	998.2	7.15	$305\times {10}^6$
gas	$18.6\times {10}^{-6}$	derived	1.4	287.0
	$\mathbf{\sigma }$[$\mathrm{kg}{\mathrm{m}}^{-\mathbf{2}}$]	${\mathit{p}}_{\infty }$[Pa]	${\mathbf{T}}_{\mathbf{0}}$[K]
	0.0725	101315	293.15

.

The discretized equations are solved implicitly in a segregated manner. Therefore, the solver allows for large time steps under certain circumstances. In the present work, the flow patterns are of interest that persist long after the bubble oscillations have ceased. The simulation therefore has to pass a long period of time where the dynamics are not violent anymore. Hence, a suitable time stepping is important for the different phases of the bubble dynamics for a trade-off between numerical stability and computational effort. The Courant number Co is calculated to determine the time step size $\mathsf{\Delta }t$ for adaptive time stepping. It is defined as the spatial maximum of the ratio of the local flow speed $\mathbf{U}\left(\mathbf{x}\right)$ at location $\mathbf{x}$ to the maximum resolvable flow speed by the spatial and temporal discretization ($\mathsf{\Delta }x\left(\mathbf{x}\right)/\mathsf{\Delta }t$) [51]:

$$\begin{array}{c}\hfill \mathrm{Co}=\underset{\mathbf{x}}{\max }\left(\frac{\mathsf{\Delta }t}{\mathsf{\Delta }x\left(\mathbf{x}\right)}·\left|\mathbf{U}\right|\left(\mathbf{x}\right)\right).\end{array}$$

(7)

We also use the acoustic Courant number, where the velocity $\left|\mathbf{U}\right|\left(\mathbf{x}\right)$ is replaced by $(c+\left|\mathbf{U}\right|)\left(\mathit{x}\right)$, with c representing the speed of sound. The implicit solver allows for acoustic Courant numbers $\gg 1$. The time step size is adaptive, and a supervisor script controls the maximum allowed values of the Courant criteria over time. Details are given in Appendix B.

2.2. Vorticity

The vorticity $\mathsf{\Omega }$ is defined as the curl of the velocity field [52]:

$$\begin{array}{c}\hfill \mathsf{\Omega }=\nabla \times \mathbf{U}.\end{array}$$

(8)

Vorticity maps are given for selected time instants together with flow-field direction arrows to characterize the vortex strength in the liquid. Due to the constraint of axial symmetry, ${\mathsf{\Omega }}_{\mathrm{axisymm}}\propto {\mathbf{e}}_{\theta }$, with ${\mathbf{e}}_{\theta }$ representing the unit vector in the azimuthal direction. When taking the cut through the bubble on the plane $z=0$, only the value of the z-component is plotted:

$$\begin{array}{c}\hfill {\left.{\mathsf{\Omega }}_{\mathrm{axisymm}}\right|}_{z=0}={\left.\left(\frac{\partial U^{\mathit{x}}}{\partial y}-\frac{\partial U^y}{\partial \mathit{x}}\right)\right|}_{z=0}{\mathbf{e}}_z.\end{array}$$

(9)

2.3. Selection of Time Instants for Vortex Visualization

In order to have an idea of the bubble volume versus time (bubble oscillation) compared to the time of vortex generation and existence, the equivalent radius, $R_{\mathrm{equiv}}$, versus time is given for the normalized distance of the bubble to the boundary $D^{*}=0.4$ and 1.8 in Figure 2. The equivalent radius is defined as $R_{\mathrm{equiv}}=$ ${\left(\frac{3}{4\pi }{V}_{\mathrm{bubble}}\right)}^{1/3}$, where $V_{\mathrm{bubble}}$ is the corresponding volume of the bubble.

It is seen in Figure 2 that with the bubble model applied (Equations (1)–(6)), the bubble exhibits rebounds with different strengths. However, in any case, the bubble oscillations have ceased after about 550 µs. Nevertheless, the liquid flow induced persists and may form liquid flow structures. The most prominent ones are liquid vortices of different kinds. As the vortices fully develop only long after the bubble oscillations have ceased, most vortex plots are given at $t_{\mathrm{map}}$ = 550 µs after bubble generation. To decide about the type of vortex in critical cases, the vortices at $t_{\mathrm{map}}$ = 1000 µs are given. The decision on the vortex type at 1 ms is called the 1 ms criterion.

2.4. Flow-Field Visualization

As indicated in Section 1, there exist different methods to experimentally visualize liquid flows. With the PTV and/or PIV method, not only the flow can be visualized, but also the velocity field can be measured. The particles used in PTV and PIV can either be reflective or fluorescent. Their location and motion is followed photographically or holographically. For reflecting particles applied in the flow around bubbles, see, e.g., Lauterborn and Vogel [11], Vogel et al. [13] and Lauterborn and Kurz [53], and for fluorescent particles, see, e.g., Kröninger et al. [26] and Reuter et al. [30]. From consecutive illumination flashes, velocity vectors can be derived from the illuminated particles.

A more compact method is dye advection applied to experiments, where the liquid is at rest in the beginning: An ink-colored liquid layer is put into the liquid sample at the start of the measurement, which is then distorted by the flow to give a flow-field map when photographed at times wanted.

When post-processing the data of a time-varying velocity field, gained by either experiment or numerical simulation, the dye advection method can be set up a posteriori by digitally coloring the liquid and advecting the color with the corresponding, calculated tracer trajectories. When the dye is injected also at different times $t\ne 0$, this is termed dye injection [31]. In the case of a laser-induced cavitation bubble it is, however, more suitable to inject the color only at time zero, since the dynamics is very fast.

A similar method is used here with digitally colored liquid layers, stratified parallel to the flat surface of the solid in a quiet liquid at bubble generation, whereby the liquid flow is not obtained from measurements but from numerical calculations. The numerical dye or, in OpenFOAM-terminology, passive scalar $\mathsf{\Gamma }$ is advected passively with the velocity field by an advection equation [51]:

$$\begin{array}{c}\hfill \frac{\partial \mathsf{\Gamma }}{\partial t}+\nabla ·\left(\mathsf{\Gamma }\mathbf{U}\right)-\left(\nabla ·\mathbf{U}\right)\mathsf{\Gamma }=0.\end{array}$$

(10)

$\mathsf{\Gamma }$ is initiated varying only in direction orthogonal to the solid surface and thus as stratified color layers when applying an adequate color map to the values of $\mathsf{\Gamma }$. The validity of the method was investigated in Koch [8]. The initial data of $\mathsf{\Gamma }$ are set to:

$$\begin{array}{c}\hfill \mathsf{\Gamma }=\left\{\begin{array}{cc}\frac{{\displaystyle y}}{{\displaystyle 1.5·R_{\max }^{\infty }}}\hfill & y<1.5·R_{\max }^{\infty }\hfill \\ 1\hfill & \mathrm{else}\hfill \end{array}\right.,\end{array}$$

for the center of the bubble in $\mathit{x}=(0,0)$at $t=0$. With these initial data, the color layers extend from the solid surface, $y=-D_{\mathrm{init}}$, to $y=1.5·R_{\max }^{\infty }$. For different $D_{\mathrm{init}}$, hence different $D^{*}$, the thickness of the color layers is adapted.

Figure 3 shows an example of one type of a liquid vortex after the expansion, collapse and rebounds of a bubble close to a solid boundary. The vortex is visualized from the simulation with the dye advection method via Equation (10). The stratified colors from bottom to top are: whitish, dark blue, black, green, light blue and purple. The bubble oscillation has long gone, but the liquid set into motion still persists and forms the large vortex shown at 550 µs after bubble generation. The initial bubble is modeled by a small sphere of high pressure, which is given by the small, white sphere at the center of the large, transparent sphere in the cut through the axis of symmetry. The transparent sphere with radius $R_{\max }^{\infty }$is copied in to show the maximum extension the bubble would reach in an unbounded liquid. It gives an impression of the normalized distance $D^{*}$ (smaller than one, here $D^{*}$ = 0.8) by crossing the solid boundary and lets us imagine that the actual bubble must have undergone a (rotationally symmetric) distortion during expansion by the existence of the solid surface.

For this type of diagram, also, the motion direction of the vortex can be inferred from the color it carries. In the given case of $D^{*}=0.8$, the motion is upwards for the large vortex. The near-boundary liquid is carried upwards into the free liquid. The colors in the large ring vortex from inside of the vortex ring to outside are dark blue, black, green, light blue and violet.

Figure 4 shows the same vortex as in Figure 3 but in a different way with the vorticity $\begin{array}{c}\hfill {\mathsf{\Omega }}^z=\end{array}(\nabla \times \mathit{U})·{\mathbf{e}}_z$ (z-direction pointing toward the reader) and superimposed flow-field direction arrows to show the actual flow field at $t_{\mathrm{map}}=550\mu \mathrm{s}$. In this case, the flow of the large vortex is upwards (away from the solid boundary) at the axis of symmetry and redirected back downwards at larger distances from the axis to form a large-scale, dipole-like vortex flow. Due to the axial symmetry, the rotation is clockwise (in this example) in the right half of the cut and anti-clockwise in the left half. From the flow-field direction arrows and the vorticity colors, it is seen that there are several small-scale ring vortices present at the solid boundary and one with opposite turning direction between the solid surface and the large-scale vortex.

The two forms of flow-field visualization as presented in this section will be used in the following to characterize the vortices found in laser-generated bubbles. They give different insights as to the properties of the vortices. The dye advection maps incorporate the flow up to the plot time, and the vorticity maps with superimposed flow-field direction arrows are a snapshot of the vortex at plot time. Thus, they carry different information and add to one another.

3. Results

The color-layer mixing dynamics of the dye advection is explained by means of a time series. Afterwards, two types of vortices in the case of a bubble close to a solid boundary are introduced with their different rotation and translational motion directions (Section 3.2). They coexist throughout the interval of $D^{*}$ investigated. In a transition interval, they exchange their dominating strength, as given by the modulus of the vorticity, |$\nabla \times \mathbf{U}$|, and the spatial extent. As the settling time for the flow field becomes larger, when approaching the transition threshold, the respective simulations were extended to 1 ms since bubble generation. Finally, the numerically found transition interval is compared to the experimentally determined one by Reuter et al. [30] (Section 3.4).

3.1. Color-Layer Mixing Dynamics

The color mixing—as seen in Figure 3 for the bubble at $D^{*}=0.8$ at $t_{\mathrm{map}}$ = 550 µs—shows a complicated structure explained as a large vortex with additional details. The green needle along the axis of symmetry, for instance, can be attributed to the motion of the bubble and its surrounding liquid from the green layer toward the solid and the squeezing back of the green liquid toward the axis of symmetry in the course of collapses and rebounds. This will become clear when the development of the flow is depicted in a time series of dye advection maps in Figure 5. The first bubble collapse produces a liquid microjet toward the solid. The jet then drives liquid toward the solid boundary. Liquid from the green layer around the axis of symmetry is transported downwards by the jet (frames 2–6, Figure 5). Liquid from the light blue and black color layers are sucked into the jet as well. In frame three of Figure 5, a ring sheet (in 3D) of liquid shoots into the bubble. This phenomenon, termed annular nanojet, has been described in Lechner et al. [43] and the reader is referred to this reference for details. After impact at the solid surface, liquid from the jet spreads along the surface, as seen in frame 4. The reason as to why the final motion of the liquid is directed into the liquid bulk is the inversion of motion by the second bubble collapse and only small rebound (frames 5–9). The second collapse (frame 5) starts at the outer rim of the gas torus, and the flow is directed inwards toward the axis of symmetry (frames 6 and 7). The rebound (frame 7) of the gas phase then takes place from near the solid boundary with a layer of liquid in between. The second rebound is of small amplitude (see Figure 2) and is not able to stop the inflow of liquid from the second collapse. The former downwards-oriented part of the light-blue and green layer is now squeezed to the axis by the inward liquid motion from the second collapse and necessarily pushed upwards together with all layers except the whitish one. The whitish layer has been pushed so far away from the axis of symmetry by the surface flow from other layers that it no longer takes part in the essential dynamics.

Additionally, several small-scale vortices are seen along the surface of the solid (see also Figure 4, the vorticity diagram). Prominently seen is the small bubble vortex following the jet impact during first collapse. It spreads along the surface by the ongoing inflow of liquid (frame 4 at time 109.61 µs in Figure 5). The complete set of vortices found is described next.

3.2. The Two Types of Vortices and Their Origins

In axial symmetry, there are two types of bubble-jet-induced liquid vortices close to a solid boundary. Either the liquid flow through the vortex center along the axis of symmetry is upwards, away from the boundary—“up-vortex” (motion of the forming vortex also upwards as can be guessed from the color advection or seen from a succession of frames), or downwards, toward the solid boundary—“down-vortex” (motion of the forming vortex downwards as demonstrated later). An example of an up-vortex is given in Figure 6, and an example of a down-vortex is given in Figure 7. The figures are a combined version of a flow and a vorticity map, as shown separately in Figure 3 and Figure 4. Because of the axial symmetry, a half part of the cut is sufficient.

The two types of vortices can further be distinguished by the second type of visualization, as presented in Figure 4 by a toroidal rotation criterion. When in the cut only the right half part of the flow-direction map is considered, then the flow around an up-vortex is “clockwise”. Figure 6 shows an example of a vortex with clockwise rotation of the liquid, i.e., an up-vortex. Note that the rotation direction is only valid for the right part of the cut through the torus vortex. It is known that a ring vortex in an unbounded liquid propagates in the direction of the flow along its axis (see, e.g., Fraenkel [54]). When the vortex axis coincides with the axis of symmetry in an axially symmetric setting, then for a single vortex, the translational motion is upwards—using the terms introduced in the present work—when rotating clockwise, hence downwards when rotating counter-clockwise. However, when vortices come in larger numbers and are densely packed, the rotation criterion is best suited to characterize the type of vortex in an axially symmetric setting.

In Figure 7, an example of a down-vortex is given with its counter-clockwise rotation for $D^{*}$ = 1.4. It is located near the solid surface. A small, locally confined region with clockwise rotation is seen at the axis of symmetry near the location of bubble generation (the small white circle), too.

The temporal dynamics of the down-vortex generation are shown in the montage in Figure 8 for a bubble at $D^{*}=1.4$. The main difference in the dynamics compared to the ones for the up-vortex in Figure 5 is the sequence of the bubble shape in the first rebound and second collapse. Here again, flow focusing isthe main driving force [5,55]. Flow focusing is the phenomenon generating the strongest acceleration where interface curvatures are highest. When comparing frame 5 ofFigure 5 to frame 6 of Figure 8, it is noted that the curvature of the torus bubble during the second bubble maximum volume and onset of the second collapse is highest at the outer rim for the up-vortex case, and the curvature is highest in the inner region for the down-vortex case. Therefore, the fluid motion is accelerated highest toward the axis of symmetry in case of the up-vortex, and it is accelerated outwards in case of the down-vortex during the second bubble collapse.

The vortices are dynamical objects that may move and/or change their form. Up-vortices move upwards away from the solid boundary, while down-vortices move toward the solid boundary. Reuter et al. [30] use the motion criterion for classification and call the upwards moving vortices “free vortices” and the downwards moving vortices “wall vortices”. The authors introduced a third status as well, tagged “no vortex”, when neither an up-vortex nor down-vortex was observed within experimental uncertainties. However, in the numerical simulations presented here, it becomes evident that the vortices are always present after the bubble oscillations have ceased. The difference is only in the intensity of the vorticity. In fact, both types of vortices are generally present. It is the ratio of their respective intensities that determines the macroscopic diagnosis made by an observer. An example is given in Figure 9 for $D^{*}=0.4$. The clockwise rotation direction indicates an up-vortex, while the vortex does not seem to have moved upwards up to $t_{\mathrm{map}}$ = 550 µs but is located away from the boundary. There, the flow rather rolls up, as seen in the left part of the plot. The extension to $t_{\mathrm{map}}$ = 1 ms indicates a vortex with very slow motion upwards, confirming the classification of an up-vortex by its rotation direction from one frame only.

3.3. Survey of Vortex Structures in Dependence on $D^{*}$

The aim is to find the (dominating) vortices for 0.4 $<D^{*}<$ 1.8 and whether they are up-vortices or down-vortices. The vortices for the values $D^{*}=0.4,0.8$ and 1.4 were already given in Figure 6, Figure 7 and Figure 9. The flow maps for the $D^{*}$-values of 0.6, 1.0 and 1.2 are given in Figure 10, Figure 11 and Figure 12, respectively, and in Figure 13 for $D^{*}$ = 1.6 and 1.8. It is seen that up to $D^{*}$ = 1.2, the up-vortices dominate (clockwise rotation). For the values of $D^{*}$ higher than or equal to 1.4, the turning direction of the dominating vortex changes by fading of the up-vortex and by the increasing in strength and extent of the down-vortex. The result is a down-vortex with its counter-clockwise rotation and a location close to the solid surface (Figure 7 and Figure 13).

3.4. The Transition Interval

To better locate the change in the vortex rotation direction from clockwise (up-vortex) to counter-clockwise (down-vortex), the interval 1.2 $<D^{*}<$ 1.4 was scanned in steps of $\mathsf{\Delta }{D}^{*}=0.01$. The result is presented in Figure 14, Figure 15, Figure 16 and Figure 17 for the transition interval from dominating clockwise rotation to dominating counter-clockwise rotation. The simulations have been extended to $t_{\mathrm{map}}$ = 1 ms, as a unique rotation direction obviously needs time to settle near the transition threshold (1 ms criterion). It is seen that the dominating rotation direction is clockwise at $D^{*}$ = 1.27 and up to $D^{*}$ = 1.29 and counter-clockwise at $D^{*}$ 1.35 and for further values of $D^{*}$. A quite sharp transition interval$D_{\mathrm{tr}}^{*}$ is found between $D^{*}$ = 1.30 and 1.34. At $D^{*}=1.31$, it seems that both rotation directions are of equal strength, and thus, there is no overall flow field rotation, and the vortex region is locally confined.

Up- and down-vortices of about equal strength coexist. As $D^{*}$ increases, the up-vortices gradually become slower in their upward motion. A closer inspection of the vortex structure reveals that a coexistence of up- and down-vortices is always the case for the $D^{*}$-values studied. It is the large-scale rotation direction that changes between $D^{*}$ = 1.30 and 1.34. The region of the fading up-vortex has a surrounding clockwise flow field of fading extent superimposed on the counter-clockwise overall flow field from $D^{*}=1.34$ onwards.

3.5. Comparison with Experiments

The numerical findings are compared to the experimental findings of Reuter et al. [30] in Figure 18 and Figure 19.

Reuter et al. [30] use a different definition for the normalized distance, $\gamma$, from the one used here, $D^{*}$. In the experiment, it is impossible to obtain the maximum radius in the unbounded liquid of the bubble just generated if the current setup contains restrictive objects or dimensions. Therefore, instead of using $R_{\max }^{\infty }$ (refer to Equation (1)) corresponding to the maximum radius of the same bubble in an unbounded liquid, the apparent maximum radius from the measurement was taken, resulting in a $R_{\max }^{\mathrm{eq}}$, meaning the radius of a sphere with equivalent volume. The precise algorithm of how $D_{\mathrm{init}}$ and $R_{\max }^{\mathrm{eq}}$ were extracted from the measurement data is not clear from the publication. Here, it was decided to plot the data from Reuter et al. [30] as if the relation $\gamma \approx D^{*}$ would hold, which is assumed to be justified for $D^{*}\ge 0.4$, where fast-jet dynamics are not present yet. The data points from the experiment are plotted with a larger diameter to account for the uncertainty of that assumption.

Reuter et al. [30] discern between three types of vortices: the free vortex, the wall vortex and “no vortex” according to their translational motion. No vortex then means no discernible vortex motion. Here, only two categories of vortices are introduced: the up-vortex and the down-vortex according to their rotation direction. Usually, a free vortex is an up-vortex and a wall vortex is a down-vortex. There is no correspondence to a no-vortex. Using the derived knowledge from the present work, the terminology was adapted to the notation of the rotation direction, up- or down-vortex and whether the vortex propagation velocity nearly vanishes (“several, locally confined”).

When comparing the experimental data points with the numerical ones, an excellent agreement is found. In particular, the numerical transition interval ($1.30<D_{\mathrm{tr}}^{*}<1.34$) fits perfectly into the experimental transition interval $1.30<D_{\mathrm{tr}}^{*},\exp <1.36$, where Reuter et al. [30] could find no definite signs for a free or wall vortex.

4. Discussion

An oscillating bubble is an object with long-range interaction forces. In the spherically symmetric case of a free bubble in an infinite, incompressible liquid, it is the breathing mode of the oscillation that leads to a spherically symmetric, quadratically decaying liquid motion. Thus, in principle, the interaction range is infinite. For numerical calculations in single-bubble dynamics, this means that the boundaries must be sufficiently far away from the bubble. It has been found that the boundaries should be at least about 100 times the maximum extension of a spherical cavitation bubble away in order to not significantly alter the dynamics [37] and for numerical convergence reasons [8].

In the non-spherical case with boundaries nearby, the dynamics is totally different with its liquid jet and vortex formation. However, even in a free liquid, jet formation occurs when the liquid is not in a spherically symmetric condition, as for instance with pressure gradients in a gravitational field [56,57] or in a standing acoustic field [58]. As non-spherical conditions prevail in nature and applications, the study of the flow in non-spherical cases is of general interest. There are an infinite number of non-spherical cases. The most important one from a technical application standpoint (cleaning and erosion) is the case of a bubble near a surface or surfaces with the special case of a plane solid surface. There are only a few experimental studies on the flow in this case, in particular on the flow vortices. The, at the time being, most advanced experimental study on vortex formation is from Reuter et al. [30]. Since then, numerics has improved significantly and now is capable of modeling and simulating practically any non-spherical case.

In the present work, the liquid motion of a bubble near a plane solid surface is studied numerically for several distances of the bubble from the solid boundary. A bubble with a radius of 500 µm in an unbounded liquid has been chosen. The focus is on the vortices—generated by the liquid jet, torus bubble collapse, rebound and the further decaying oscillations—and their long-term behavior with respect to translational motion and rotation direction. The experimentally found two types of vortices could be confirmed numerically, which are here called the up-vortex and down-vortex. Both types are present simultaneously with one type being dominant in strength and extent. In a transition interval, one type ceases and the other type takes over. The transition interval is located at $1.30<D_{\mathrm{tr}}^{*}<1.34$. It lies inside the experimentally found transition interval of $1.30<D_{\mathrm{tr},\exp }^{*}<1.36$ [30]. The two vortex types can be discerned by the rotation direction of the flow (see Section 3.3 and Figure 6 and Figure 7). The flow needs some time to develop and settle to a dominant type. Therefore, the numerical transition interval was determined at 1 ms after bubble generation, i.e., a long time after jet formation compared to the collapse time of about 50 µs (see Figure 2).

The origin for the two main vortex propagation directions, upwards and downwards, has been determined to be flow focusing during the second collapse (refer to Section 3.3 and Figure 5 and Figure 8). High bubble interface curvatures are accelerated faster than low curvatures. The torus–bubble cross-sectional shape (volume) after rebound determines the fluid acceleration direction afterwards.

It is assumed that there must be a connection to the dynamics of the bubble and therewith to the erosion of the surface. The erosion pattern strongly alters with the normalized distance in a characteristic manner as does the bubble dynamics. Both show a strong correlation: erosion is only observed for bubbles collapsing directly in contact with the solid. In the interval $1.2<D^{*}<1.4$ around the transition, the erosion is strongly increased (Figure 21 [59]). Already at $D^{*}\approx 1.1$ and higher values, the first collapse does not take place directly at the solid boundary anymore; however, in this aforementioned interval, the second collapse is from a torus bubble directly sitting on the surface. Correspondingly, the erosion pattern is ring-like as observed. Wang et al. [60], with a modified BIM model (boundary integral method), was able to numerically follow the bubble oscillations to the second collapse, when it becomes a torus bubble sitting on the surface. The torus size coincided with the damage ring found by Philipp and Lauterborn [59]. Moreover, at $D^{*}$ = 1.2, the bubble lifts off the surface after the second collapse, and the third collapse takes place away from the surface just like the first collapse before and like the following ones (Figure 2f [59]). This is observed in the simulation of the present work for the same $D^{*}$, too. This is in correspondence with the final up-vortex, here carrying a larger gas fraction with it.

Now, when looking at the bubble dynamics in the experiment for $D^{*}$ = 1.4, it is seen that after the second collapse, the torus bubble does not lift up but rather stays at the surface and spreads there. Thus, the flow of liquid connected with it is downwards, and it is in correspondence with the down-vortex. The transition $1.30<D_{\mathrm{tr}}^{*}<1.34$ is well in between 1.2 < $D^{*}$ < 1.4. Tomita and Shima [61] measured the maximum collapse pressure for the second collapse at $D^{*}$ = 1.3. This is an extremely interesting experimental finding in view of the numerical findings here. At $D^{*}$= 1.3, the transition from a dominating up-vortex to a dominating down-vortex takes place. Simultaneously, Philipp and Lauterborn [59] found the maximum erosion around this $D^{*}$. Thus, there seems to exists a correlation between all four quantities/events: maximum pressure, maximum erosion, second collapse (directly at the solid surface for $D^{*}$= 1.3), and the transition. The erosion is increasingly diminished for growing $D^{*}$> 1.4 by the growing distance of the relevant events (jet impact pressure and bubble collapse pressure) from the solid surface. At about $D^{*}$= 2.0, erosion ceases.

At first sight, it could be expected that the jet flow dominates the more the nearer the bubble collapses to the surface, i.e., at small $D^{*}$. However, the example of $D^{*}=1.2$ shows that it is not the case: the final flow is upwards. The reason is that there is a competition between the jet flow from the first collapse and the subsequent rebound and second collapse. The time when the bubble touches the solid boundary is always between the first bubble volume minimum and the second maximum for $D^{*}>1.1$ in the investigated interval of $D^{*}$. The higher the value for $D^{*}$, the later the bubble touches the solid boundary during the first rebound. It seems to matter how much of the rebound momentum takes place along the solid boundary in order to create a flow reversal during second collapse against the inflow from the jet from the first collapse.

For values below $D^{*}$ = 0.4, the transition to the fast jet is entered with the fast-jet region below $D^{*}$≈ 0.2. These values are left for further studies in view of their special dynamics and jet production mechanism.

5. Conclusions and Outlook

Vortex dynamics following the jet formation of a bubble expanding and collapsing close to a plane solid boundary is studied numerically. Two types of vortices, up- and down-vortices, exist simultaneously with one type dominating in strength and extent. The interval of $D^{*}$ from 0.4 to 1.8 has been studied. In the interval$1.30<D^{*}<1.34$, a transition from one type to the other takes place with the up-vortices below and the down-vortices above the transition value. The numerically found transition value lies well in the experimentally found transition interval.

The benefit of the comparison of the numerical results to the experimental results is threefold. Firstly, the stability of the solver with its extensions (see Appendix A–Appendix C) is demonstrated for long calculation times. Secondly, the perfect agreement with the experiments validates the solver for the cases studied, in particular long-term vortex formation. Thirdly, it is shown that the main bubble and vortex dynamics are essentially axially symmetric, even if this is not evident from the experiments (2D projections).

Moreover, by comparison of different works, a correlation between erosion, bubble dynamics, pressure at the solid and vortex dynamics has been established in the range 0.4 < $D^{*}$ < 1.8, which is valid for a bubble with a maximum radius of 500 µm. Only bubbles collapsing in direct contact with the solid cause erosion.

With these results, new tasks come into reach: for instance, long-term studies of acoustic bubble oscillations with extreme bubble surface deformations near solids. Further bubble mass control methods may be necessary for this case. Then, studies as to the growth to a stable state or to instability in dependence on the acoustic pressure amplitude can be conducted and compared with experiments. Long-term bubble interaction studies may profit from the knowledge of the large spatial extent of the liquid vortex motion and the longevity of the vortices. Thus, the present work opens up a large area for future research.