Morphological Analysis of a Collapsing Cavitation Bubble near a Solid Wall with Complex Geometry

Featured Application

Collapsing cavitation bubbles can be used in material surface cleaning, the medical field, and so on. By adjusting the micro-jet intensity of the collapsing bubbles, the cavitation phenomenon can be employed to clean irregular material surfaces, such as sections, cracks, and vegetable leaves. In the medical field, cavitation bubbles can be used as microbubble contrast agents for ultrasound diagnostic imaging or vehicles for drug or gene delivery. The growth and violent collapse of cavitation bubbles can also be employed in sterilization or killing viruses such as COVID-19.

Abstract

The interaction mechanism between the cavitation bubble and a solid wall is a basic problem in bubble collapse prevention and application. In particular, when bubble collapse occurs near solid walls with arbitrarily complex geometries, it is difficult to efficiently establish a model and quantitatively explore the interaction mechanism between bubbles and solid walls. Based on the advantages of the lattice Boltzmann method, a model for cavitation bubble collapse close to a solid wall was established using the pseudopotential multi-relaxation-time lattice Boltzmann model. Solid walls with arbitrarily complex geometries were introduced in the computational domain, and the fractal dimension was used to quantify the complexity of the solid wall. Furthermore, owing to the lack of periodicity, symmetry, spatial uniformity and obvious correlation in this process, the Minkowski functionals-based morphological analysis method was introduced to quantitatively describe the temporal evolution of collapsing bubble profiles and acquire effective information from the process. The interaction mechanism between the bubble and solid wall was investigated using evolutions of physical fields. In addition, the influences of the solid walls’ surface conditions and the position parameter on collapsing bubbles were discussed. These achievements provide an efficient tool for quantifying the morphological changes of the collapsing bubble.

1. Introduction

Cavitation occurs when the local pressure becomes less than a liquid’s saturated vapor pressure. The bubble collapse that happens near a solid wall can severely damage the neighboring solid surface owing to the resulting high velocities, pressure, and temperature [1]. In addition, the bubble experiences inward deformation, micro-jet formation, collapse, or multiple collapses owing to the retarding effect of the solid wall. The potential applications of cavitation have been extensively explored in material surface cleaning [2], ultrasonic therapy [3], and environmental protection [4].

The interaction mechanism between a cavitation bubble and a solid wall has an important role in various applications. In particular, when cavitation occurs near a solid wall involving a complex geometry, the profile evolution of the collapsing bubble can be more varied. The effect of the surface geometry of the solid wall on the dynamic behavior of the collapsing bubble has been experimentally investigated [5,6,7]. A complex boundary can significantly alter the dynamic behavior of cavitation bubbles, such as the bubble shape and micro-jet. However, it is difficult to investigate the interaction mechanism between bubbles and solid walls with arbitrarily complex geometries because of the limitations of experimental methods, e.g., experimental conditions cannot be changed at any time. Numerical methods have been adopted to efficiently simulate cavitation bubble collapse near a solid wall and investigate the dynamic bubble collapse process. Conventional numerical simulation methods include the finite volume method (FVM) [8,9], the boundary element method (BEM) [10], and the finite element method [11]. Yu et al. [9] investigated the dynamic process of a vapor bubble collapse near a rectangular groove wall using the FVM and the volume of fluid (VOF) method. Saleki-Haselghoubi and Dadvand [10] adopted the BEM to investigate the dynamics of a spark-generated bubble that oscillated near a circular aperture in a curved rigid plate. However, these macroscopic numerical methods require other approaches to track or capture gas–liquid interfaces, such as the VOF [12] and level set method [13], which may reduce computational efficiency. In implementing the complex geometry boundaries, it is difficult for these macroscopic models to choose approximate functions in each discrete unit to characterize changes in arbitrarily complex boundaries.

Over the past few decades, the lattice Boltzmann method (LBM), has been developed as a powerful tool for simulating multiphase flow problems [14,15,16,17,18,19]. The LBM based on the discrete Boltzmann equation with the simplified collision operator can model phase segregation and interfacial dynamics in multiphase flows with collision and streaming steps [15]. The mesoscopic approach has multiple unique advantages such as high efficiency, completely parallel algorithms, and the simplicity of programming. Generally, the multiphase models in the LBM inlude four kinds of different models, i.e., the color-gradient [20], pseudopotential [21,22,23], free-energy [24], and phase-field models [25]. The pseudopotential model proposed by Shan and Chen, is also named as Shan–Chen (SC) model. In this model, an interparticle potential is proposed to simulate the interactions between the fluids, from which a non-monotonic equation of state (EOS) can be obtained. As a result, the interface can be formed automatically, and the pressure can be calculated from the EOS efficiently. Sukop and Or [26] first introduced the Shan–Chen LBM into the cavitation field and concluded that applications to more complex problems should be possible. Later on, several scholars devoted themselves to investigating the cavitation phenomenon with the pseudopotential LBM [27,28,29,30,31]. Meantime, the LBM has good spatiotemporal flexibility in dealing with boundary conditions, which enables the easy introduction of solid boundaries with complex geometries into numerical models. Shan et al. [32,33] simulated the collapse process of a single bubble near a plane solid wall and between two parallel solid walls. Mao et al. [34] and He et al. [35] both used a multi-relaxation-time (MRT) pseudopotential model to simulate dual-bubble collapse near a rigid boundary. It can be found that the relative position of two bubbles plays an important role in the bubble near-wall collapse process. Xue et al. [36] simulated the cavitation bubble near a curved wall using the MRT-SC model and discussed the influence of bubble–wall distance on bubble collapse. Shan et al. [37,38] investigated the influence of a rectangular bulge solid wall and hydrophobic concave corner on cavitation bubble collapse from the evolution of physics fields. The geometry of the solid wall significantly influences bubble collapse. However, previous studies have considered the complex geometries of solid walls as periodic shapes with equal widths [9,37]. This cannot account for certain characteristics of actual complex solid walls, such as sections and cracks. Thus, it is necessary to accurately model actual complex solid walls.

In contrast to the solid walls in previous works, the arbitrarily complex geometry solid walls have more complex geometries. To reflect the differences between the solid walls, the fractal dimension is adopted to quantitatively describe the complexity of a solid wall, and it increases with the complexity. Various processes occur when a cavitation bubble collapses near a solid wall with a complex geometry. Further research is required on the methods for describing the cavitation bubble near-wall collapse process and the effect of the solid wall on the collapsing bubble considering the large amount of data obtained from the bubble collapse process, various complex configurations, and evolutions of physical fields. The Minkowski functionals-based morphological analysis method (Min-MAM) [39] has been introduced to investigate complicated physical fields. Sofonea and Mecke [40,41] introduced the Minkowski functionals to describe the morphological evolution of homogeneous phases during spinodal decomposition. Gan et al. [42] characterized the density and velocity fields of the phase separation process in a complex fluid using the Minkowski functionals. They clarified the effect of the temperature field on phase morphology and spinodal decomposition. Xu et al. [43] conducted a detailed review of the Minkowski functionals and adopted them to explore the complex procedure of shock wave reaction on porous material. The Min-MAM is an efficient tool for analyzing and acquiring effective information in complicated physical fields [44,45,46]. Therefore, the Min-MAM can potentially be adopted to quantify the profile evolution of collapsing bubbles.

In this study, we introduce a solid wall with an arbitrarily complex geometry and investigate the dynamic bubble collapse process. The pseudopotential MRT-LBM is adopted to establish the collapse model of a cavitation bubble in the vicinity of a solid wall with an arbitrarily complex geometry. Furthermore, the fractal dimension and Min-MAM are used to quantitatively describe the complexity of the solid wall and the evolution of bubble morphology. The interaction mechanism between a cavitation bubble and a solid wall is investigated by changing the surface complexity, geometric features of the solid wall, and position offset parameter. In general, this study systematically investigates the effect of solid walls with arbitrarily complex geometries on the bubble collapse process. The morphological evolution of the collapsing bubble is described with quantitative measures, which provides an alternative tool for analyzing the cavitation bubble near-wall collapse process. The remainder of this paper is organized as follows. The pseudopotential MRT-LBM and Min-MAM are presented in Section 2. The numerical simulation model and dynamic bubble collapse process are described in Section 3. The morphological characterization of the collapsing bubble in the vicinity of a solid wall with a complex geometry is discussed in Section 4. Finally, the conclusions are presented in Section 5.

2. Methods

2.1. Pseudopotential MRT-LBM

The pseudopotential MRT-LBM is adopted to simulate the flow field. By using the MRT collision operator, the evolution of the density distribution function can be expressed as [47]

$$f_{\alpha }\left(\mathbf{r}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)=f_{\alpha }\left(\mathbf{r},t\right)-{\left({\mathbf{M}}^{-1}\mathbf{\Lambda }\mathbf{M}\right)}_{\alpha \beta }\left(f_{\beta }-f_{\beta }^{eq}\right)+{\delta }_t{F^{\prime }}_{\alpha }$$

(1)

where $f_{\alpha }(\mathbf{r},t)$ is the particle density distribution function, $f_{\beta }^{eq}(\mathbf{r},t)$ is the corresponding equilibrium distribution function, $\mathbf{r}$ is the particle position, ${\delta }_t=1$ is the time step, ${\mathbf{e}}_{\alpha }$ is the discrete velocity along the $\alpha \mathrm{th}$ direction, ${F^{\prime }}_{\alpha }$ is the component of forcing term $F^{\prime }$ along the $\alpha \mathrm{th}$ direction, $\mathbf{M}$ is an orthogonal transformation matrix, and $\mathbf{\Lambda }$ is the relaxation diagonal matrix.

For the D2Q9 lattice model, ${\mathbf{e}}_{\alpha }$ is defined as [48]

$${\mathbf{e}}_{\alpha }=\{\begin{array}{ll}\left(0,0\right)& \alpha =0\\ \left(\cos \left[(\alpha -1)\frac{\pi }{2}\right],\sin \left[(\alpha -1)\frac{\pi }{2}\right]\right)& \alpha =1,2,3,4\\ \sqrt{2}\left(\cos \left[(\alpha -5)\frac{\pi }{2}+\frac{\pi }{4}\right],\sin \left[(\alpha -5)\frac{\pi }{2}+\frac{\pi }{4}\right]\right)& \alpha =5,6,7,8\end{array}$$

(2)

$\mathbf{M}$ is given by

$$\mathbf{M}=\left(\begin{array}{ccccccccc}1& 1& 1& 1& 1& 1& 1& 1& 1\\ -4& -1& -1& -1& -1& 2& 2& 2& 2\\ 4& -2& -2& -2& -2& 1& 1& 1& 1\\ 0& 1& 0& -1& 0& 1& -1& -1& 1\\ 0& -2& 0& 2& 0& 1& -1& -1& 1\\ 0& 0& 1& 0& -1& 1& 1& -1& -1\\ 0& 0& -2& 0& 2& 1& 1& -1& -1\\ 0& 1& -1& 1& -1& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& -1& 1& -1\end{array}\right)$$

(3)

and $\mathbf{\Lambda }$ is expressed as

$$\mathbf{\Lambda }=diag\left({\tau }_{\rho }^{-1},{\tau }_e^{-1},{\tau }_{\xi }^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_j^{-1},{\tau }_q^{-1},{\tau }_{\nu }^{-1},{\tau }_{\nu }^{-1}\right)$$

(4)

where ${\tau }_{\rho }^{-1}$, ${\tau }_e^{-1}$, ${\tau }_{\xi }^{-1}$, ${\tau }_j^{-1}$, ${\tau }_q^{-1}$, and ${\tau }_{\nu }^{-1}$ are the relaxation time factors for different moments. In the simulation, these relaxation times are selected as ${\tau }_{\nu }^{-1}={\tau }_{\rho }^{-1}={\tau }_j^{-1}=1.0$, ${\tau }_e^{-1}={\tau }_{\xi }^{-1}=0.8$, and ${\tau }_q^{-1}=1.1$.

Using $\mathbf{m}=\mathbf{M}f$ and ${\mathbf{m}}^{eq}=\mathbf{M}{f}^{eq}$, $f_{\alpha }$ and $f_{\beta }^{eq}$ can be projected onto the moment space. Equation (1) can be transformed to [49]

$$\mathbf{m}*=\mathbf{m}-\mathbf{\Lambda }\left(\mathbf{m}-{\mathbf{m}}^{eq}\right)+{\delta }_t\left(\mathbf{I}-0.5\mathbf{\Lambda }\right)\mathbf{S}$$

(5)

where $\mathbf{m}*=({\mathbf{m}}_0^{\ast },{\mathbf{m}}_1^{\ast },\dots ,{\mathbf{m}}_8^{\ast })$, $\mathbf{I}$ is the unit tensor, and $\mathbf{S}$ is the forcing term in the moment space with $(\mathbf{I}-0.5\mathbf{\Lambda })\mathbf{S}=\mathbf{M}{F}^{\prime }$. For the D2Q9 lattice, ${\mathbf{m}}^{eq}$ can be expressed as

$${\mathbf{m}}^{eq}=\rho {\left(1,-2+3|v{|}^2,1-3|v{|}^2,v_x,-v_x,v_y,-v_y,v_x^2-v_y^2,v_x{v}_y\right)}^T$$

(6)

where $\rho ={\displaystyle {\sum }_{\alpha }{f}_{\alpha }}$ is the macroscopic density, $\mathbf{v}$ is the macroscopic velocity and ${\left|\mathbf{v}\right|}^2=v_x^2+v_y^2$.

Then, the streaming process can be given as

$$f_{\alpha }\left(\mathbf{r}+{\mathbf{e}}_{\alpha }{\delta }_t,t+{\delta }_t\right)=f_{\alpha }^{*}\left(\mathbf{r},t\right)$$

(7)

where $f^{*}={\mathbf{M}}^{-1}{\mathbf{m}}^{*}$.

$\mathbf{v}$ in Equation (6) can be obtained from

$$v=\frac{{\displaystyle {\sum }_{\alpha }{\mathbf{e}}_{\alpha }{f}_{\alpha }+{\delta }_t\mathbf{F}/2}}{\rho }$$

(8)

where $\mathbf{F}=(F_x,F_y)$ is the total force acting on the fluid system in the two-dimensional (2D) space. $\mathbf{F}$ includes the fluid–fluid interactive, fluid–solid interactive, and bulk forces. According to the criteria described in Ref. [50], the influence of gravity can be neglected in the investigation of the cavitation bubble collapse problem. The fluid–fluid interaction force is adopted in the pseudopotential LBM, which is expressed as [51]

$${\mathbf{F}}_{\alpha }(\mathbf{r})=-G\psi (\mathbf{r}){{\displaystyle {\sum }_{\alpha }\omega }}_{\alpha }\psi \left(\mathbf{r}+{\mathbf{e}}_{\alpha }\right){\mathbf{e}}_{\alpha }$$

(9)

where $G$ is the overall interaction strength, $\psi (\mathbf{r})$ is the pseudopotential, and ${\omega }_{\alpha }$ is the weight coefficient, which is given by

$${\omega }_{\alpha }=\{\begin{array}{c}1/3,\\ 1/12,\end{array}\begin{array}{c}|{\mathbf{e}}_{\alpha }{|}^2=1\\ |{\mathbf{e}}_{\alpha }{|}^2=2\end{array}$$

(10)

In this study, $\psi (\mathbf{r})$ is considered as [51]

$$\psi (\mathbf{r})=\sqrt{\frac{2(p_{EOS}-\rho c_s^2)}{G{c}^2}}$$

(11)

where $p_{EOS}$ is the prescribed EOS, $c=1$ represents the lattice constant, and $c_s=c/\sqrt{3}$ is the lattice sound speed. The Carnahan–Starling (C–S) EOS is adopted according to Ref. [52], which is expressed as

$$p_{EOS}=\rho RT\frac{1+b\rho /4+{(b\rho /4)}^2-{(b\rho /4)}^3}{{(1-b\rho /4)}^3}-a{\rho }^2$$

(12)

where $a=0.4963{\left(R{T}_c\right)}^2/p_c$ and $b=0.1873{\left(R{T}_c\right)}^2/p_c$ with $T_c$ and $p_c$ being the critical temperature and pressure, respectively. In our simulation, $a=1$, $b=4$, and $R=1$ [52].

For the D2Q9 lattice, the forcing term $\mathbf{S}$ proposed by Li et al. is expressed as [49]

$$\mathbf{S}=\left[\begin{array}{c}0\\ 6(v_x{F}_x+v_y{F}_y)+\frac{0.75\epsilon |\mathbf{F}{|}^2}{{\psi }^2{\delta }_t({\tau }_{\mathbf{e}}-0.5)}\\ -6(v_x{F}_x+v_y{F}_y)-\frac{0.75\epsilon |\mathbf{F}{|}^2}{{\psi }^2{\delta }_t({\tau }_{\xi }-0.5)}\\ F_x\\ -F_x\\ F_y\\ -F_y\\ 2(v_x{F}_x-v_y{F}_y)\\ v_x{F}_y+v_y{F}_x\end{array}\right]$$

(13)

where $\epsilon$ is an adjustable coefficient used to tune the mechanical stability condition. In the present study, when $\epsilon$ is set as 1.86 in the case of a cavitation bubble, the coexistence curve simulated by the pseudopotential MRT-LBM agrees well with the solution of the Maxwell equal-area construction over a wide temperature range [33].

In addition, all physical quantities use the corresponding lattice units under the framework of the LBM, as shown in

Table 1. Basic physical quantity units.

Name	Variable	Unit
Time	${\delta }_t$	$\mathrm{ts}$
Length	${\delta }_x$	$\mathrm{lu}$
Mass	$m$	$\mathrm{mu}$
Density	$\rho$	$\mathrm{mu}/{\mathrm{lu}}^3$
Pressure	$p$	$\mathrm{mu}/\left({\mathrm{ts}}^2\mathrm{lu}\right)$
Velocity	$v$	$\mathrm{lu}/\mathrm{ts}$

.

2.2. Min-MAM

The Min-MAM is an effective quantitative method for characterizing shocked porous materials [43] and patterns in the phase separation of complex fluids [42]. Therefore, as the cavitation field is a typical multiphase flow system, the Min-MAM can be applied to quantify the profile evolution of cavitation bubble.

For a pixelized density map $\rho (\mathbf{r})$, the first step is to set a threshold value ${\rho }_{th}=0.25465 \mathrm{mu}/{\mathrm{lu}}^3$ and form a binary map, which can be expressed as

$${\rm P}(\mathbf{r},{\rho }_{th})=\{\begin{array}{c}0,\\ 1,\end{array}\begin{array}{c}\rho (\mathbf{r})<{\rho }_{th}\\ \rho (\mathbf{r})\ge {\rho }_{th}\end{array}$$

(14)

where $\mathbf{r}$ is the position, $\rho (\mathbf{r})$ is the density at $\mathbf{r}$, and ${\rm P}(\mathbf{r},{\rho }_{th})$ denotes the corresponding binary map. The grid node at $\mathbf{r}$ is regarded as a black vertex when $\rho (\mathbf{r})$ is less than ${\rho }_{th}$. Otherwise, it is regarded as white. A region containing connected black (white) pixels is defined as the black (white) domain [42]. The value of the black (white) domain is 0 (1), as shown in Figure 1. In the 2D space, the three Minkowski functionals geometrically correspond to the fractional area $F({\rho }_{th})$, boundary length $U({\rho }_{th})$, and Euler characteristic $X({\rho }_{th})$.

$F({\rho }_{th})$ represents the black domain in a binary map and can be given by

$$F({\rho }_{th})=\frac{N-N_{\rho th}}{N}=\frac{1}{N}\left[N-{\displaystyle \sum _{\mathbf{r}}{\rm P}\left(\mathbf{r},{\rho }_{th}\right)}\right]$$

(15)

where $N_{\rho th}$ is the number of white pixels and $N$ is the total number of pixels whose function is to normalize $N_{\rho th}$. The process of a cavitation bubble collapse near a solid wall is accompanied by the inward deformation of the bubble. Thus, $F({\rho }_{th})$ decreases from the maximum value to zero in a certain curve form.

The second morphological parameter of the Minkowski functionals is $U({\rho }_{th})$, which is expressed as

$$U({\rho }_{th})=\frac{B_{\rho th}}{N}=\frac{1}{N}{\displaystyle \sum _{{\rm P}(\mathbf{r},{\rho }_{th})=0}{\displaystyle \sum _{i=1}^4{\rm P}\left(\mathbf{r}+{\mathbf{e}}_i,{\rho }_{th}\right)}}$$

(16)

where $B_{\rho th}$ denotes the pixels that separate the black and white domains. $U({\rho }_{th})$ decreases to zero with $F({\rho }_{th})$. However, the two curves are significantly different because of the depression at the top of the cavitation bubble, causing the acceleration of $U({\rho }_{th})$ to decrease before the first collapse.

$X({\rho }_{th})$ can directly reflect the number of bubbles in the bubble collapse process and is given by

$$X({\rho }_{th})=N_{\rho th}^b$$

(17)

where $N_{\rho th}^b$ is the number of connected black domains.

When the first collapse occurs, $X({\rho }_{th})$ increases from one to two or more and decreases to zero after the second collapse. The vanishing $X({\rho }_{th})$ represents a complete collapse in the dynamic process of bubble near-wall collapse.

3. Cavitation Bubble Collapse near a Solid Wall with Complex Geometry

3.1. Numerical Model

The computational domain for the collapsing bubble near a solid wall with a complex geometry is established in a $512 \mathrm{lu}\times 512 \mathrm{lu}$ lattice system. The content inside the bubble is only vapor without other non-condensable gas. As shown in Figure 2, $p_v$ and $p_{\infty }$ are the internal and external pressure. $R_0$ denotes the maximum radius of the initial spherical bubble. $h$ represents the distance between the bubble center and the bottom of the solid wall. $\lambda =h/R_0$ is defined as the dimensionless position offset parameter. A complex geometry boundary with a bounce-back boundary condition [53] is applied in the bottom region. The periodic boundary is defined in the vertical direction, and a pressure boundary by the Zou–He scheme [54] is applied in the top region.

The process from the initial spherical bubble with the maximum radius to complete disappearance is studied to investigate the interaction mechanism between the bubble and solid wall. The maximum radius of the initial spherical bubble and the non-dimensional temperature are set as $R_0=80 \mathrm{lu}$ and $T/T_c=0.7$, respectively. The density field is initialized as [55]

$$\rho (x,y)=\frac{{\rho }_l+{\rho }_v}{2}+\frac{{\rho }_l-{\rho }_v}{2}\times \tan h\left[\frac{2\left(\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}-R_0\right)}{W}\right]$$

(18)

where ${\rho }_l$ and ${\rho }_v$ denote the densities of the liquid outside the bubble and the vapor inside the bubble, respectively, $(x_0,y_0)$ denotes the coordinate of the bubble center, $\tan h(x)=({\mathbf{e}}^{2x}-1)/({\mathbf{e}}^{2x}+1)$ is the hyperbolic tangent function, and $W=5 \mathrm{lu}$ represents the prescribed width of the phase interface.

In simulating the process of cavitation bubble collapse, we can artificially tune the initial density of the liquid to achieve a positive pressure difference $\Delta p=p_{\infty }-p_v$ based on the C–S EOS.

3.2. Solution of Boundary Conditions

The fractal dimension, used to quantitatively describe the surface complexity of the solid wall with an arbitrarily complex geometry, is expressed as [56]

$$D=\frac{\ln N\left(d\right)}{\ln \left(L/d\right)}$$

(19)

where $D$ is the fractal dimension which fluctuates between 1 and 2, $N(d)$ is the scale number which is the total number of steps of the curve measured with the scale of step $d$, and $L$ is the length of the straight line connecting the start and end points of the curve. $N(d)$ and $L$ are the lengths of the red polyline and green straight line in Figure 3, respectively.

In the previous works [29,31,32,37], the fluid–solid interactive force is zero, and the interaction between the fluid and solid wall is considered as fully rebounded. However, in contrast to a plane solid wall, a solid wall with arbitrarily complex geometry exhibits irregular geometric shapes and a solid wall orientation, causing the fluid–solid interface to undulate. Thus, an algorithm is proposed to automatically differentiate between the fluid and the solid wall. In the simulation, a label array $r({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)$ is adopted to signify the adjacent grid node that ${\mathbf{x}}_1$ reaches in a time step ${\delta }_t$ along a certain direction ${\mathbf{e}}_{\alpha }$, where ${\mathbf{x}}_1$ represents the position of a certain lattice node in the first layer fluid above the solid wall. Therefore, the value of $r({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)$ can be used to identify the fluid and solid wall. The grid node at location ${\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t$ is viewed as a solid wall node when the value of $r({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)$ is one. Otherwise, it is regarded as a fluid node.

In this study, the interaction between the fluid and solid wall is considered as fully rebounded, which can be expressed as

$$f_{-\alpha }({\mathbf{x}}_1,t+{\delta }_t)=f_{\alpha }({\mathbf{x}}_1,t) \mathrm{when} r\left({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t\right)=1$$

(20)

where $f_{\alpha }$ and $f_{-\alpha }$ denote the density distribution functions along the directions toward and away from the solid wall, respectively.

When the value of $r({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)$ is one, the interaction potential $\psi ({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)$ at location ${\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t$ is numerically the interaction potential at position ${\mathbf{x}}_1$. Therefore, the pseudopotential force ${\mathbf{F}}_{\alpha }$ can be expressed as

$$F_{\alpha }({\mathbf{x}}_1)=\{\begin{array}{c}-G\psi ({\mathbf{x}}_1){\displaystyle {\sum }_{\alpha }{\omega }_{\alpha }\psi ({\mathbf{x}}_1){\mathbf{e}}_{\alpha },}\hfill \\ -G\psi ({\mathbf{x}}_1){\displaystyle {\sum }_{\alpha }{\omega }_{\alpha }\psi ({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t){\mathbf{e}}_{\alpha },}\hfill \end{array}\begin{array}{c}r({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)=1\\ r({\mathbf{x}}_1+{\mathbf{e}}_{\alpha }{\delta }_t)=0\end{array}$$

(21)

The boundary condition and interaction potential of the interface between the liquid and rigid wall are discussed in this subsection. If the lattice node is fluid, it normally executes the collision and streaming processes, then enters the next round of evolution.

3.3. Verification of the Collapsing Bubble near the Solid Wall

In previous work, the validity of the pseudopotential MRT-LBM in simulating the cavitation bubble collapse process has been verified from the perspectives of coexistence curves and Laplace law verification [33]. As shown in Figure 4, the experimental and LBM simulation results are presented for comparison. Figure 4 depicts the profile evolutions of the collapsing bubble near a solid wall with complex geometry ($D=1.03476$). In the experiments, a metal block is attached to the solid wall sample to prevent the solid wall from moving, because the mass of the sample is light. So the bottom of the solid wall sample needs to be sketched in Figure 4a. It is clearly shown that the numerical result is in agreement with the experiment. These two groups of figures can present a similar dynamic process of the collapsing bubble. The validity of the numerical model to simulate the cavitation bubble near-wall collapse process is verified from the comparisons with the experimental photos. Owing to the solid wall with complex geometry, the bubble deformation is asymmetric. The evolutions of the physical fields will be analyzed in detail in the following subsections.

In addition, the grid independence of the numerical model is verified. As shown in Figure 5, three different lattice systems, i.e., $312 \mathrm{lu}\times 312 \mathrm{lu}$, $412 \mathrm{lu}\times 412 \mathrm{lu}$, and $512 \mathrm{lu}\times 512 \mathrm{lu}$, are chosen to carry out a more precise comparison and present asymmetry of the collapsing bubble near the solid wall with complex geometry. According to the same proportion principle, the maximum radii of the cavitation bubble are set as $R_0=50 \mathrm{lu}$, $R_0=65 \mathrm{lu}$, and $R_0=80 \mathrm{lu}$, respectively, and the offset parameters are both $\lambda =1.5$. Moreover, a non-dimensional time $T^{\ast }$ is defined to represent the interval between the initial state and the first collapse. Therefore, Figure 5 extracts the bubble profiles at 0, $0.2T^{\ast }$, $0.4T^{\ast }$, $0.6T^{\ast }$, $0.8T^{\ast }$, and $T^{\ast }$ with the acquired coordinate points normalized in the x and y directions.

The normalized profiles of the collapsing bubble near the solid wall with complex geometry almost coincide with one another under different grid resolutions. In other words, the collapse processes of the cavitation bubble near the solid wall under different grid resolutions are the same. Given the efficiency and accuracy of the numerical model, the computational domain is established in a $512 \mathrm{lu}\times 512 \mathrm{lu}$ lattice region in this manuscript.

3.4. Evolution of Density Field

Figure 6 shows the evolution of the density field, where the bubble profiles can be visually distinguished. This figure depicts the dynamic process of cavitation bubble collapse near two different solid walls, i.e., the plane solid wall and solid wall with a complex geometry ($D=1.03476$). The offset parameters are set as $\lambda =1.25$ and $\lambda =1.5$ in these two cases, respectively. The morphological evolution of the collapsing bubble in the vicinity of the solid wall with a complex geometry is analyzed in detail.

The initial spherical bubble starts to collapse under the pressure difference $\Delta p=0.0135994 \mathrm{mu}/\left({\mathrm{ts}}^2\mathrm{lu}\right)$. The collapse velocity of the bubble is low in its initial stage, and the bottom bubble wall exhibits a significant deformation at $t=100 \mathrm{ts}$. The bubble is elongated into an ellipsoid before $t=550 \mathrm{ts}$, owing to the retarding effect of the solid wall. With the higher density area diffusing, the upper bubble wall is flattened, then it becomes concave from $t=550 \mathrm{ts}$. There appears a denser density area in the concave region, which accelerates the change of the bubble profiles. Consequently, a micro-jet is formed because of the interaction between the bubble and the solid wall. When the upper and bottom bubble walls collide, the cavitation bubble splits into small bubbles. This phenomenon can be referred to as the first collapse. The bubble deformation is asymmetric, and the deformed bubble tends to collapse into two asymmetric small bubbles owing to the irregular geometry of the solid wall. Subsequently, the small bubbles continue to collapse, and the visible bubble surface completely disappears after the second collapse.

The processes of bubble collapse near the plane solid wall and solid wall with a complex geometry exhibit a similar dynamic process. However, the solid wall with a complex geometry causes the asymmetric deformation of the bubble. The effect of the solid wall with complex geometry on the collapsing bubble can be reflected in the profile evolution. Meantime, the validity and feasibility of the numerical model to simulate the bubble collapse process are again verified.

3.5. Evolution of Pressure and Velocity Fields

In contrast to the density field evolution, the evolution of the pressure and velocity fields can clearly and intuitively reflect the driving force behind the dynamic process of the collapsing bubble close to a solid wall with a complex boundary and reveal the corresponding physical mechanism. The simulated pressure and velocity fields are shown for several representative moments in Figure 7.

Initially, the cavitation bubble is spherical. Thereafter, it deforms into an elliptical bubble under the pressure difference $\Delta p=0.0135994 \mathrm{mu}/\left({\mathrm{ts}}^2\mathrm{lu}\right)$. A micro-jet is formed, and the micro-jet in the concave area is considerably stronger than that in other areas. The micro-jet velocity increases as the upper bubble wall sags. The micro-jet breaks through the bubble wall when the first collapse occurs at $t=705\mathrm{ts}$. The cavitation bubble splits into two asymmetric small bubbles owing to the retarding effect of the solid wall with a complex geometry. The smaller right bubble collapses first and forms a radial pressure propagation trend at the collapse point. After the second collapse at $t=770\mathrm{ts}$, the cavitation bubble completely collapses, and there appears a high-pressure shock wave near the surface of the solid wall. Then, the shock wave is driven by the micro-jet and spread to the surrounding solid wall until it disappears at the collapse point.

4. Morphological Analysis Method Exerted on the Cavitation Bubble Collapse Model

4.1. Morphological Verification

The phenomenon examined in this study belongs to multiphase flow. However, it is necessary to verify the validity and feasibility of the Min-MAM in quantitatively describing the evolution of collapsing bubble profiles. Therefore, Figure 8 depicts the time evolution of the Minkowski measures for the collapsing cavitation bubble near a solid wall with a complex geometry. This figure also includes the density field at several representative moments and the extended dot–line diagrams of the Minkowski measures for the period before the second collapse.

As shown in Figure 6, the collapsing bubble near the solid wall generally goes through the stages of inward collapse, top depression, micro-jet formation, first collapse, and second collapse. Therefore, as shown in Figure 8, the entire evolution process can be divided into the following parts: the first deformation from $t=0$ to $t=550\mathrm{ts}$, the second deformation before $t=705\mathrm{ts}$, the first collapse at $t=705\mathrm{ts}$, the circular bubble formation from $t=705\mathrm{ts}$ to $t=770\mathrm{ts}$, and the second collapse at $t=770\mathrm{ts}$. In the first deformation part, the initial spherical bubble deforms into an elliptical bubble under the pressure difference $\Delta p$. The fractional area $F$ and boundary length $U$ decrease during this period, and their slopes increase with time $t$. This indicates that the inward collapse velocity of the bubble increases gradually. In addition, the fluid velocity increases toward the solid wall, thereby forming an inceptive micro-jet. In the second deformation part, $F$ and $U$ continue to decrease, but their corresponding accelerations begin to increase. From $t=550\mathrm{ts}$, the bubble top wall begins to sag downwards, which partly compensates for the effect of the inward collapse on the bubble fractional area and boundary length. When the upper and bottom bubble walls collide at $t=705\mathrm{ts}$, the cavitation bubble splits into two small asymmetrical bubbles. After the first collapse, the visible bubble surface is composed of two asymmetrical bubbles in the 2D space and an annular bubble in the three-dimensional (3D) space. At this stage, the smaller bubble collapses earlier than the larger bubble, causing the slope of the $U(t)$ curve to increase significantly. After the second collapse, the values of $F$ and $U$ simultaneously decrease to zero.

As the third Minkowski measure, the Euler characteristic $X$ can directly reflect the number of bubbles in the bubble collapse process. $X$ remains at 1 before the first collapse. $X$ increases to 2 at $t=705\mathrm{ts}$ and remains at 2 before the smaller bubble collapses. $X$ decreases to 1 when the smaller bubble collapses. $X$ decreases to 0 after $t=770\mathrm{ts}$, indicating that the cavitation bubble completely collapses. In the curves for $U(t)$ and $X(t)$, the red circles in the figure can be described as “the first collapse point” and “the second collapse point”.

The results demonstrate that the three Minkowski measures can describe the morphological characteristics of the collapsing bubble in the vicinity of a solid wall with a complex geometry from different perspectives. The three curves can clearly distinguish between several typical moments in the cavitation bubble collapse process. Hence, we successfully verify the validity and feasibility of the Min-MAM in quantifying the morphological evolution of the collapsing bubble.

4.2. Evolution of Bubble Profiles

In the simulation, solid walls with arbitrarily complex geometries are synthesized, and their complexities are quantified by the fractal dimension. It is found in processing a fractal wall that two solid walls with different geometric features may have the same fractal dimension. Therefore, similar complex solid walls with different fractal dimensions and different complex solid walls with the same fractal dimension are designed as rigid boundaries to comprehensively analyze the effect of the fractal wall on the profile evolution of the collapsing bubble. Meanwhile, given that the bubble’s relative position to the solid wall has an important significance for the retarding effect of the solid wall, the influence of the position offset parameter on the bubble collapse process is also analyzed.

4.2.1. Complexity of Solid Walls

The first factor involves the complexity of the solid wall, which can be quantitatively described using the fractal dimension. As shown in Figure 9, there are five solid walls with different fractal dimensions but similar orientations. Figure 10 depicts the morphological metric curves of the Minkowski functionals for $\lambda =1.5$ and $D=1.03476$, $D=1.06289$, $D=1.08107$, $D=1.10391$, and $D=1.13124$.

The three Minkowski measures for the collapsing bubble near a solid wall with a complex geometry depict similar trends regardless of the fractal dimension. It can be obviously observed from Figure 9 that the surface of the solid wall with a higher fractal dimension is rougher than the one of the solid wall with a lower fractal dimension. The maximum scale difference of the solid wall surface is about $5\mathrm{lu}$ in these five solid walls, which is much less than the radius of the bubble $R_0=80\mathrm{lu}$. As a result, it takes more time to reflect the differences in the retarding effect of the solid wall on the collapsing bubble. Therefore, from the beginning to $t=550\mathrm{ts}$, the variation curves of the fractal dimension and boundary length almost coincide with each other. It is illustrated that the dynamic processes of the collapsing bubbles have very high similarity. In addition, there is a clear distinction between these curves from $t=550\mathrm{ts}$. This is because the differences in the retarding effect of the solid wall on the collapsing bubble begin to take effect as the time goes on.

There is no doubt that the retarding effect of the solid wall becomes stronger as the fractal dimension increases. The velocity of the bubble’s inward collapse will be different because of the variation in the retarding effect of the solid wall, which directly influences the size of the collapsing bubble at $t=550\mathrm{ts}$. In Figure 10b,c, as $D$ increases, the first and second collapses are delayed, and the interval between the first and second collapses evidently increases. These changes are also closely associated with the size of the collapsing bubble. The smaller the bubble, the faster it collapses. Therefore, the retarding effect of the solid wall on the collapsing bubble can be indicated by the occurrence time of the first and second collapses. In conclusion, the retarding effect of the solid wall is proportional to its complexity when the solid walls with complex geometries have similar orientations.

Figure 11 depicts the jet velocity evolution of the collapsing bubble from $t=550\mathrm{ts}$ to the first collapse. The figure also shows the corresponding boundary length of the collapsing bubble under five different fractal dimensions to clearly interpret the variations in the cavitation bubble jet velocity. The evolution of the jet velocity almost follows a trend, that is, it gradually increases to a peak and then slightly decreases with time. Moreover, the intervals at which the jet velocity reaches its peak are the same. The variation range of the jet velocity decreases as the fractal dimension increases. The trend in the boundary length indicates that an increase in the fractal dimension influences the depression speed of the collapsing bubble and causes the depression degree of the bubble top wall to increase slowly. When the depression degree of the bubble top wall increases with time, the surface tension increases correspondingly, causing the acceleration of the jet velocity to decrease. Therefore, the higher the fractal dimension, the slower the peak jet velocity. Based on the above analysis, an increase in the fractal dimension can also affect the evolution of the jet velocity.

4.2.2. Geometric Features of Solid Walls

As shown in Figure 12, three solid walls with the same fractal dimension are selected as rigid boundaries to further analyze the effect of the surface geometry of the solid wall on the collapsing bubble. From Figure 12, it can be found that the actual distance between the bubble and the solid wall’s surface has a great difference in these three solid walls. The difference between the solid walls a and b is up to $25\mathrm{lu}$. Figure 13 depicts the morphological metric curves of the Minkowski functionals for $\lambda =1.5$ and $D=1.03274$.

The three Minkowski measures for the collapsing bubble near a solid wall with a complex geometry no longer show similar trends, as shown in Figure 13. In the initial stage of bubble collapse, the actual distance from the bubble to solid wall’s surface plays a major role in the morphological evolution of the collapsing bubble. As a result, the curves of $F(t)$ and $U(t)$ have a certain difference before $t=300\mathrm{ts}$. Except for the actual distance between the bubble and the solid wall’s surface, the variations in the surface geometry of the solid wall also influence the dynamic process of the collapsing bubble. For instance, the first collapse near the solid wall b occurs before that near the solid wall c, but the corresponding second collapse occurs after that near the solid wall c. This is because the surface geometry of the solid wall influences its retarding effect on bubble collapse and leads to differences in the size of the collapsing bubble at the first collapse and the length of the collapse time. Therefore, the geometric feature of the solid wall plays a significant part in the evolution of the bubble profiles, thereby affecting the bubble collapse process.

Figure 14 depicts the evolution of the boundary length and jet velocity of the collapsing bubble near the three different solid walls from $t=550\mathrm{ts}$ to the first collapse. The trend of the jet velocity shows clear differences before the first collapse due to the differences in the surface geometry of the solid wall. For instance, the evolution of the jet velocity of the collapsing bubble near the solid wall a is similar to that in Figure 11. However, the jet velocities near solid walls b and c increase almost throughout the process. Moreover, the evolution of the corresponding boundary length of the collapsing bubble appears to have no apparent connection. Therefore, the surface geometry of the solid boundary affects the bubble collapse process, and there are no evident trends in the jet velocity evolution.

4.2.3. Dimensionless Position Offset Parameter

In this subsection, the effect of the position offset parameter on the collapsing bubble is discussed. Five offset parameters are selected, and the corresponding morphological metric curves of the Minkowski functionals are obtained under two different fractal dimensions. The position offset parameter influences the deformation speed of the collapsing bubble (Figure 15), causing the curves of the Minkowski functionals to clearly differ from $t=550\mathrm{ts}$. It is apparent that the deformation speed increases with $\lambda$. The retarding effect of the solid wall on the collapsing bubble is more suitable with the vapor bubble in the lower region of the liquid ($1.4<\lambda \le 1.5$), making the interval between the first and second collapses more obvious.

Figure 15I,II demonstrates that the absence of the second collapse is related to $\lambda$ and the geometric characteristics of the solid wall. The profile evolution of the collapsing bubble illustrates that the bubble deformation process is significantly affected by the solid wall. There is a significant difference between the shrinkage speeds of the left and right bubble walls because of the geometric characteristics of the solid wall, determining the size of the bubble at the first collapse. When the bubble at the first collapse is extremely small, the corresponding second collapse disappears during the entire process. Therefore, it is necessary to comprehensively consider various parameters that affect the deformation of bubbles to reduce the instability of the collapse process.

5. Conclusions

The interaction mechanism between a cavitation bubble and a solid wall with a complex geometry is investigated using the pseudopotential MRT-LBM. In this study, we introduce an arbitrarily complex geometry and propose an effective boundary processing method. A model of gas–liquid–solid three–phase coupling with adjustable solid boundary characteristics is established. The dynamics of the collapsing bubble are investigated using the evolution of the density, pressure, and velocity fields. The Min-MAM is used to quantitatively describe the evolution of the bubble profiles. Several representative moments of bubble collapse can be clearly distinguished from the morphological metric curves of the Minkowski functionals. Furthermore, the effect of different surface conditions of the solid walls on the dynamic characteristics of the collapsing bubble is systematically analyzed. The results indicate that the shift in the bubble profiles is inseparable from the surface geometry of the solid wall. The conclusions of this study are summarized below.

(1) The validity and feasibility of the numerical model to simulate the bubble collapse process are verified by comparing the LBM simulation with experimental results. The collapse process of a vapor bubble close to a solid wall with an arbitrarily complex geometry is accurately simulated. This provides theoretical support for investigating the dynamics of collapsing bubbles in fractal cracks.

(2) The Min-MAM is adopted to quantitatively describe the evolution of the bubble profiles and acquire effective information from the collapse process. The morphological metric curves of the Minkowski functionals can characterize several typical moments, including the top depression, first collapse, and second collapse of the bubble.

(3) The inherent properties of the solid wall, such as its complexity and surface geometry, should be considered when the solid wall is generalized from a plane solid wall to a solid wall with a complex geometry. Specifically, the influence of the solid wall on collapsing bubble profiles can be analyzed by tuning different parameters, e.g., the complexity of solid walls, geometric features of solid walls, and the position offset parameter. The results demonstrate that the bubble collapse becomes slower as the complexity of the solid wall increases. The geometric features of the solid wall can influence its retarding effect on the collapsing bubble, making the bubble collapse process more varied. The effect of the geometric features of the solid wall on the collapsing bubble is stronger than that of the surface complexity. In addition, the length of the collapse time is affected by the position offset parameter. Bubble collapse becomes faster as the offset parameter increases. However, the effect of the position offset parameter on the collapsing bubble is also associated with the geometric features of the solid wall, which influences the length of the collapse time and the occurrence of the second collapse.

Last but not least, based on the advantages of the pseudopotential MRT-LBM, the model in this paper can simulate the bubble dynamics near a solid wall with arbitrarily complex geometry. Several typical solid walls are introduced in this study to investigate the effect of fractal walls on the bubble collapse process, and general conclusions can be drawn. In the future, we will further study the temporal evolution of the collapsing bubble profiles near more solid walls with complex geometry or in a fractal crack. Moreover, the effect of the viscosity and surface tension on the bubble collapse model will be discussed in the subsequent study. Then we will attempt to extend the present 2D model into a 3D one and explore the dynamic characteristics of the collapsing bubble near a solid wall with complex geometry in the 3D space.