Study on the Collapse Process of Cavitation Bubbles Including Heat Transfer by Lattice Boltzmann Method

Abstract

In this study, an improved double distribution function based on the lattice Boltzmann method (LBM) is applied to simulate the evolution of non-isothermal cavitation. The density field and the velocity field are solved by pseudo-potential LBM with multiple relaxation time (MRT), while the temperature field is solved by thermal LBM-MRT. First, the proposed LBM model is verified by the Rayleigh–Plesset equation and D2 (the square of the droplet diameter) law for droplet evaporation. The results show that the simulation by the LBM model is identical to the corresponding analytical solution. Then, the proposed LBM model is applied to study the cavitation bubble growth and collapse in three typical boundaries, namely, an infinite domain, a straight wall and a convex wall. For the case of an infinite domain, the proposed model successfully reproduces the process from the expansion to compression of the cavitation bubble, and an obvious temperature gradient exists at the surface of the bubble. When the bubble collapses near a straight wall, there is no second collapse if the distance between the wall and the bubble is relatively long, and the temperature inside the bubble increases as the distance increases. When the bubble is close to the convex wall, the lower edge of the bubble evolves into a sharp corner during the shrinkage stage. Overall, the present study shows that this improved LBM model can accurately predict the cavitation bubble collapse including heat transfer. Moreover, the interaction between density and temperature fields is included in the LBM model for the first time.

1. Introduction

Cavitation is a normal phenomenon in hydraulic and marine engineering. For example, the cavitation bubble always exists around the water turbines, hydrofoils and ship propellers [1]. On one hand, the collapse of a cavitation bubble will generate extremely high pressure and temperature and this characteristic will cause great harm to the surrounding structures. On the other hand, the intense power of the cavitation bubble collapse can be utilized to increase the drilling rate of oil. So, it is very important to show the collapse process of the cavitation bubble and its interaction with the structures.

Cavitation refers to the complicated multiphase flow process, which includes formation, growth, shrinkage and collapse of the cavitation bubble. When the cavitation bubble collapses close to a solid wall, the bubble cannot maintain its spherical shape and a micro-jet is formed, causing considerable pressure. In this process, shock wave and the micro-jet are supposed as major reasons for mechanical damage. Kornfeld and Suvorov raised the idea of the micro-jet first and thought that the micro-jet is the major reason for the damage to the structure [2]. Observed by experimental methods, the cavitation bubble has been investigated for many years. There are three main ways of generating cavitation bubbles by experiments: spark [3,4], laser [5,6] and acoustic wave [7,8]. Naude observed the micro-jet during the experimental process [9]. Kling and Hammitt investigated the collapse of the cavitation bubble, and they also measured the harm caused by cavitation to the aluminum by experiments [10]. Vogel et al. measured the micro-jet and the counterjet velocity accurately and described the collapse of the cavitation bubble in different positions by high-speed camera technology [5,11]. Tomita and Kodama studied bubble collapse close to a composite surface [6]. Furthermore, the thermodynamic effects connected to the bubble collapse are also a vital factor. Dular and Coutier-Delgosha observed the temperature evolution during the growth and collapse stages of the bubble [12].

However, due to the complexity of cavitation collapse, it is extremely hard to obtain information on the whole flow field by experiments, especially for the final stage of collapse. Numerical simulation is another way to investigate the cavitation bubble. Traditionally, in order to simulate the multiphase flow, not only the discretization of the N-S (Navier-Stokes) equations are needed, but some methods that are applied to track the interface between different phases will also be coupled with the N-S equations. However, these numerical methods (such as the boundary element method, finite element method and the finite volume method) require interface capture [13], which adds additional difficulty to the simulation and decreases the efficiency of calculation. Furthermore, most of above methods calculate the pressure field by solving the Poission equation, which adds difficulty to the simulation.

Recently, the lattice Boltzmann method has been used by an increasing number of scholars to simulate the multiphase flows. The advantage of the LBM includes simple programming, good parallelism and easy implementation of the boundary condition. The collapse of the cavitation bubble can be simulated using the LBM pseudo-potential model proposed by Shan and Chen [14]. In the LBM pseudo-potential model, the pressure can be calculated by a non-ideal equation of state (EOS) [15]. Moreover, the interface between the phases would be captured automatically. The original LBM pseudo-potential model has been improved using the multi-relaxation time (MRT) [16] and an improved force scheme [17].

To date, the LBM has become a powerful method to simulate the cavitation bubble [18,19,20,21,22,23,24,25]. Sankaranarayanan et al. first investigated the bubble’s evolution with rise velocity by LBM [18]. Shan et al. and Mao et al. investigated the cavitation close to a solid wall by the pseudo-potential single-relaxation-time (SRT) model coupled with the exact differential method (EDM) force scheme [21,23]. Sukop and Or studied the homogeneous and heterogeneous cavitation by the LBM pseudo-potential BGK (Bhatnagar-Gross-Krook) model [19]. Shan et al. also investigated the homogeneous and heterogeneous cavitation by the pseudo-potential multi-relaxation-time (MRT) model coupled with an improved force scheme [21]. Peng et al. studied the evolution of a cluster of cavitation bubbles by the pseudo-potential model [24]. Liu and Peng studied cavitation bubble collapse close to a concave wall by the MRT-LBM model [25]. Su et al. verified the pseudo-potential SRT model by energy barrier theory [22] Mishra et al. studied the hydrodynamics of cavitation bubbles in the evolution of collapse coupled with chemical reactions [20]. It should be noted that these studies focus on the collapse mechanism of cavitation bubbles based on the isothermal model.

Because the collapse of the cavitation is instantaneous, it is assumed to be an adiabatic process in general. However, thermodynamic effects are also of great importance during the evolution, especially for the stage of collapse. So, it is necessary to investigate the temperature field in this process. To date, some simulations of bubble boiling and evaporation have verified that the LBM is also a mature method to investigate the heat transfer problem [26,27,28,29]. Generally, three strategies are utilized to couple the flow field and the temperature field with LBM: the Multi-Speed (MS) approach [30,31], Double-Distribution Function (DDF) approach [32,33,34,35] and the hybrid approach [36,37]. The first approach (MS) is an extension of the pseudo-potential LBM model, which is numerically unstable. The DDF approach describes the velocity field and the temperature field, respectively, by corresponding LBM models. For the hybrid approach, it solves the velocity field by LBM and temperature field by other methods. For the DDF method, the density field and the temperature field are solved based on the lattice Boltzmann equation. So, it is easy to implement the complicated boundaries compared to the hybrid approach.

Yang et al. studied the temperature field in the evolution of cavitation bubble collapse by the DDF approach without considering the feedback of the temperature field to the density field [38]. In the present study, the DDF approach improved by Li et al. [39] is used to simulate the evolution (non-isothermal) of the cavitation bubble in an infinite domain, near a straight wall and a convex wall. The density field and temperature field can be obtained by corresponding LBM models. Moreover, the interaction between the temperature field and the density field has been included in the simulation for the first time.

2. Material and Methods

In this section, a double distribution function (the density distribution function and the temperature distribution function (DDF)) based on LBM is introduced for simulation of the non-isothermal cavitation. The pseudo-potential LBM-MRT model was used for density and velocity fields, while thermal LBM-MRT model was used for the temperature field.

2.1. Pseudo-Potential LBM-MRT Model

For the density and velocity fields, a multiple relaxation time [14,40] was coupled with the LBM model (LBM-MRT), and a force scheme improved by Li et al. [17] was included into the model. The governing equation [17,41] is as follows:

$$f_{\alpha }\left(x+e_{\alpha }{\delta }_t,t+{\delta }_t\right)=f_{\alpha }(x,t)-{\overline{\Lambda }}_{\alpha }{}_{\beta }(f_{\alpha }(x,t)-f_{\alpha }{}^{eq}(x,t))+{\delta }_t(S_{\alpha }-0.5{\overline{\Lambda }}_{\alpha }{}_{\beta }{S}_{\beta })$$

(1)

where $f_{\alpha }$ denotes the local distribution function, $e_{\alpha }$ represents the lattice velocity, $\overline{\Lambda }=M^{-1}\Lambda M$, $f_{\alpha }$ represents the equilibrium distribution function and S denotes the source term.

M is defined as follows [42]:

$$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)$$

(2)

The diagonal matrix $\Lambda$ is defined as:

$$\Lambda =diag({\tau }_{\rho }{}^{-1},{\tau }_{\mathrm{e}}{}^{-1},{\tau }_{\zeta }{}^{-1},{\tau }_j{}^{-1},{\tau }_q{}^{-1},{\tau }_j{}^{-1}{\tau }_q{}^{-1},{\tau }_{\upsilon }{}^{-1},{\tau }_{\upsilon }{}^{-1})$$

(3)

By introducing the matrix M, the Equation (1) can be expressed as follows [43]:

$$m_{\alpha }^{*}=m_{\alpha }-{\Lambda }_{\alpha \beta }(m_{\beta }-m_{\beta }^{e\mathrm{q}})+{\delta }_t\left[I-\frac{{\Lambda }_{\alpha \beta }}{2}\right]{\overline{S}}_{\beta }$$

(4)

$I$ represents the unit matric; $m^{eq}$ is given by:

$$m^{e\mathrm{q}}=\rho (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)$$

(5)

$$\rho ={\displaystyle \sum _{\alpha }{f}_{\alpha }}$$

(6)

$$v={\displaystyle \sum \left(\frac{f_{\alpha }{e}_{\alpha }}{\rho }\right)}+\frac{{\delta }_{t}F}{2\rho }$$

(7)

where $\rho$ denotes the density and $v$ represents the macroscopic velocity. In addition, F denotes the total force, which includes the solid–fluid interaction, fluid–fluid interaction and the volume force. The fluid–fluid interaction can be expressed [44,45] as follows:

$$F_m=-G\varphi (x){\displaystyle \sum _{\alpha =1}^N\omega ({\left|e_{\alpha }\right|}^2)}\varphi (x+e_{\alpha })e_{\alpha }$$

(8)

where G denotes the strength of interaction, $\varphi (x)$ represents the interaction potential [46]. $\omega$ expresses the weight coefficient, In the D2Q9 model (For DdQq model, where d is the spatial dimensions and q is the number of velocities), $\omega (1)=1/3,\omega (2)=1/12$;$\varphi (x)$ can be expressed as follows:

$$\varphi (x)=\sqrt{2(p_{eos}-\rho c_s{}^2)/G{c}^2}$$

(9)

Similarly, the fluid–solid interaction can be expressed as follows [47]:

$$F_{ads}=-G_w\varphi (x){\displaystyle \sum _{\alpha =1}^N\varpi ({\left|e_{\alpha }\right|}^2)}S(x+e_{\alpha })e_{\alpha }$$

(10)

Here, $G_w$ denotes the interaction between the solid and the fluid. $\varpi ({\left|e_{\alpha }\right|}^2)=c_s{}^2\omega ({\left|e_{\alpha }\right|}^2)$, $c_s=c/\sqrt{3}$ represents the sound speed; c denotes the lattice velocity. $S(x+e_{\alpha })=\varphi (x)s(x+e_{\alpha })$, $s(x+e_{\alpha })$ is a switch function, which has different values for different phases.

In the present research, the pressure $p_{eos}$ was obtained by C-S EOS (Carnahan-Starling Equation of State) [48] as follows:

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

(11)

R represents the gas constant, T denotes the temperature and $a=0.4963R^2{T}_c{}^2/p_c,b=0.18727R{T}_c/p_c$. Because the product of $b$ and $\rho$ is the order of 0.01, the denominator cannot be zero. As for the negative value of pressure p, we know physically reasonable cases of flows when pressure is transformed to be negative. According to the results of the profound work of Sedov [49], in the vicinity of or inside the bubbles during the boiling in the fluid flows, pressure appears to be negative, but these are very special conditions for fluid flow. This is because such fluids are incapable of taking tensile loads (at negative pressure); such physical mechanism may be treated as a phenomenon of boiling of a liquid as soon as its pressure is lowered locally below zero.

$\overline{S}$ in the Equation (4) is defined [17] as follows ($\sigma$ affects thermodynamic consistency):

$$\overline{S}=\left(\begin{array}{c}0\\ 6(v_x{F}_x+v_y{F}_y)+\frac{12\sigma {\left|F_m\right|}^2}{{\varphi }^2{\delta }_t({\tau }_e-0.5)}\\ -6(v_x{F}_x+v_y{F}_y)-\frac{12\sigma {\left|F_m\right|}^2}{{\varphi }^2{\delta }_t({\tau }_{\zeta }-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)$$

(12)

where $F=F_m+F_{ads}$. The element in the diagonal matrix, including ${\tau }_e$ and ${\tau }_{\zeta }$, must be greater than 0.5 and less than 2. So, the denominator in $\overline{S}$ is a nonzero number.

The streaming process can be expressed as

$$f_{\alpha }\left(x+e_{\alpha }{\delta }_t,t+{\delta }_t\right)=f_{\alpha }{}^{\ast }(x,t)=M^{-1}{m}^{*}(x,t)$$

(13)

2.2. Thermal LBM-MRT Model

By introducing the temperature equation [34], the temperature can be obtained as follows:

$$\frac{\partial T}{\partial t}+\nabla \cdot (vT)=\nabla \cdot (\alpha \nabla T)+\Psi$$

(14)

where T is the temperature, $\alpha =k/(\rho c_v)$ denotes the thermal diffusivity, k represents the thermal conductivity and $c_v$ represents the specific heat. $\Psi$ is the source term that denotes the energy change due to the phase change, which can be written as [39]:

$$\Psi =\frac{1}{\rho c_v}\nabla \cdot (k\nabla T)-\nabla \cdot (\alpha \nabla T)+T\left[1-\frac{1}{\rho c_v}{\left(\frac{\partial p_{eos}}{\partial T}\right)}_{\rho }\right]\nabla \cdot v$$

(15)

In Equation (15), the subscript $\rho$ expressing the partial derivative of $p_{eos}$ with respect to T denotes the temperature and the density are independent variables of each other. The dissipation term in the velocity field of fluid flow was also taken into account [50]. The temperature equation (14), can be solved by the thermal LBM-MRT model, which is proposed by Li et al. [39], given by:

$$g_{\alpha }\left(x+e_{\alpha }{\delta }_t,t+{\delta }_t\right)=g_{\alpha }(x,t)-{\overline{\Omega }}_{\alpha }{}_{\beta }(g_{\alpha }(x,t)-g_{\alpha }{}^{eq}(x,t))+{\delta }_t{Q}_{\alpha }(x,t)$$

(16)

where $g_{\alpha }$ represents the temperature distribution function, $Q_{\alpha }$ is the source term and $\overline{\Omega }$ denotes the relaxation matrix. Like the density field, the above equation could be simplified as follows:

$$n_{\alpha }^{*}=n_{\alpha }-{\Omega }_{\alpha \beta }(n_{\beta }-n_{\beta }^{e\mathrm{q}})+{\delta }_t{Q}_{\alpha }^{’}$$

(17)

where $n_{\alpha }=M_{\alpha \beta }{g}_{\beta }$, $n_{{}_{\alpha }}^{eq}=M_{\alpha \beta }{g}_{{}_{\beta }}^{eq}$, $\Omega$ is the diagonal matrix and $Q_{\alpha }^{’}$ represents the source term. $\Omega$, $n^{eq}$ and $Q_{\alpha }^{’}$ are given by [51]:

$$\Omega =\left(\begin{array}{ccccccccc}{\tau }_0{}^{-1}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& {\tau }_1{}^{-1}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& {\tau }_2{}^{-1}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& {\tau }_3{}^{-1}& \left(\frac{{\tau }_3{}^{-1}}{2}-1\right){\tau }_4{}^{-1}& 0& 0& 0& 0\\ 0& 0& 0& 0& {\tau }_4{}^{-1}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& {\tau }_5{}^{-1}& \left(\frac{{\tau }_5{}^{-1}}{2}-1\right){\tau }_6{}^{-1}& 0& 0\\ 0& 0& 0& 0& 0& 0& {\tau }_6{}^{-1}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& {\tau }_7{}^{-1}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& {\tau }_8{}^{-1}\end{array}\right)$$

(18)

$$n^{e\mathrm{q}}=T(1,-2,2,v_x,-v_x,v_y,-v_y,0,0)$$

(19)

$$Q_{\alpha }^{’}=(Q_0,0,0,0,0,0,0,0,0)$$

(20)

where $Q_0=\Psi +0.5{\delta }_t{\partial }_t\Psi$. In the improved thermal LBM-MRT model, ${\partial }_t\Psi =\left[\Psi \right(t)–\Psi (t–{\delta }_t\left)\right]/{\delta }_t$.

The streaming process for the temperature distribution function is the same as for the density field:

$$g_{\alpha }\left(x+e_{\alpha }{\delta }_t,t+{\delta }_t\right)=g_{\alpha }{}^{\ast }(x,t)$$

(21)

The temperature can be calculated by:

$$T={\displaystyle \sum _{\alpha }{g}_{\alpha }}$$

(22)

It is clear that the macroscopic velocity and density calculated by the pseudo-potential LBM model influences the thermal LBM model by the source term in Equation (15). The temperature obtained by the thermal LBM model can also influence the pseudo-potential LBM model by EOS.

3. Validation of the DDF-LBM Model

3.1. Comparisons of Simulation with the Analytical Solution

First, the collapse of a bubble in an infinite domain was simulated to verify the used DDF-LBM model. At the beginning, due to the pressure difference between the inside and outside the bubble, the bubble shrinks and then disappears.

In the present simulation, the computational domain is 801 × 801 grids.

Table 1. The variable units in the study.

Variable	Unit
mass	$\mathrm{mu}$
length	$\mathrm{lu}$
time	$\mathrm{tu}$
density	$\mathrm{mu}\cdot \mathrm{l}{\mathrm{u}}^{-3}$
pressure	$\mathrm{mu}\cdot \mathrm{l}{\mathrm{u}}^{-1}\mathrm{t}{\mathrm{u}}^{-2}$
speed	$\mathrm{lu}\cdot \mathrm{t}{\mathrm{u}}^{-1}$

shows the units in the present study:

The temperature unit is dimensionless, and the lattice speed is $c=\mathrm{l}\mathrm{u}/\mathrm{t}\mathrm{u}$. The properties of fluid at the interface can be obtained as follows [52]:

$$\chi ={\chi }_l\cdot \frac{\rho -{\rho }_v}{{\rho }_l-{\rho }_v}+{\chi }_v\cdot \frac{{\rho }_l-\rho }{{\rho }_l-{\rho }_v}$$

(23)

where $\chi$ represents $c_v$ or k. The periodic boundary is applied to the four sides. The bubble locates at the center of the domain and shrinks under the pressure difference.

For this case, the modified Rayleigh–Plesset equation [53] can be used to obtain the analytical solutions:

$$\ln \left(\frac{r_{\infty }}{R}\right)\cdot \left({\stackrel{·}{R}}^2+\stackrel{··}{R}\right)-\frac{{\stackrel{·}{R}}^2}{2}+\frac{2{\mu }_l}{{\rho }_{l}R}\stackrel{·}{R}-\frac{\sigma }{{\rho }_{l}R}=\frac{p_v-p_{\infty }}{{\rho }_l}$$

(24)

where $r_{\infty }$ is the boundary position, $R$ denotes the radius of the bubble, $\stackrel{·}{R}=dR/dt$ and $\stackrel{··}{R}={d^{2}R/dt}^2$. ${\mu }_l=c_s{}^2\left(2{\tau }_v-1\right)/2$ is the viscosity of fluid. $\sigma$ is the surface tension and can be calculated by Laplace’s law. Although the model is non-isothermal, $\sigma$ can be obtained by the isothermal model, which does not include the temperature function. The initial bubble radius $R_0=40\mathrm{lu}$. The initial temperature $T_0=0.5T_c$, ${\mu }_l=0.1$, $\sigma =0.0105\mathrm{m}\mathrm{u}\cdot \mathrm{t}{\mathrm{u}}^{-2}$ and the pressure difference $p_{\infty }–p_v=0.002647848 \mathrm{m}\mathrm{u}\cdot \mathrm{l}{\mathrm{u}}^{-1}\mathrm{t}{\mathrm{u}}^{-2}$. Because the Equation (24) is a second-order equation, another initial condition is required: At $t=0$, $\stackrel{·}{R}=0$. In all simulations, the diagonal matrix ${\tau }_0={\tau }_3={\tau }_4={\tau }_5={\tau }_6={\tau }_8=1.0,{\tau }_1={\tau }_2=1.1$ ${\tau }_{\rho }={\tau }_e={\tau }_j={\tau }_q=1.0,{\tau }_{\zeta }={\tau }_j={\tau }_q=1.1,{\tau }_{\upsilon }=0.8$. The parameters in the EOS are set as follows: a = 0.5, b = 1. These parameters are constant in this paper. Because the temperature in the interface between phases is nearly constant, Equation (24) can be solved conveniently by the fourth order Runge–Kutta method.

Figure 1 shows the comparison of the radius evolution between the analytical solution and the results by the DDF-LBM model. It is clear that the collapse processes by the analytical method and DDF-LBM model are basically identical. The bubble radius is equal to 0 at $t_{\max }$. The results show that the DDF-LBM model can predict cavitation bubble collapse reasonably.

3.2. D2 Law for Droplet Evaporation

Second, the DDF-LBM model was validated further by the D2 (the square of the droplet diameter) law for droplet evaporation in an infinite domain, which indicated that the relationship between the square of the droplet diameter and the time is linear [54]. The law is based on the following assumptions: the liquid and the vapor phases are both quasi-steady; the viscous heat dissipation and buoyancy are neglected. Furthermore, the thermophysical properties are constant, such as k and $c_v$. In this case, we located the droplet at the center of a domain. The density field can be described by:

$$\rho (x,y)=\frac{({\rho }_l+{\rho }_g)}{2}–\frac{({\rho }_l-{\rho }_g)}{2}\tanh \left[\frac{2\left(\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}-R_0\right)}{W}\right]$$

(25)

where g and l represent the vapor and the liquid, respectively. W represents the width of the interface.

Similarly, the temperature field can be initialized in the same way:

$$T(x,y)=\frac{(T_{\infty }+T_b)}{2}+\frac{(T_{\infty }-T_b)}{2}\tanh \left[\frac{2\left(\sqrt{{\left(x-x_0\right)}^2+{\left(y-y_0\right)}^2}-R_0\right)}{W}\right]$$

(26)

where $T_{\infty }=0.5T_c$ and $T_b=T_c$ are the temperature outside and inside of the bubble in this section, respectively.

In the numerical simulation, the computational domain is 401 × 401 grids. $R_0=50\mathrm{l}\mathrm{u}$. Due to the temperature difference at the interface between the phases, the droplet evaporation happens. $c_v=10, k=0.1$.

Figure 2 shows that the results of the DDF-LBM model are identical to those of the D2 law. However, at the initial stage, there are some errors, which may be because the initial density field is not quasi-steady. Overall, the DDF-LBM model is reliable for the simulation of cavitation.

4. Results and Discussions

After the above verification, three typical cases of cavitation (in an infinite domain, close to a solid wall and a convex wall) were simulated by the improved DDF-LBM model. In the following, the changes of velocity and temperature information during the evolution of the cavitation bubble are discussed in detail.

4.1. Evolution of the Cavitation in an Infinite Domain

In this section, the evolution of cavitation in an infinite domain is simulated, which includes the growth, the shrinkage and the depression of the cavitation bubble. A $501\times 501 {\mathrm{lu}}^2$ square domain was used for the simulation. The periodic boundary is adopted to the four sides. The initial bubble radius $R_0=30 \mathrm{lu}$, the density field and the temperature field are initialized by Equations (25) and (26), where ${\rho }_l,{\rho }_g$ is 0.42 and 0.0006 $\mathrm{m}\mathrm{u}\cdot \mathrm{l}{\mathrm{u}}^{-3}$, leading to a pressure difference of $p_{\infty }–p_v=0.011 \mathrm{m}\mathrm{u}\cdot \mathrm{l}{\mathrm{u}}^{-1}\mathrm{t}{\mathrm{u}}^{-2}$. $T_{\infty }=0.5T_c$ and $T_b=T_c$.

Figure 3 shows the experimental results [55] and simulated results by the DDF-LBM model. It can be found that the simulation by the DDF-LBM model agrees with the experimental results. At first, the bubble expands due to the pressure difference between the inside and outside of the bubble. When the bubble reaches the maximum volume, because the additional pressure at the boundaries transfers to the liquid near the bubble, the ambient pressure is stronger than that in the bubble. So, it starts to shrink and collapse. The results illustrate that there is an obvious different temperature at the surface of the bubble. The temperature field is displayed in Figure 3c. It is clear that the temperature inside the bubble decreases in the process of expansion. Then, it keeps increasing during the shrinkage process and reaches the peak when it collapses.

4.2. Evolution of the Cavitation Near a Straight Wall

In this section, the evolution of cavitation near a straight wall is discussed in detail. Different from the evolution of the cavitation in an infinite area, the bubbles in this case do not remain spherical when it collapses close to a straight wall. The domain is also set as 501 × 501 $\mathrm{l}{\mathrm{u}}^2$. The periodic boundary is adopted for the right and left sides. In addition, the open boundary [56,57,58] is implemented in the top of the domain. A bounce-back [58,59] boundary is applied to the bottom side. The initial bubble radius $R_0=70 \mathrm{l}\mathrm{u}$, the initial density field and the initial temperature field are obtained by the same method; ${\rho }_l,{\rho }_g$ is 0.485 and 0.0006 $\mathrm{m}\mathrm{u}\cdot \mathrm{l}{\mathrm{u}}^{-3}$, $T_{\infty }=T_b=0.5T_c$. As shown in Figure 4, another parameter $\lambda =h/R_0$ is proposed, which represents the distance between the wall and the bubble center.

Figure 5 shows the evolution of the cavitation by the experimental method [60] and the simulation results by the DDF-LBM model when λ = 1.5. The simulation results are identical to the experimental results. At the beginning, due to the pressure difference between the inside and outside of the bubble, the bubble shrinks. The bottom wall of the bubble moves slowly due to the blockage of the wall. Compared to the contraction in the longitudinal direction, the lateral contraction is more obvious, which leads to the bubble shaping into an ellipse. The depression gradually deepens at the upper bubble edge as the pressure difference increases. When the lower bubble edge collides with the upper bubble edge, the collapse occurs.

Figure 6 shows the evolution of the cavitation in different positions close to a straight wall. When λ is relatively small (λ = 1.2, 1.8), the evolution process of the bubble is similar to Figure 5. First, it gradually shrinks under the effect of the ambient pressure, and then it collapses and disappears under the pressure generated by the upper boundary. When λ = 2.4, the collapse of the bubble is not consistent with the other two cases. The cavitation bubble shrinks in a spherical shape under the influence of ambient pressure, instead of shrinking into an elliptical shape. The volume of the cavitation bubble is already very small when it begins to depress, and the bubble collapses without a full development of depression. It can be shown from these cases that the bottom edge speed of the bubble gradually increases as λ increases, resulting in a gradual increase in the shrinkage of the bubble, but the occurrence time of the depression is almost the same, so the degree of depression gradually decreases; the cavitation bubble in Figure 6c did not collapse twice.

Figure 7 shows the velocity and the temperature fields during the cavitation bubble collapse in different positions. From the previous analysis, it is known that the effect of the wall on bubble reduces with the increase in distance between the bubble and wall, which mainly affects whether the cavitation bubble has secondary collapse or not. It can be seen that the velocity at the bottom of the bubble is nearly zero in Figure 7a because λ is small. As λ increases, the velocity of the cavitation bubble edge also increases, and the maximum velocity appears on the upper edge of the cavitation bubble. As the cavitation bubble continues to shrink, the peak of the temperature inside the bubble also appears in the upper half. When the bubble starts to depress, the velocity of the upper bubble edge increases rapidly. Therefore, the process from cavitation bubble depression to collapse has a short duration, and the location of the highest temperature in the bubble also moves downward. The first collapse produces a pressure wave. With the combined effects of this pressure wave and the micro-jet, the remaining annular bubble eventually collapses. During this process, the location of the highest temperature also moves from the center to the inside of the annular bubble, which is the second collapse.

The bubble in Figure 7b is far from the solid wall. The shape of the bubble during the collapse is different from that in Figure 7a. First, due to the small obstruction by the wall, the lower edge of the bubble has a certain velocity, but it is lower than the velocity of the upper bubble edge. The cavitation bubble shrinks approximately in a spherical shape, and there is a low velocity area in the bubble. As the pressure above the cavitation bubble gradually rises, the velocity of the upper bubble edge gradually increases, and the temperature of the upper portion of the bubble gradually increases, so that a depression appears. Since the volume of the cavitation bubble is already small when it is depressed, it collapses before the depression is fully developed, and there is no second collapse process. In the process of depression, the location of the temperature peak in the bubble gradually moves downward. The reason why there is no second collapse in this case is that λ is large enough, and the lower bubble edge has a certain velocity in the collapse process, so that the micro-jet cannot be fully developed. So, the cavitation bubble disappears during depression.

Figure 8 displays the max temperature evolution in different cases. If λ ≥ 2.0, the peak of the temperature increases as the λ increases. It can be found from the Figure 8 that the temperature peak in the cavitation bubble keeps rising during the shrinkage stage. In the depression stage, the temperature changes of various cases are not consistent. When λ < 1.8, the temperature peak of the cavitation bubble in the initial stage of the depression decreases slightly. It can be seen from Figure 7a that the position of the highest temperature in the bubble at this stage gradually moves down, with a small increase. When the position moves to the lower bubble wall, it stops moving downward. During the later period of the depression, the micro-jet has been fully developed, and the compression of the bubble causes the rise in temperature in the bubble. Then, it does not decrease until the second collapse of the cavitation bubble occurs. When λ ≥ 1.8, the maximum temperature in the cavitation bubble shows a trend of increase in the shrinkage and depression stages. Since there is no second collapse, the temperature peak position moves to the lower part.

4.3. Evolution of the Cavitation Near a Convex Wall

In this section, a single cavitation bubble’s growth and collapse near a convex wall is described in detail. In practical engineering, the cavitation often happens near an irregular boundary, so it is vital to investigate the cavitation in complicated boundary conditions. There are also some experiments on cavitation bubble collapse near a convex wall [61], which can be used to verify the simulation. The computational layout for this case is shown in Figure 9. For the density and velocity fields, the open boundary is implemented for the top, right and left boundaries. Additionally, a half-way bounce-back boundary is adopted to the bottom of the domain. As for the temperature field, an isothermal boundary is applied to these four sides. The initial bubble radius $R_0$ is 40 $\mathrm{l}\mathrm{u}$. ${\rho }_l,{\rho }_g$ is 0.422 and 0.0006 $\mathrm{m}\mathrm{u}\cdot \mathrm{l}{\mathrm{u}}^{-3}$, $T_{\infty }=0.5T_c,T_b=T_c$. The computational domain is 501 × 501 grids. The curvature of solid wall is controlled by a function c (x,y) [61]:

$$c(x,y)=\sqrt{x^2+{\left(h-y\right)}^2}\left({\xi }^2\left(y+\xi h\right)+\left(1-{\xi }^2\right)\sqrt{\left(x^2+{\left(y+\xi h\right)}^2\right)}\right)-\left(h-y\right)\sqrt{\left(x^2+{\left(y+\xi h\right)}^2\right)}$$

(27)

where h is the height of the wall, l represents the distance from the bubble center to the wall. $\xi$ is in inverse proportion to the curvature of the wall, which indicates the curvature increase as $\xi$ decreases. The curvature at point M can be defined as [61]: $3(\xi -1)/(4\xi )$.

Figure 10 illustrates the evolution of the bubble by experiments [61] and simulation by the DDF-LBM model. It can be shown that the numerical results are identical to those of the experiments. At the beginning, the external pressure of the bubble is lower than the internal pressure. Under this pressure difference, the cavitation bubble expands and grows. After the growth reaches a certain stage, the lower part of the cavitation bubble develops slowly due to the inhibition of the convex wall. When the boundary pressure spreads around the bubble, it begins to shrink. In the shrinking process, compared with the flat wall, the lower bubble edge of the cavitation bubble does not remain smooth, but a sharp corner is evolved. The upper half of the bubble is similar to the shrinking process for the case with a flat wall.

Figure 11 shows the whole process of a bubble’s growth and collapse close to a convex wall with two curvatures. It can be seen that the bubble is basically not affected by the curvature of the wall during the growth stage. The cavitation bubble has a certain outward movement in the initial stage of growth, and the temperature inside the bubble gradually decreases. However, because the lower part of the bubble is suppressed by the wall, it would stop developing when close to the wall. When the boundary pressure spreads around the cavitation bubble, it starts to shrink. Moreover, due to the change of the curvature, the evolution of the velocity of the bubble edge is also different. When the wall curvature is large, in other words, $\xi$ is relatively small, the area occupied by the wall is smaller, resulting in a sharper shape of the lower part of the bubble. However, the velocity of the lower part of the bubble is still lower than that of the upper part. In the initial stage of shrinkage, the location of maximum velocity is always close to the center of the upper bubble wall. In the late stage of shrinkage, it can be found that the upper bubble edge also forms a sharp angle, and the maximum velocity appears on both sides of the upper bubble edge. The bubble temperature continues to rise during the shrinking process. If the curvature is small, the area affected by the wall becomes larger, and the velocity of the lower part of the cavitation bubble becomes smaller, so there is no obvious sharp corner. During the depression of the cavitation bubble, the top edge speed increases sharply, the temperature inside the bubble also increases and the position of peak temperature moves downward accordingly. It can be seen that if the curvature is small, the second collapse occurs. However, if the curvature is large, the lower sharp corner is obvious, and the second collapse does not exist.

Figure 12 shows the evolution of the maximum temperature inside the bubble and the change of pressure in the center of the convex wall during the growth and collapse processes of the cavitation bubble under different curvatures. During the growth stage of the cavitation bubble, the maximum temperature in the bubble shows some fluctuations to a certain extent, but it is close to the ambient temperature. When the bubble starts to shrink, the highest temperature inside the bubble rises sharply, and reaches the extremum when it collapses. Among these three cases, the largest temperature appears in the $\xi =0.22$ case. In this case, compared to the lower part, the bubble has a larger velocity on the upper part of the bubble. The compression speed of the bubble is the fastest. The pressure at the center of the convex wall is close to 0 before the collapse of the bubble. When the pressure wave induced by collapse reaches the wall, the pressure increases sharply. The maximum pressure on the convex wall comes from the $\xi =0.75$ case. This is because the cavitation bubble has a second collapse, and the time between the two collapses is short. The pressure waves generated by the two collapses superimpose and impact the convex wall.

5. Conclusions

In this study, an improved DDF-LBM model was applied for the simulation of the non-isothermal cavitation. The density and velocity fields were simulated by the pseudo-potential model, and the temperature field was simulated by the thermal LBM model. After the proposed LBM model was verified, it was used to study the cavitation bubble collapse in three typical boundaries, namely, an infinite domain, a straight wall and a convex wall. In view of the simulation, the following conclusion can be obtained:

(1)

The LBM model successfully reproduces the process from the expansion to compression of the cavitation bubble for the case of an infinite domain. In the initial growth stage, the temperature inside the bubble decreases, and the high pressure and the temperature are found at the final collapse stage due to the bubble compression.

(2)

When a cavitation bubble is near the straight wall, λ determines whether the second collapse exists or not. During the shrinkage stage, the temperature inside the bubble keeps rising, and the position of the maximum temperature in the bubble moves downwards in the depression stage. If the second collapse exists, the temperature inside the bubble has a slight fluctuation during the depression stage and the position of the maximum temperature has lateral movement following the annular bubble after the first collapse. When λ is large enough (λ = 2.2, 2.4), the temperature inside the bubble is obviously higher compared with the other cases.

(3)

When a bubble is near a convex wall, the curvature of the wall can affect the behavior of the bubble collapse. In the shrinkage stage, the lower edge of the bubble evolves into a sharp corner due to the compression of both sides of the bubble’s lower part for the case of $\xi =0.22$. However, the sharp corner does not exist when $\xi =0.75$, and the bubble has a second collapse, which is similar to the bubble near the straight wall. However, the temperature for the case of $\xi =0.22$ in the evolution is the largest. When $\xi =0.75$, the pressure at the center of the convex wall is the largest, which indicates that the wall pressure reduces as the $\xi$ decreases.

(4)

Overall, the improved LBM model can accurately predict cavitation bubble collapse, including the heat transfer. Moreover, the interaction between density and temperature fields has been included in the LBM model for the first time.