f α ( x + e α δ t , t + δ t ) = f α ∗ ( x , t ) , (7) where f ∗ = M − 1 m ∗ . The macroscopic velocity of the Eq. 5 can be calculated by ρ v = ∑ α e α f α + δ t 2 F , (8) where the fluid-fluid interactive force is F = ( F x , F y ) in 2D space, which is given by (Shan, 2006; Shan, 2008) F ( x ) = − G ψ ( x ) ∑ α = 1 N ω ( | e α | 2 ) ψ ( x + e α ) e α , (9) where ω ( | e → α | 2 ) denotes the weight coefficient with ω ( 1 ) = 1/3 and ω ( 2 ) = 1/12, G represents the strength of interaction between fluids. ψ ( x ) is the interaction potential, which is given by (Sbragaglia et al., 2007) ψ ( x ) = 2 ( p E O S − ρ c s 2 ) / G c 2 , (10) where c S = c / 3 is the lattice sound speed, p E O S is a prescribed non-ideal equation of state (EOS), coupling the Peng-Robinson (P-R) EOS, which is given by (Yuan and Schaefer, 2006) p E O S = ρ δ T 1 − b ρ − a ϑ ( T ) ρ 2 1 + 2 b ρ − b 2 ρ 2 , (11) where a = 0.45724 R 2 T c 2 / p c , b = 0.0778 R T c / p c , and δ is the gas constant. ϑ ( T ) = [ 1 + ( 0.37464 + 1.54226 ω − 0.26992 ω 2 ) ( 1 − T T c ) ] 2 , in which ω = 0.344 is the acentric factor, T _c and p _c indicate the critical temperature and pressure. The constants a = 2/49, b = 2/21 and δ = 1 are used in the present study (He et al., 1998).

For the D2Q9 lattice, Li et al.’s (Li et al., 2012; Li et al., 2013) forcing scheme can be given by S = [ 0 6 ( v x F x + v y F y ) + 0.75 ε | F | 2 ψ 2 δ t ( τ e − 0.5 ) − 6 ( v x F x + v y F y ) − 0.75 ε | F | 2 ψ 2 δ t ( τ ς − 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 ) ] . (12)

The adjustable coefficient ε , seen in Eq. 12, which can make the mechanical stability solution and the thermodynamic consistency requirement approximately the same in a wide temperature range.

Now turns to solving the temperature equation, the target temperature equation is written as (Guo et al., 2002; Li et al., 2016) ∂ T ∂ t + ∇ ⋅ ( v T ) = ∇ ⋅ ( α ∇ T ) + ϕ , (13) where the thermal diffusivity is α = k ρ c v , k indicates the thermal conductivity, ρ denotes the macroscopic density and c v is the specific heat at constant volume. The source term ϕ can be written as ϕ = T [ 1 − 1 ρ c v ( ∂ p E O S ∂ T ) ρ ] ∇ ⋅ v . (14)

The discrete equation of temperature distribution function can be given by g α ( x + e α δ t , t + δ t ) − g α ( x , t ) = − Λ ¯ α β ( g β − g β e q ) | ( x . t ) + δ t Q α ' ( x , t ) , (15) where g α ( x , t ) is the temperature distribution function at position x and time t, Q α ' represents the source term of the discrete velocity space, and Λ ¯ α β = ( M − 1 Λ M ) α β is the collision matrix.

Eq. 15 can be converted as follows n ∗ = n − Λ ( n − n e q ) + δ t Q , (16) where n = M g with n = ( n 0 ∗ , n 1 ∗ , ... , n 8 ∗ ) Τ and n e q = M g e q , in which g = ( g 0 , g 1 , ... , g 8 ) Τ and the equilibrium distribution g e q = ( g 0 e q , g 1 e q , ⋯ , g 8 e q ) Τ . In the D2Q9 model, n e q are given by n e q = T ( 1 , − 2 , 2 , v x ′ − v x ′ v y , − v y ′ 0 , 0 ) T . (17)

Q is the source term in the momentum space, which is shown as Q = ( Q 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ) T , (18) where Q 0 = ϕ + 0.5 δ t ∂ t ϕ . According Ref. (Chopard et al., 2009), in the numerical implementation, Eq. 18 is calculated with ∂ t ϕ ≈ [ ϕ ( t ) − ϕ ( t − δ t ) ] δ t .

The streaming process of temperature function is similar to Eq. 7 g α ( x + e α δ t , t + δ t ) = g α ∗ ( x , t ) , (19) where g ∗ = M − 1 m ∗ . It should be emphasized that the unit used in this work is the LBM lattice unit. The length’s unit is lu, time step’s unit is ts and the unit of mass is mu. The units of velocity, density, pressure, and viscosity are given as lu/ts, mu/lu ³ , mu/(ts ² lu) and lu ² /ts, respectively.

Laplace’s law states that in the case of cavitation bubbles, the pressure difference between inside and outside the bubble is inversely proportional to the bubble radius, which is shown as Δ p = p i n − p o u t = σ R , (20) where p i n indicates the pressure inside the bubble, p o u t indicates the pressure outside the bubble, σ represents the surface tension, and R is the equilibrium radius of the bubble. A spherical bubble in the liquid domain is located in the center of a periodic field without gravity. The initial density field formula is (Shan et al., 2016) ρ ( x , y ) = ρ i n + ρ o u t 2 + ρ i n − ρ o u t 2 × t a n ⁡ h [ 2 ( ( x − x 0 ) 2 + ( y − y 0 ) 2 − R ) W ] , (21) where tan ⁡ h ( x ) = ( e 2 x − 1 ) / ( e 2 x + 1 ) is the hyperbolic tangent function. ρ i n represents the density inside the bubble and ρ o u t represents the density outside the bubble. The width of the phase interface W = 5. ( x 0 , y 0 ) is the central coordinates of the 401 × 401 lattice system. At dimension temperature 0.85T _c, the initial vapor and liquid density are approximately equal to the 0.34 and 6.6. The initial bubble radius R ₀ is set as 20, 30, 40, 50 and 60. It can be clearly observed from Figure 1 that the direct proportionality of the pressure difference Δp to the inverse radius 1/R, indicating that the simulation results are in good agreement with Laplace’s law.

As is well known, the bubble dynamics equation in a viscous fluid—i.e., the R–P equation—is (Brennen, 1973) R R ¨ + 3 2 R ˙ 2 = 1 ρ l [ ( p ∞ + 2 σ R − p v ) ( R 0 R ) 3 γ + p v − p ∞ − 4 μ R ˙ R − 2 σ R ] , (22) where p v is the vapor pressure in the bubble, p ∞ is the liquid pressure. Dots represent time derivative: R ˙ = d R / d t , R ¨ = d 2 R / d t 2 . γ is the adiabatic coefficient. μ = ( τ v − 1 2 ) c s 2 δ t is dynamic viscosity. In order to compare this simulation result with the numerical solution and R-P equation solution, the corresponding LBM parameters are substituted into Eq. 22 to obtain the solution of the R-P equation. The pressure difference Δ p = p ∞ − p v =   0.0853 obtained is used to trigger the collapse of the bubble. The initial densities of vapor and liquid are ρ v ≈ 0.34 and ρ l ≈ 6.6 at 0.85T _c , respectively. The size of the numerical simulation domain is 401 × 401, and the periodic boundary conditions are used in all directions. The time is dimensionless throughout the simulation time step. The bubble radius is normalized to the initial radius R ₀. As shown in Figure 2, the simulation results are consistent with the numerical results of the R-P equation. The correctness of the improved pseudopotential thermal MRT LB model is verified.

The temperature formula of the cavitation bubble, deduced from the R–P equation, is T = T ∞ ( R 0 R ) 3 ( γ − 1 ) , (23) where T ∞ indicates the ambient temperature of the liquid, and γ = 1.3 is a constant representing the adiabatic coefficient. In our simulation, periodic boundary conditions are adopted in all directions. The liquid temperature T ∞ = 0.85T _c is considered uniform. The pressure difference Δ p = p ∞ − p v . The instantaneous temperature T is normalized, and the t ^* is dimensionless throughout the simulation time step.

The current theoretical results show that the maximum collapse temperature of the bubble is infinite. What prevents such a runaway increase is heat dissipation, viscous dissipation and other uncontrollable factors in experiments or numerical simulations. As shown in Figure 3, the simulation results closely track the theoretical solution results, supporting the reliability of the present model. The thermal LB model can therefore be used to predict and analyze the collapse of cavitation bubble near solid wall.

As illustrated in Figure 4, numerical simulation is carried out in the two-dimensional x-y plane, and the computational domain l _x × l _y . The nonslip boundary scheme (Mohamad, 2017) is applied at the wall that is defined at y = 0 and periodic boundaries are utilized in the horizontal direction. The top of the computational domain is set as an open boundary, and the anti-bounce-back approach is applied. In our simulation, the initial radius of the bubble is R ₀, b denotes the distance from the center of the bubble to the solid wall, and λ = b R 0 is defined as the offset parameter. The initial equilibrium temperature is taken as T _∞ = 0.85T _c , corresponding to the coexistence densities ρ v ≈ 0.34 and ρ l ≈ 6.63.

Figure 5 shows the comparison between the current bubble profile and the results of an underwater pulsed discharge experiments published in the literature (Shan et al., 2019). Due to the difference in the size of the cavitation bubbles in the experimental study and numerical simulation, their collapse time cannot be directly compared and needs to be normalized. Comparing the Figures 5A,B, it can be found that the morphological evolution process of the bubble is consistent, and the dimensionless time obtained after normalization is the same. The result verifies the correctness of the LBM numerical experiment from the appearance change of the bubble.

In addition, in order to make a more accurate comparison and prove the symmetry of the collapse of cavitation bubble near the solid wall, the grid independence verification is carried out in this section. First the coordinates of the bubble surface need to be obtained. The criterion for the interface between the vapor and the liquid is that the density of the interface point is equal to ( ρ v + ρ l ) / 2 . We chose 301 × 301, 401 × 401 and 501 × 501 grids, according to the principle of the same scale, the initial bubble radii are 60, 80, and 100, respectively. The dimensionless time T ^* is defined, which represents the time from the initial state to the first collapse. The bubble contour points of 0, 0.2T ^*, 0.4T ^*, 0.6T ^*, 0.8T ^*, and T ^* under three grid resolutions are extracted, and the coordinate points obtained in the x and y directions are dimensionless, as shown in Figure 6. The result shows that the normalized contours of cavitation bubble collapse near the wall almost coincide with each other of the three grid resolutions, which indicates that the numerical model of cavitation bubble collapse near the wall is independent of the grid size. Considering the computational efficiency and the observability of the results, we choose 401 × 401 grids for the following numerical experiments.

Figure 7 shows the evolution of the temperature, pressure and fluid velocity fields of a single cavitation bubble collapse near the solid wall (prior to complete collapse). Due to the left-right symmetry of the entire system, in the evolution result image, the left side is the pressure field and the fluid velocity field, and the right side is the temperature field. As shown in Figures 7A,B, the spherical cavitation bubble is deformed into an elongated bubble in the direction perpendicular to the solid wall. The fluid velocity above the bubble is relatively high, which leads to a higher temperature in the liquid area and the area above the bubble. The near-wall fluid velocity is lower, causing the temperature in the bottom area of the bubble to decrease. As can be seen from Figures 7C,D, the micro-jet formed above the bubble causes local liquid pressure fluctuations and makes the local fluid particles vibrate rapidly. A thermal boundary layer forms at the bubble surface. Figure 7E shows the first collapse of the bubble, where the micro-jet impinges on the bottom edge of the bubble and the bubble temperature rises instantaneously. After the first collapse, a shock wave is emitted toward the solid wall and forms a hot spot on the wall, as has been discussed in Refs. (Yang et al., 2020) and (Yang et al., 2021). As shown in Figure 7F, the bubble becomes toroidal and its internal initial temperature is low. Then the toroidal bubble collapses rapidly under the driving action of the micro-jet, that is, the second collapse of the bubble. It can be seen from the time steps of Figure 7F-(i) that the second collapse of the bubble is very fast. The toroidal bubble temperature increases sharply and reaches a maximum when the bubble collapses completely.

From previous work, we know that the bubble temperature affected by the nondimensional offset parameter λ, initial pressure difference Δp and initial bubble radius R ₀. In order to understand the evolution results of the bubble temperature, the following sections discuss and analyze the above-mentioned factors in detail.

In this section, the effect of the different offset parameters on the bubble temperature will be discussed. Since the temperature distribution in the bubble is not uniform, we take a certain point inside the bubble as the temperature measurement point. For example, when λ = 1.2, as shown in Figure 8A, the coordinates of the measurement point are (201, 27). This is a point located at the bottom of the bubble on the vertical bisector of the x axis. For other values of λ, the temperature measurement point can be obtained by the same way. Figure 9 shows the bubble temperature curves with offset parameters of 1.2, 1.4 and 1.6. In order to explain the changes of the bubble temperature more clearly, the maximum jet velocity of the perpendicular bisectors of the x axis under three offset parameters is also shown in Figure 8. The initial conditions and parameter settings of the three cases are the same. In Figure 9A, the bubble temperature decreases slightly during the bubble shrinkage stage, which is caused by the lower pressure of the bottom liquid and lower velocity. With the formation of the micro jet, the maximum jet velocity accelerates gradually, and the bubble temperature also increases. Then, the jet velocity decreases after reaching its maximum, which is impeded by the deforming bubble. Meantime, the bubble temperature drops briefly and rises sharply at the first collapse. In Figure 9B, the trend of the bubble temperature curve is similar to that of Figure 9A, and the jet velocity decreases slightly before the first collapse. The maximum bubble temperature is higher than in Figure 9A. However, as shown in Figure 9C, the jet velocity keeps increasing in the early stages of the collapse, and the acceleration of the final stage of the collapse decreases. The bubble temperature increases as λ increases and reaches the greatest value when λ = 1.6. It can be inferred that the larger is λ, the higher is the maximum bubble temperature. In addition, we find that the cavitation bubble will become toroidal when λ = 1.2 after the first collapse. The bubble collapses directly without becoming toroidal when λ = 1.4 or 1.6 (i.e., when it is farther from the wall). As shown in Figure 8B, the coordinates of the temperature measurement point of the toroidal bubble are (183, 39). Figure 10 shows the evolution of the toroidal bubble temperature and the maximum jet velocity around the bubble from the second collapse to complete collapse. Due to the driving effect of the jet and the water hammer effect, the jet velocity increases rapidly, and the toroidal bubble collapses and produces high temperature almost instantaneously.

In this section, the influences of the initial driving pressure and bubble radius are studied. Because the collapse time is different for different Δp and different initial radius R0, the time step is normalized in the following section for comparison. Figure 11 shows the bubble temperature vs. time during the first collapse for three different values of initial Δp. The result shows that there is no obvious dependence of temperature on Δp at the bubble contraction stage. As Δp increases, the bubble temperature increases gradually. Then the bubble temperature drops briefly, as shown in more detail in Figure 9A, the bubble temperature is related to the jet velocity. Due to the impeding effect of the bubble, the jet velocity decreases briefly, causing the bubble temperature to drop. One may note that, the smaller Δp, the more significant the temperature drop. When the bubble collapses for the first time, its temperature increases dramatically. For Δp = 0.0897, the bubble temperature rises fastest and to the greatest value. It can be concluded that the large pressure difference has large potential energy, which releases more energy when collapsing leading to a higher temperature. The bubble becomes toroidal after the first collapse in all three cases. Figure 12 shows the time-dependence of the toroidal bubble temperature. The initial toroidal temperature is low; it increases gradually as the bubble collapses and eventually reaches a maximum. The collapse temperature is almost the same as that of the three Δp. It can be inferred that the pressure difference has little effect on the toroidal bubble temperature during the second collapse.

The evolutions of the bubble temperature during the cavitation bubble collapse with initial radii of 80, 85 and 90 and with an initial driving pressure Δp = 0.0853 are shown in Figure 13. It can be seen that the bubble temperature rises gradually after a brief cooling. When R ₀ = 90, the larger value, the bubble temperature drops significantly. The bubble temperature then increases sharply when the bubble collapses for the first time. The larger the bubble radius, the greater is the internal energy, and the more energy is released when the bubble collapses, which leads to higher bubble temperature.

The bubble becomes toroidal after the first collapse. In Figure 14, the bubble temperature of the newly formed toroidal bubble is the lowest among these three examples when R ₀ = 90, but for the value of R0 the final temperature is the highest, with the most drastic bubble temperature change. For larger initial radius, the maximum bubble temperature is therefore higher at the second collapse.