Study on the Collapse Process of Cavitation Bubbles Near the Concave Wall by Lattice Boltzmann Method Pseudo-Potential Model

: In this paper, the lattice Boltzmann pseudo-potential model coupled the Carnahan–Starling (C-S) equation of state and Li’s force scheme are used to study the collapse process of cavitation bubbles near the concave wall. It mainly includes the collapse process of the single and double cavitation bubbles in the near-wall region. Studies have shown that the collapse velocity of a single cavitation bubble becomes slower as the additional pressure reduces, and the velocity of the micro-jet also decreases accordingly. Moreover, the second collapse of the cavitation bubble cannot be found if the additional pressure reduces further. When the cavitation bubble is located in di ﬀ erent angles with vertical direction, its collapse direction is always perpendicular to the wall. If the double cavitation bubbles are arranged vertically, the collapse process of the upper bubble will be quicker, as the relative distance increases. When the relative distance between the bubbles is large enough, no second collapse can be found for the upper bubble. On the other hand, when two cavitation bubbles are in the horizontal arrangement, the suppression e ﬀ ect between cavitation bubbles decreases as the relative distance between the bubbles increases and the collapse position of cavitation bubbles moves from the lower part to the upper part.


Introduction
Cavitation is common in hydraulic and marine engineering, such as in water turbines and ship propellers [1]. It may cause great damage on the surrounding structures when the cavitation bubbles collapse in the near-wall region [2]. Therefore, it is vital to reveal the process of cavitation-bubble collapse and the effect of its interaction with the structures. For one thing, it can reduce the damage to the physical structure; for another, it can also be used to clean the surface of the structures.
The process of cavitation bubble collapse has been studied in the past years. Kornfeld and Suvorov firstly proposed the concept of micro-jet and believed that micro-jet was the main cause of structural damage [3]. There are many experimental results about the cavitation-bubble collapse near a solid wall. Naude' and Ellis discovered the micro-jet process in the experiment when the bubble in the near wall region [4]. Kling and Hammitt described the dynamic process of cavitation-bubble collapse in detail and the damage to the aluminum by experimental methods [5,6]. Vogel et al. utilized high-speed camera technology to accurately measure the micro-jet and the counter-jet velocity and carefully observed the evolution of the cavitation bubble in different distances between the bubble and the wall during the collapse [7]. Tomita and Kodama investigated a laser-induced cavitation-bubble collapse near a composite surface [8]. These experiments are mainly about the bubble collapse near a

Pseudo-Potential LBM-MRT
In the present study, the pseudo-potential LBM combining a multiple relaxation time (LBM-MRT) [15,24] and the external force scheme proposed by Li et al. Ref. [17] will be used. The governing equation [17,25] can be expressed as follows: where f α is the density distribution function, e α is the discrete velocity in the α direction, ∧ = M −1 ∧ M, is the relaxation matrix, M is the transformation matrix, ∧ is the diagonal matrix, f eq is the equilibrium distribution function, and S is the source term in the velocity space. For D2Q9 model, M is given by [26]: The diagonal matrix can be written as: Multiplying the both sides of the equation by M, the Equation (1) can be rewritten [27] as follows: I is the unit matric, m α = M αβ f β , m eq = M αβ f eq β , S = MS. The expression of m eq is given by: where ρ = α f α is the density and v = α f α e α /ρ + δ t F/2ρ is the macroscopic velocity. F is the total force, including the fluid-fluid interaction force, the solid-fluid interaction force, and the volume force. Because the key point of the present work is to study the interaction between bubbles, the fluid-solid interaction is not included; only the fluid-fluid interaction force is considered and is given by [28] as follows: where G is the interaction strength, φ(x) is the interaction potential [29]. ω is the weight coefficient. In the D2Q9 model, ω(1) = 1/3, ω(2) = 1/12: Energies 2020, 13, 4398 4 of 18 in which c = 1 is the lattice constant and c s = c/ √ 3 is the lattice sound speed. In the present study, the pressure p eos is calculated by Carnahan-Starling (C-S) equation of state, which is given by [30] as follows: R is the ideal gas constant, T is the temperature, b = 0.18727RT c /p c , a = 0.4963R 2 T 2 c /p c , T c is the critical temperature, and p c is the critical pressure. T c can be calculated by the parameter a, b and R.
The source term S in the Equation (4) is given by [17] as follows with a parameter σ, which effects the thermodynamic consistency and the stability of the model: The streaming process is given by

Model Validation
Firstly, the verification of the model is conducted through the Laplace's law. The pressure difference between the interior and the exterior of the droplet could be obtained by the Laplace's law: where p in is the pressure in the droplet and p out is the pressure out the droplet. σ is the surface tension and R is the radius of the droplet. In the simulation, a droplet is initially suspended at the center of the computational domain with initial radius, RR. The domain has a 501 × 501 lattice mesh system. The saturated temperatures of the droplet are 0.5 T c , 0.6 T c and 0.7 T c , respectively. The periodic boundary is implemented for the four sides of the domain. The terminal value of the R at the equilibrium state is slightly different from the RR. A series of RRs is selected in the calculation, so a series of relationships can be obtained between R and ∆p. The results are shown in Figure 1. Secondly, this LBM-MRT will be validated by the cavitation-bubble collapse process in the near-wall region. In this case, the initialized density field can be given by: where ρ l , ρ g are the liquid and the vapor density, respectively. x 0 , y 0 are the coordinates of the bubble's center. RR indicates the radius of the bubble. W denotes the width of the interface, which is set as 5 in all cases. The units in the present study are all denoted by the lattice units. For example, the mass unit is mu, the length unit is lu, the time unit is tu, the corresponding speed unit is lu·tu −1 , the density unit is mu·lu −3 , the pressure unit is mu·lu −1 tu −2 , and the lattice speed is c = lu/tu. Secondly, this LBM-MRT will be validated by the cavitation-bubble collapse process in the nearwall region. In this case, the initialized density field can be given by: where l g  are the liquid and the vapor density, respectively. 00 x , y are the coordinates of the bubble's center. RR indicates the radius of the bubble. W denotes the width of the interface, which is set as 5 in all cases. The units in the present study are all denoted by the lattice units. For example, the mass unit is mu , the length unit is lu , the time unit is tu , the corresponding speed unit is -1 lu tu  , the density unit is -3 mu lu  , the pressure unit is -1 -2 mu lu tu  , and the lattice speed is c=lu/tu .
In the numerical simulation, RR = 70 lu, the distance between the center of the cavitation bubble and the wall is 1.5 × RR. In the diagonal matrix: The parameters in the equation of state are set as follows: a = 0.2, b = 1, and T = 0.5 Tc. The initial pressure difference between the inside and outside of the cavitation bubble is ∆ = 0.00543 mu • lu −1 tu −2 . The above parameters are used for all cases in this paper. Figure 2 illustrates the comparison of the cavitation-bubble collapse process between the experimental [31] and the numerical results using the LBM-MRT. It can be seen that the simulation agrees well with the experimental results. The cavitation bubble shrinks under the effect of pressure difference in the initial stage. Obstructed by the straight wall, the speed of bubble bottom wall is small, and the lateral contraction is greater than that in longitudinal direction, shaping into an ellipse. As the pressure difference between inside and outside the bubble gradually increases, the depression appears and gradually deepens at the top of the cavitation bubble. Then, the first collapse occurs until the upper bubble wall contacts the lower bubble wall. In the numerical simulation, RR = 70 lu, the distance between the center of the cavitation bubble and the wall is 1.5 × RR. In the diagonal matrix: τ ρ = τ e = τ j = τ q = 1.0, τ ζ = τ j = τ q = 0.91, τ υ = 0.57. The parameters in the equation of state are set as follows: a = 0.2, b = 1, and T = 0.5 Tc. The initial pressure difference between the inside and outside of the cavitation bubble is ∆p = 0.00543 mu·lu −1 tu −2 . The above parameters are used for all cases in this paper. Figure 2 illustrates the comparison of the cavitation-bubble collapse process between the experimental [31] and the numerical results using the LBM-MRT. It can be seen that the simulation agrees well with the experimental results. The cavitation bubble shrinks under the effect of pressure difference in the initial stage. Obstructed by the straight wall, the speed of bubble bottom wall is small, and the lateral contraction is greater than that in longitudinal direction, shaping into an ellipse. As the pressure difference between inside and outside the bubble gradually increases, the depression appears and gradually deepens at the top of the cavitation bubble. Then, the first collapse occurs until the upper bubble wall contacts the lower bubble wall.    Figure 3a illustrates the shape of the bubble near a concave wall in the collapse process by boundary integral method [9]. The distance between the bubble center and the wall is 1.7 RR, and the origin of the coordinate is at X* (= x/RR) = 0, Y* (= y/RR) = 0. It can be observed that the LBM results are basically consistent with the BIM results, which are similar as cavitation bubbles collapse near a flat wall.   Figure 3a illustrates the shape of the bubble near a concave wall in the collapse process by boundary integral method [9]. The distance between the bubble center and the wall is 1.7 RR, and the origin of the coordinate is at X* (= x/RR) = 0,    Figure 3a illustrates the shape of the bubble near a concave wall in the collapse process by boundary integral method [9]. The distance between the bubble center and the wall is 1.7 RR, and the origin of the coordinate is at X* (= x/RR) = 0, Y* (= y/RR) = 0. It can be observed that the LBM results are basically consistent with the BIM results, which are similar as cavitation bubbles collapse near a flat wall.

Collapse Process of the Cavitation Bubble Near the Concave Wall
After the model validation, the collapse process of a single cavitation bubble with different angles and additional pressure will be studied in this section. Moreover, double cavitation bubbles with horizontal and vertical layouts will be simulated. The evolution of the pressure field, characteristic point pressure, and micro-jet will be discussed.

Evolution of the Single Bubble under Different Additional Pressures
The simulated layout is shown in Figure 4. The pressure boundary condition [32] is used for the upper boundary, and the half-way bounce-back boundary [33] is adopted for the concave wall.
In this section, the effect of additional pressure on the collapse process for a single cavitation bubble will be discussed. The bubble is located at λ = 1.6 × RR. Four cases will be studied, and the parameters are shown in Table 1.

Collapse Process of the Cavitation Bubble Near the Concave Wall
After the model validation, the collapse process of a single cavitation bubble with different angles and additional pressure will be studied in this section. Moreover, double cavitation bubbles with horizontal and vertical layouts will be simulated. The evolution of the pressure field, characteristic point pressure, and micro-jet will be discussed.

Evolution of the Single Bubble under Different Additional Pressures
The simulated layout is shown in Figure 4. The pressure boundary condition [32] is used for the upper boundary, and the half-way bounce-back boundary [33] is adopted for the concave wall.  . Simulated layout for Cases 1-4 (RR is the initial radius of bubble, λ is the distance between the center of the bubble and the concave wall, Pv is the vapor pressure, P∞ is the ambient pressure, θ is the angel to the vertical, and P is the center point at the concave wall).
Due to the irregular boundary at the bottom, some special treatments are needed to obtain the information at the bottom boundary. It should be illustrated that the concave wall in the simulation is not a smooth semicircle, which is composed by a series of polylines. The curvature is determined by a certain expression of coordinates, such as y = x 2 . A key problem for the fluid nodes that are near the boundary is the in all directions could not be known by streaming process. . Simulated layout for Cases 1-4 (RR is the initial radius of bubble, λ is the distance between the center of the bubble and the concave wall, Pv is the vapor pressure, P∞ is the ambient pressure, θ is the angel to the vertical, and P is the center point at the concave wall). In this section, the effect of additional pressure on the collapse process for a single cavitation bubble will be discussed. The bubble is located at λ = 1.6 × RR. Four cases will be studied, and the parameters are shown in Table 1.
Due to the irregular boundary at the bottom, some special treatments are needed to obtain the information at the bottom boundary. It should be illustrated that the concave wall in the simulation is not a smooth semicircle, which is composed by a series of polylines. The curvature is determined by a certain expression of coordinates, such as y = x 2 . A key problem for the fluid nodes that are near the boundary is the f α in all directions could not be known by streaming process. For the convenience, three different kinds of nodes will be introduced in the discussion, which are fluid nodes, boundary nodes, and solid nodes. The detail has been shown in Figure 5. Simulated layout for Cases 1-4 (RR is the initial radius of bubble, λ is the distance between the center of the bubble and the concave wall, Pv is the vapor pressure, P∞ is the ambient pressure, θ is the angel to the vertical, and P is the center point at the concave wall).
Due to the irregular boundary at the bottom, some special treatments are needed to obtain the information at the bottom boundary. It should be illustrated that the concave wall in the simulation is not a smooth semicircle, which is composed by a series of polylines. The curvature is determined by a certain expression of coordinates, such as y = x 2 . A key problem for the fluid nodes that are near the boundary is the in all directions could not be known by streaming process. For the convenience, three different kinds of nodes will be introduced in the discussion, which are fluid nodes, boundary nodes, and solid nodes. The detail has been shown in Figure 5. Unlike the fluid nodes, the boundary nodes could not get all populations in the streaming process. It is necessary to determine the unknown population firstly. For the center node, 4 , 7 , 3 could be known by streaming process according to Equation (10). However, 1 , 2 , 5 , 6 , 8 could not be known, because the solid nodes are not involved in calculation. Therefore, the additional constraints should be included in these nodes. For halfway bounce-back boundary, the unknown population f can be determined by the following: Unlike the fluid nodes, the boundary nodes could not get all populations in the streaming process. It is necessary to determine the unknown population firstly. For the center node, f 4 , f 7 , f 3 could be known by streaming process according to Equation (10). However, f 1 , f 2 , f 5 , f 6 , f 8 could not be known, because the solid nodes are not involved in calculation. Therefore, the additional constraints should be included in these nodes. For halfway bounce-back boundary, the unknown population f can be determined by the following: where i is the direction of unknown population and i is the opposite direction of i. f * is the population before streaming. the first collapse is about to occur. Obviously, with the decrease in the additional pressure around the bubble, the pressure at the top of the bubble also reduces, and the velocity of collapse becomes slower accordingly. In Case 1, the collapse process only takes about 900 tu. While in Case 4, the collapse process takes about 2000 tu. It is found that when the distance between the bubble and the wall is fixed, the pressure that is greater than a critical value is required to produce the second collapse. For example, there is no second collapse of the bubble for Case 4. The first and second collapses will generate pressure waves and both of them are much larger than the additional pressure. Generally speaking, as the additional pressure around the bubble decreases, the process of cavitation-bubble collapse is almost the same except for the increase in the collapse time. When ∆p is smaller than a critical value, the second collapse of the cavitation bubble cannot be found.   6 and 7 show the collapse process of the cavitation bubble under different additional pressures. Stage 1 refers to the beginning of the depress process, and Stage 2 refers to the moment that the first collapse is about to occur. Obviously, with the decrease in the additional pressure around the bubble, the pressure at the top of the bubble also reduces, and the velocity of collapse becomes slower accordingly. In Case 1, the collapse process only takes about 900 tu. While in Case 4, the collapse process takes about 2000 tu. It is found that when the distance between the bubble and the wall is fixed, the pressure that is greater than a critical value is required to produce the second collapse. For example, there is no second collapse of the bubble for Case 4. The first and second collapses will generate pressure waves and both of them are much larger than the additional pressure. Generally speaking, as the additional pressure around the bubble decreases, the process of cavitation-bubble collapse is almost the same except for the increase in the collapse time. When ∆ is smaller than a critical value, the second collapse of the cavitation bubble cannot be found.   Figure 8 shows the evolution of the pressure at point P under different additional pressures. The pressure generated by the collapse spreads in all directions and will arrive at the concave wall finally. As the pressure around the bubble decreases, the collapse velocity reduces. Correspondingly, it will spend more time spreading to the wall for the pressure wave. Furthermore, the pressure peak decreases from 0.035 mu • lu −1 tu −2 in Case 1 to 0.018 mu • lu −1 tu −2 in Case 4. It reaches the peak almost instantaneously when the pressure acts on the wall and then begins to fall.  finally. As the pressure around the bubble decreases, the collapse velocity reduces. Correspondingly, it will spend more time spreading to the wall for the pressure wave. Furthermore, the pressure peak decreases from 0.035 mu·lu −1 tu −2 in Case 1 to 0.018 mu·lu −1 tu −2 in Case 4. It reaches the peak almost instantaneously when the pressure acts on the wall and then begins to fall.  Figure 8 shows the evolution of the pressure at point P under different additional pressures. The pressure generated by the collapse spreads in all directions and will arrive at the concave wall finally. As the pressure around the bubble decreases, the collapse velocity reduces. Correspondingly, it will spend more time spreading to the wall for the pressure wave. Furthermore, the pressure peak decreases from 0.035 mu • lu −1 tu −2 in Case 1 to 0.018 mu • lu −1 tu −2 in Case 4. It reaches the peak almost instantaneously when the pressure acts on the wall and then begins to fall.  Figure 9 shows the evolution process of micro-jet when cavitation-bubble collapse under different additional pressures. In the present study, the measurement of the micro-jet by Plesset and Chapman is adopted [34]. This method can effectively avoid the effect of virtual velocity in the pseudo-potential model by treating the velocity of the bubble wall at the top of the cavitation bubble as the micro-jet velocity. When the additional pressure around the bubble decreases, the velocity of  Figure 9 shows the evolution process of micro-jet when cavitation-bubble collapse under different additional pressures. In the present study, the measurement of the micro-jet by Plesset and Chapman is adopted [34]. This method can effectively avoid the effect of virtual velocity in the pseudo-potential model by treating the velocity of the bubble wall at the top of the cavitation bubble as the micro-jet velocity. When the additional pressure around the bubble decreases, the velocity of the bubble collapse also reduces due to the micro-jet's velocity decreases and the peak of the micro-jet decreases accordingly. However, the peak occurs during the depression of the bubble rather than the period of the first collapse. Moreover, there is a slight attenuation after micro-jet velocity reaches the maximum. In Case 4, due to the lower additional pressure around the bubble, the bubble collapse pattern is different from the other cases. For example, there is no extremum for the micro-jet.
Energies 2020, 13, x FOR PEER REVIEW 10 of 20 the bubble collapse also reduces due to the micro-jet's velocity decreases and the peak of the microjet decreases accordingly. However, the peak occurs during the depression of the bubble rather than the period of the first collapse. Moreover, there is a slight attenuation after micro-jet velocity reaches the maximum. In Case 4, due to the lower additional pressure around the bubble, the bubble collapse pattern is different from the other cases. For example, there is no extremum for the micro-jet.

Evolution of the Single Bubble in Different Angles with Vertical Direction
In this section, the effect of the cavitation bubble angle on the evolution process of collapse will be discussed and the used parameters for the Cases 5 and 6 are shown in Table 2. Table 2. Parameters of single bubble in Cases 5-6.

Evolution of the Single Bubble in Different Angles with Vertical Direction
In this section, the effect of the cavitation bubble angle on the evolution process of collapse will be discussed and the used parameters for the Cases 5 and 6 are shown in Table 2.  Figures 10 and 11 illustrate the collapse process when the cavitation bubble is in different angles with vertical direction. In the initial stage, the external pressure of the cavitation bubble is larger than the pressure inside the bubble. Therefore, the cavitation bubble shrinks overall under the pressure difference, and a low-pressure area is formed in the near-wall zone. The moving speed at the bottom of the bubble is lower than other positions on the interface. It evolves subsequently into an ellipse due to the blockage of the wall (t = 500 tu). After that, the pressure on the upper surface of the bubble continues to increase, causing the bubble to depress and forming a micro-jet. Unlike the straight wall, the bubble depression direction in Cases 5 and 6 is not vertical but perpendicular to the concave wall. Under the effect of the pressure difference, the bubble forms the first collapse at t = 910 tu, and the collapse pressure is generated. At this time, the cavitation bubble evolves into a ring shape. The collapse pressure is much greater than the surrounding pressure, and it could reach 0.045 mu·lu −1 tu −2 . The ring-shaped bubble finally collapses under the combined influence of the collapse pressure and the surrounding pressure, which is called the second collapse. The pressure caused by the second collapse is smaller than the first, which is only 0.03 mu·lu −1 tu −2 . The direction of the micro-jet velocity changes due to the blockage of the wall surface after the second collapse, resulting in the formation of a vortex around the ring-shaped bubble. Moreover, the vortex lasts for quite a long time, and its existence leads to a low-pressure zone.

Evolution of the Double Bubbles with Vertical Arrangement
In this section, the collapse process of double cavitation bubbles with vertical arrangement is studied, and the effect of the relative distance between the double cavitation bubbles on the collapse process is discussed. Figure 12 illustrates the computation layout. The parameters and the computational domain are the same as in the previous cases, some of which are shown in Table 3.

Evolution of the Double Bubbles with Vertical Arrangement
In this section, the collapse process of double cavitation bubbles with vertical arrangement is studied, and the effect of the relative distance between the double cavitation bubbles on the collapse process is discussed. Figure 12 illustrates the computation layout. The parameters and the computational domain are the same as in the previous cases, some of which are shown in Table 3.      Figure 13 displays the cavitation-bubble's collapse for Cases 7-9. It can be observed that the upper bubble is far away from the rigid wall, so the effect of the wall on it can be negligible. The upper bubble shrinks under the surrounding additional pressure and then sags and collapses. This process is similar to the collapse of the single bubble near the rigid wall. As λ 2 increases, the pressure above the upper bubble increases, and the velocity of collapse also raises accordingly. In Case 9, no second collapse of upper cavitation bubble can be found. From these three cases, it can be seen that the effect of the lower bubble on the upper cavitation bubble is similar to the rigid wall. With the increase in λ 2 , the interval between the first and the second collapses of the upper bubble becomes shorter; eventually there is no occurrence of the second collapse. In Case 7, there is no obvious displacement at the top and bottom of the lower bubble, and the lateral velocities on both sides are relatively large. The lower cavitation bubble is elongated along the vertical direction. It seems that the influence of the upper bubble and the wall on it reaches approximate balance. The lower bubble tends to sag from the middle (t = 1110 tu), which is similar to the collapse of the bubble with both rigid walls on the top and the bottom boundary [35]. It can be observed from Case 7 that a vortex forms after the upper bubble collapses. and the pressure in the  In Case 7, there is no obvious displacement at the top and bottom of the lower bubble, and the lateral velocities on both sides are relatively large. The lower cavitation bubble is elongated along the vertical direction. It seems that the influence of the upper bubble and the wall on it reaches approximate balance. The lower bubble tends to sag from the middle (t = 1110 tu), which is similar to the collapse of the bubble with both rigid walls on the top and the bottom boundary [35]. It can be observed from Case 7 that a vortex forms after the upper bubble collapses. and the pressure in the vortex zone is negative. In other cases, the pressure in that zone is relatively close to that in the surrounding environment. The reason for that is mainly because the relative distance between the bubbles is rather small in Case 7, and the upper bubble is obviously affected by the obstruction of the lower one. Therefore, the direction of the micro-jet changes accordingly, resulting in a vortex zone, which generates negative pressure later.
After the upper bubble collapses, the collapse process of the lower one in each case is not the same. For example, the mutual effect between cavitation bubbles decreases as λ 2 increases. Besides, the shrinking center of the lower bubble shifts from its middle to the upper. Figure 14 shows the evolution of pressure at point P in Cases 7-9. It can be seen that the evolution of pressure at point P for Cases 7 and 8 is similar. For example, the peaks are all round 0.03 mu·lu −1 tu −2 , and the peaks appear at almost the same time. As the distance between the upper cavitation bubble and the wall increases, the time of arriving at the wall for the pressure wave becomes shorter as well. In Case 9, the pressure peak is rather smaller than that in other cases, which is only 0.018 mu·lu −1 tu −2 . It is found that the pressure waves generated by the collapse of two bubbles counteracts. It can be seen that when the relative distance of cavitation bubbles is large enough, the wall surface can be protected.
Energies 2020, 13, x FOR PEER REVIEW 14 of 20 lu −1 tu −2 . It is found that the pressure waves generated by the collapse of two bubbles counteracts. It can be seen that when the relative distance of cavitation bubbles is large enough, the wall surface can be protected.  Figure 15 shows the evolution of the cavitation bubble micro-jet in Cases 7-9. In the initial stage, the cavitation bubble is in the contraction stage, and the velocity of the bubble's upper surface is the same in all cases. It has a slight difference in the stage of the depression when a micro-jet generates. Around t = 850 tu, as the relative distance of the cavitation bubbles increases, the growth rate of the micro-jet velocity also increases, and its peak ascends accordingly. After the micro-jet velocity reaches the peak, a slight downward trend can be found, and that is similar to the single bubble.  Figure 15 shows the evolution of the cavitation bubble micro-jet in Cases 7-9. In the initial stage, the cavitation bubble is in the contraction stage, and the velocity of the bubble's upper surface is the same in all cases. It has a slight difference in the stage of the depression when a micro-jet generates. Around t = 850 tu, as the relative distance of the cavitation bubbles increases, the growth rate of the micro-jet velocity also increases, and its peak ascends accordingly. After the micro-jet velocity reaches the peak, a slight downward trend can be found, and that is similar to the single bubble.  Figure 15 shows the evolution of the cavitation bubble micro-jet in Cases 7-9. In the initial stage, the cavitation bubble is in the contraction stage, and the velocity of the bubble's upper surface is the same in all cases. It has a slight difference in the stage of the depression when a micro-jet generates. Around t = 850 tu, as the relative distance of the cavitation bubbles increases, the growth rate of the micro-jet velocity also increases, and its peak ascends accordingly. After the micro-jet velocity reaches the peak, a slight downward trend can be found, and that is similar to the single bubble.

Evolution of the Double Bubbles with Horizontal Arrangement
In this section, the collapse process of double bubbles with horizontal arrangement will be simulated, and the effect of the relative distance between the bubbles on the collapse process will be discussed. The sketch of computational layout is shown in Figure 16, and the used parameters are displayed in Table 4.

Evolution of the Double Bubbles with Horizontal Arrangement
In this section, the collapse process of double bubbles with horizontal arrangement will be simulated, and the effect of the relative distance between the bubbles on the collapse process will be discussed. The sketch of computational layout is shown in Figure 16, and the used parameters are displayed in Table 4.     Figure 17 shows that the collapse process of double cavitation bubbles with horizontal arrangement. It can be seen that cavitation bubbles shrink under the effect of environmental pressure at the beginning, and then cavitation bubble sags diagonally under the joint action of the wall and the adjacent bubble, generating a micro-jet. As the relative distance of cavitation bubbles increases, the position of the Energies 2020, 13, 4398 15 of 18 depression gradually moves upward. With the deepening of the depression, the effect of the relative distance between the bubbles on the collapse process can be clearly observed. In Case 10, the bubble collapses from the lower part at first, and it disappears rapidly under the effect of the collapse pressure generated by the lower part. Then a low-pressure zone generates where a vortex is formed. In Case 11, the cavitation bubbles form the first collapse from the middle, and it indicates that the suppression effect between cavitation bubbles is nearly equal to the inhibition effect of the wall on the bubble. In Case 12, the upper part collapses first, and the bubbles disappear under the collapse pressure. In Case 10, the pressure in the vortex zone is extremely low, and it is only about 0.001 mu·lu −1 tu −2 . The pressure in these low-pressure zones for Case 12 is around 0.005 mu·lu −1 tu −2 , which is bigger than that in Case 10. Therefore, it can be seen that the pressure in the vortex zone after the collapse is rather lower when the cavitation bubbles are relatively close. In addition, as the distance between them increases, the pressure in this zone increases accordingly.

Discussion
In this section, a brief comparison is displayed between results by other authors about the cavitation bubbles and the simulation in this paper. Mao et al. investigated the single and the doublebubbles' collapse near a flat wall by single relaxation time lattice Boltzmann method (SRT-LBM) [19]. A series of collapse laws have been found by analyzing the density and pressure field in different initial conditions. Tomita studied the effect of the rigid surface curvature on the bubble behavior by

Discussion
In this section, a brief comparison is displayed between results by other authors about the cavitation bubbles and the simulation in this paper. Mao et al. investigated the single and the double-bubbles' collapse near a flat wall by single relaxation time lattice Boltzmann method (SRT-LBM) [19]. A series of collapse laws have been found by analyzing the density and pressure field in different initial conditions. Tomita studied the effect of the rigid surface curvature on the bubble behavior by BIM [9]. In addition, they observed a mushroom-shaped bubble during the bubble collapse stage and found that the boundary curvature increases the micro-jet velocity. Shervani-Tabar and Rouhollahi investigated a single-bubble collapse near a concave rigid wall by two different methods, including BIM and finite difference method [36]. They found that the micro-jet velocity increased with the decrease in the surface concavity. Xue et al. numerically studied a single bubble collapse near a convex wall in different curvatures by multi-relaxation time lattice Boltzmann method (MRT-LBM) and obtained a relationship between micro-jet velocity and the initial pressure differences [37]. In the present study, an improved MRT-LBM model has been used to investigate the collapse of single and double bubbles near a concave rigid wall. It is found that the discrepancy exists both for the bubbles' behavior in different angels and the interactions between bubbles. More details of the work in this paper can be seen in the conclusion.

Conclusions
In this paper, the LBM-MRT model is used to simulate the collapse process of cavitation bubbles near the concave wall. Firstly, the model is verified by the experiment results. After that, the model is used to simulate the collapse process of cavitation bubbles near the concave wall. It includes the cases of different additional pressures and angles with vertical direction, as well as double cavitation bubbles with vertical and horizontal arrangements. Moreover, the evolution of the pressure field, the characteristic point pressure, and the micro-jet have been discussed in detail. Based on the studies, the following conclusions can be drawn:

1.
The collapse process of cavitation bubble is affected by the pressure of the surrounding environment. When the additional pressure around the environment decreases, the velocity of cavitation-bubble collapse becomes slower, and the duration of the collapsing process increases accordingly. Moreover, no second collapse of cavitation bubble can be found when the additional pressure is lower than a critical value. When the angle of the cavitation bubble with vertical direction changes, the collapse process of cavitation bubble is similar, but the depression direction is perpendicular to the concave wall. After it collapses, the low-pressure zone is generated due to the vortex.

2.
When the double cavitation bubbles are arranged vertically, as the relative distance between cavitation bubbles increases, the pressure above the upper bubble increases, and the velocity of collapse also raises accordingly. With the increase in the relative distance, the interval between the first and the second collapses of the upper bubble becomes shorter; eventually there is no occurrence of the second collapse. After the upper bubble collapses, the collapse process of the lower one in each case is not the same. For example, the mutual effect between cavitation bubbles decreases as the relative distance increases. Besides, the shrinking center of the lower bubble shifts from its middle to the upper.

3.
When the double cavitation bubbles are arranged horizontally, the mutual effect between the cavitation bubbles gradually decreases, as the relative distance of the cavitation bubbles increases. The depression position of cavitation bubble gradually moves from the lower part to the upper part, as the relative distance increases. In Case 11, the cavitation-bubble collapse from the middle of the bubble under the interaction between cavitation bubbles and the influence of the wall on the bubbles. After the collapse, the pressure in the vortex zone increases accordingly, as the relative distance of the cavitation bubbles increases.