Simulation and analysis of collapsing vapor-bubble clusters with special emphasis on potentially erosive impact loads at walls

Cavitation is a process where a liquid evaporates due to a pressure drop and re-condenses violently. Noise, material erosion and altered system dynamics characterize for such a process for which shock waves, rarefaction waves and vapor generation are typical phenomena. The current paper presents novel results for collapsing vapour-bubble clusters in a liquid environment close to a wall obtained by computational fluid mechanics (CFD) simulations. The driving pressure initially is 10 MPa in the liquid. Computations are carried out by using a fully compressible single-fluid flow model in combination with a conservative finite volume method (FVM). The investigated bubble clusters (referred to as “clouds”) differ by their initial vapor volume fractions, initial stand-off distances to the wall and by initial bubble radii. The effects of collapse focusing due to bubble-bubble interaction are analysed by investigating the intensities and positions of individual bubble collapses, as well as by the resulting shock-induced pressure field at the wall. Stronger interaction of the bubbles leads to an intensification of the collapse strength for individual bubbles, collapse focusing towards the center of the cloud and enhanced re-evaporation. The obtained results reveal collapse features which are common for all cases, as well as case-specific differences during collapse-rebound cycles. Simultaneous measurements of maximum pressures at the wall and within the flow field and of the vapor volume evolution show that not only the primary collapse but also subsequent collapses are potentially relevant for erosion.


Introduction
Collapsing single bubbles have been extensively investigated in the past and results obtained by experiments, theory and simulations are broadly represented in the literature [1], [2]. Significantly fewer publications can be found where larger bubble clusters are considered [3][4][5]. Although there is a simple theory provided by Brennen et al [6], clusters of bubbles are difficult to investigate experimentally. One reason is that a repeatable generation of a cluster of well-defined bubbles (radii, position) is hardly achievable. Therefore, numerical simulation of collapsing bubble clusters is a suitable way to enhance physical insight into collapse processes.
To estimate the interaction among bubbles a cloud interaction parameter β was introduced [4]. This parameter is defined as α0·(1-α0)·A 2 0 /R 2 0, where α0 is initial cloud void fraction and A0 and R0 are the cloud and the bubble radius, respectively.
It was shown in previous studies [7], [8] that cloud collapse mechanisms are significantly depending on the cloud interaction parameter. For high values (β >>1) the shock develops and propagates from the outside towards the center of a cloud, therefore producing a focused collapse. At low values (β < 1) the interaction of the bubbles is weak and one might rather observe individual collapses. In this investigation the interaction parameter is chosen to be 4, 8 and 11, indicating that bubble interaction is relevant.
There are several common approaches to model liquid-vapor interfaces, for example: two-fluid approaches [9] with sharp interface treatment [1], the volume of fluid technique [10], and so called single-fluid approaches with implicit interface capturing techniques. Two-fluid methods are capable of representing the phase interphase with high precision and allow for inclusion of surface tension. However, they require the solution of two sets of governing equations and an explicit interfaceinteraction model.
In the present work a single-fluid approach coupled with a thermodynamic equilibrium model is adopted. When the grid resolution is sufficiently fine, the thermodynamic equilibrium model implicitly resolves the interface and becomes comparable to a sharp interface technique without surface tension [3]. It enables the prediction of re-evaporation processes, such as bubble rebound. However, the inclusion of surface tension is difficult and thus, the method is suitable for processes where surface tension is negligible, such as the investigated inertia dominated collapses of bubble clusters under high ambient pressure.
In this investigation, two major physical aspects, namely the cloud interaction parameter and the standoffdistance, are investigated.
First, a series of three bubble clusters at a fixed stand-off distance from a wall are generated. The radii of the approximately spherical clusters are kept constant while the common radii of the involved bubbles are varied. This leads to three different cloud interaction parameters.
In a second series of computations, the cloud interaction parameter is kept constant while four different stand-off distances are investigated.
For all configurations, the dominant collapse mechanism, rebound structures and maximum impact loads are analysed. A major finding is that the secondary collapse after the first rebound may lead to the highest loads at the solid wall. This indicates that the secondary collapse can considerably contribute to cavitation induced material erosion.

Numerical approach
The main focus of this investigation is on inertia-driven effects and wave dynamics arising during the collapse of a bubble cloud. Surface tension and viscous effects are neglected and the governing equations are the compressible Euler equations. For this investigation the flow simulation code CaTUM is used [11]. All numerical procedures have been extensively validated and the results have been published in the literature [12][13][14][15][16]. CaTUM is a density based finite volume method employing a Low-Mach-number consistent flux function and an explicit time marching procedure. The compressibility of the fluid is taken into account in order to capture shock formation and wave propagation. Hence the numerical time step is proportional to the ratio of the minimum grid size and the fastest signal speed (speed of sound in the liquid).
Phase change is modelled with equilibrium phase transition, i.e. both phases are in thermal, mechanical and phase equilibrium. In this study, the thermodynamic properties of water (in its liquid and vapor phases) are described by efficient barotropic equations of state. Pure liquid can be modelled by a Tait equation while saturated mixtures are modelled by a barotropic relation p= p (ρ), which is obtained by integration of the equilibrium speed of sound along an isentrope. The vapor volume fraction of a saturated mixture of liquid and vapor can be derived from the saturation densities of both phases.

Cloud specification and computational domain
As there is insufficient experimental data of "reference clouds" available, a random procedure proposed by Schmidt et al [3] is utilized to generate bubble clusters. The following criteria are assumed: The bubble cluster covers a spherical domain with a diameter of 30 mm. Within this domain, 150 spherical bubbles of equal radii are randomly generated. The minimum distance between neighbouring bubbles is one bubble radius. It is further assumed that the bubble number density reaches its maximum at the center of the cloud and follows a K/r rule. Here, K is a positive constant dependent on the selected bubble radius and r is the radial coordinate from the center to the outside of the bubble cluster.
Numerical tests with 100.000 bubbles showed that the random bubble generator produces the desired discrete bubble distribution consistently with the analytic model. Figure 1 shows a comparison of the distribution achieved by the random generator in comparison with the analytic model.
The generated cloud is embedded within a small computational domain of 35x35x40 mm³ (see Fig. 2) which itself is embedded in a larger domain of 4x4x2 m³. In the small domain a fine Cartesian grid is applied in order to minimize numerical errors.  All bubbles are filled with water vapor while the surrounding domain contains liquid water at an initial pressure of p∞ = 100 bar. The initial pressure inside the bubbles is equal to the vapor pressure psat = 2340 Pa at T = 293K. The velocity field is initially at rest.
The effects of surface tension, viscosity, gravity and non-condensable gas content are neglected.

Recorded simulation data
A numerically defined "pressure transducer" is used to record the static pressure on an area of 15x15 mm 2 at each time step. The numerical transducer is located directly below the bubble cluster at the center of the bottom wall. The numerical sampling frequency of the transducer is 49.3·10 6 Hz. Aside from the transducer at the wall, the maximum pressure within the entire flow field is monitored.
By monitoring the total vapor volume Vvap within the entire domain the temporal evolution of evaporation and condensation can be analysed. Instead of showing vapor volume versus time directly, it is more intuitive to present the vapor volume through an equivalent radius Requiv. This is achieved by defining Vvap: = 4/3·π·R³equiv.
Collapsing bubbles or vapor patterns arising from rebounds are tracked by application of a "Collapse detector approach" developed by Mihatsch et al [13]. This algorithm detects isolated collapses during the simulation and characterizes the strength of a collapse by recording the maximum pressure at the focal point. Additionally, the positions and time instants where collapses have been detected are stored. One application of this "Collapse detector" is to estimate the potential of cavitation erosion by analysing strong collapses close to a material surface. [13][14][15].

Problem setup
The aim of this investigation is to examine the influence of two parameters (cloud interaction parameter β and stand-off distance d) with respect to the collapse process and on maximum loads at the wall. Therefore, two setups are generated and simulated. For both setups and all computations conducted, the radius of the generated cloud is A0 = 15mm.
The first setup contains three random clouds with different initial bubble radii of R0 = 0.5 mm, 1.0 mm and 1.5 mm, respectively. As the number of bubbles equals 150 for all cases, the initial vapor volume fraction α0 and the interaction parameter β depend on the initial bubble radius only. The initial vapor volume fraction is computed by α0 = 150·Vbubble/Vcloud and the interaction parameter is β = 4, 8 and 11 respectively. As the distribution of bubbles is randomly generated for each configuration, the positions of individual bubbles are case dependent. The stand-off distance d = 3.2 mm of the cloud from the wall is kept constant.
The second setup is designed to investigate the effects of different stand-off distances to the wall. The clouds are identical in terms of bubble distribution and initial bubble radius, which is R0 = 1 mm. The stand-off distance is defined as the minimum distance between the virtual surface of the cloud and the bottom wall. In the following, stand-off distances of d = 1.2 mm, 3.2 mm, 6.0 mm and 10.0 mm are applied.
Setup details are shown in the Table 1 and visualized in Figure 3.

Results
This section describes major results where the prior interest is to analyse collapse mechanisms depending on the investigated configuration. The analysis contains the pressure field within the domain, the maximum load at the wall, and the temporal evolution of the amount of vapor. In particular, attention is paid to secondary clouds caused by rebounds and the location and intensity of their collapses.

Setup I: Variation of interaction parameter β
An increase of the bubble radius leads to an increase of the cloud interaction parameter at constant cloud radius and a constant number of bubbles. Therefore, the variation of this quantity will affect the cloud collapse behaviour, e.g. the strength of bubble-bubble interaction, focussing effects towards the center, shielding effects [4], as well as the intensities of collapses of individual bubbles and of the focused collapse at the center of the bubble cluster. Figure 4 shows the dimensionless temporal evolution of the equivalent radius Requiv as well as of the output from the pressure transducer at the wall for all three interaction parameters. Note that the Rayleigh time τ [17] used to nondimensionalize the time axis has to be computed individually for each configuration based on the case-specific vapor volume in order to obtain a similarity of the temporal evolution of the equivalent radius. It can be seen that the relative duration of the cloud collapse is about 60% of the case-specific Rayleigh time for all three configurations. However, the subsequent re-evaporation (rebound) is quite different between the three cases. The smaller the interaction parameter, the more pronounced is the delay between collapse and rebound. This tendency seems to hold for subsequent collapses of re-evaporated structures as well. As a result, the number of rebounds during unit time (nondimensionalized with the Rayleigh time) increases with increased interaction parameter.
With respect to the evolution of the pressure recorded at the wall, the following trends are found. Weak interaction leads to a low frequency increase of the pressure. Higher interaction causes high frequency pressure oscillations and larger amplitudes (see Fig. 4). However, a comparison of the maximum pressure recorded by the numerical transducers subsequent to the collapse phases show an unexpected behaviour, as In case of intermediate interaction (here for β = 8) a nonlinear trend is predicted. In this case, the most intense load at the wall occurs after the secondary collapse In the following, the physical details during the collapse and rebound phases are discussed for each configuration individually.

Weak interaction at β = 4
The first cloud contains the least amount of vapor of α0 = 0.004 and the initial bubble radius is R0 = 0.5 mm. Figure 6 shows the output of the collapse detector. Size and color of a diamond indicate maximum pressure of a collapse and its position. The cloud starts to collapse from the outer layer where individual bubble collapses with intensities of 10 4 bar. Bubbles continue to collapse independently in the whole cloud volume without pronounced interaction. All collapses are of comparable strength and there is no focussing or shielding effect visible.
The first rebound is rather weak and does not result in a pronounced secondary cloud (Fig. 7a). The collapse of the original cloud initiates a velocity field towards the center of the cloud. As the lower part of the cloud is confined by the adjacent wall, the upward motion is slightly slower than the downward motion. Thereby, the center of the rebound structure is shifted towards the wall. Consequently, the collapse of this structure takes place closer to the wall as compared to the original cloud (Fig. 7d). The pressure peak reaches 1.1·10 4 bar. The previously described processes reappear, leading to a rebound structure that finally collapses directly at the wall (Fig. 7g). The intensity of this collapse reaches 0.4·10 4 bar.
In despite of the decrease in pressure at the focal points of the collapses, the respond at the wall remains more or less the same. This explains the nonlinear behavior of the recorded maximum pressure at the wall as already shown in Figure 5. The shift of the focal point compensates the weakening of the collapse intensity, leading to almost constant loads at the wall.

Intermediate interaction at β = 8
The second cloud contains an intermediate amount of vapor of α0 = 0.03. The initial bubble radius is R0 = 1.0 mm.
Compared to the previous case with weak interaction, the intensities of individual collapses increase up to 1.5·10 4 bar. Individual collapses appear mainly in the outer region of the cloud and a pronounced focusing of the collapse towards the center is detected. The output of the collapse detector (Fig. 6b) shows two strong collapses in the center region with intensities of 1.9 and 2.1·10 4 bar.
After the primary collapse the secondary cloud appears as a result of rarefaction waves which are reflected from the wall. Thereby, the pressure reduction is amplified by the wall and re-evaporation takes place. Most of the resulting vapour structures are in contact to the wall (Fig. 7b). The collapse of these vapour structures leads to a pressure peak of 2.2·10 4 bar (Fig. 7e). This kind of impact might be highly erosive.
The following rebound appears again at the wall (Fig. 7h) and the subsequent collapse exhibits a relatively moderate magnitude of 0.8·10 4 bar.

Strong interaction at β = 11
The cloud featuring the largest amount of vapour among the considered clouds exhibits an initial vapor volume α0 = 0.08 and the initial bubble radius is R0 = 1.5mm.
Strong focussing is observed where the intensities of individual collapses are successively increasing towards the center of the cloud (Fig. 6c). The pressure peak at the focal point reaches 2.6·10 4 bar. In this case bubblebubble interaction and shielding is significant. A snapshot of rebounded single bubbles and a pronounced a -c: Pressure field and vapor iso-surface after the prior collapse before the second. In the first case no pronounced cloud appears during rebound. In the other two cases vapor rebound occurs right at the wall, a secondary cloud in the center of the domain. Clouds collapse from top to bottom in both latter cases.

d -f
: Pressure field and vapor iso-surface at the second collapse. In the first case relatively small collapse of 1.1·10 4 bar appears above the wall. In the middle case vapor at the wall collapses producing pressure peak of 1.78·10 4 bar. The last case shows the most violent collapse at the wall of 2.26·10 4 bar. This type of collapse is supposed to be highly relevant for cavitation erosion.
g -i: Pressure field and vapor iso-surface at the third collapse. In the first case the position of a collapse shifts further towards the wall and it occurs right at the wall with pressure magnitude 0.4·10 4 bar. In the middle case the third collapse of 0.8·10 4 bar is also at the wall. In opposite, the last case shows no significant impact on the wall since vapor structures are driven far above the wall and then collapse with 0.07·10 4 ·bar. Two latter clouds collapse from bottom to top.  Figure 8.
It can be seen that bubbles exhibit a toroidal shape due to liquid jet impingement during the collapse. The collapse of the secondary structure is focussed towards the wall (Fig. 7c). Thereafter, the most severe collapse occurs at the wall and reaches an intensity of 3.5·10 4 bar (Fig. 7f). The third cloud mainly consists of very weak vapor toroids. In this case, the onset of the collapse is located at the wall. The collapse front travels upwards without intensification until it terminates with a perceptible pressure peak of 0.07·10 4 ·bar (Fig. 7i). .

Setup II: Variation of stand-off distance d
Setup II consists of four identical clouds regarding the initial size and position of bubbles, which are located at four different distances from the wall: d = 1.2, 3.2, 6.0 and 10.0 mm. The distance d is defined as the minimum space between the wall and the virtual surface of the cloud. The cloud interaction parameter is β = 8 for all cases. Figure 9 shows the dimensionless temporal evolution of the equivalent radius Requiv as well as of the output from the pressure transducer at the wall for all four stand-off distances. In contrast to the corresponding analysis for setup I, the Rayleigh time τ used to nondimensionalize the time axis is in this case identical for all four configurations. Again, the duration of the collapse of the primary cloud is approximately 60% of the Rayleigh time of an equivalent bubble for all cases. It can be seen that the duration of the primary collapse is slightly shorter when the stand-off distance is increased. The wall partially confines the flow field for small stand-off distances, while it has minor effects when the cloud is located sufficiently far away from it. As β = 8, moderate interaction intensity and collapse focussing occurs for all cases. An unexpected behaviour is seen for the second and the third rebound of vapour structures. The maximum amount of vapour is inversely proportional to the initial stand-off distance during the first rebound. This is, however, not the case for the second rebound, where the tendency seems to be reversed.

Pressure on the wall
The output of the numerical pressure transducer at the wall shows that the pressure increases during the end of a collapse phase and features several maxima during collapse and rebound (Fig. 9). This observation holds for all four stand-off distances and all computed collapserebound-cycles. Figure 10 shows maximum pressures recorded at the wall after three collapse phases for all investigated stand-off distances. As expected, the intensity of the primary collapse is inversely proportional to the initial stand-off distance. This trend can be explained by the location of the focal point of the collapse, which locates more closely to the wall if the stand-off distance is smaller. Interestingly, the pressures recorded after the secondary collapses and especially subsequent to the third collapses show a nonlinear behaviour: for a standoff distance of d = 1.2 mm a decrease in pressure is predicted for subsequent collapses. In case of d = 3.2 mm to d = 10 mm the pressures recorded after the secondary collapses exceed those recorded after the primary collapses. Although the collapse intensities decrease when comparing collapses 2 and 3, an unexpected trend is found: the pressures recorded after the third collapse are now proportional to the initial stand-off distance.
In order to understand this behavior, the maximum vapor volume during rebound and the position of the focal point of the subsequent collapse are investigated. The maximum vapour volume (or the radius of the equivalent bubble) is reversed in order when comparing rebound 1 to rebound 2. Additionally, the focal points are shifted towards the wall in all cases. Both contributions affect the pressure recorded at the wall. However, the maximum pressure at the focal point of a (focussed) cloud collapse is not necessarily related to the pressure recorded at the wall. For example, in case of d = 10 mm, the pressure detected at the focal point is much higher during the primary collapse compared to the subsequent collapses.

Impact on the wall
The "numerical pressure transducer" is used in order to obtain the force acting on the sensor area of 15 x 15 mm². Since this area is relatively large compared to the size of individual bubbles, pressure signals caused by single bubble collapses are blurred over the sensor area. Although this averaging effect of a large sensor area is insufficient to identify the intensity of strong collapses that occur locally, it allows for analysis of the overall temporal evolution of the pressure at the wall. As an example, Figure 6c shows a strong collapse right at the wall, while the pressure transducer (Fig. 5, blue line) records just a medium peak, which is orders in magnitude lower than the value at the focal point. Thus, it is important to analyse simultaneously outputs of various detectors to obtain sufficient insight into the collapse dynamics.
However, the data obtained by the collapse detector can be directly used to investigate collapse focussing.
For setup II the identical initial bubble configuration is applied for all four cases, which allows to compare directly the output of the numerical transducer at the wall. For the same reason the collapse detector would show nearly identical results among different stand-off distances and is not used for comparison. This is, however, strictly true for the primary collapses, only.

Pressure magnitude comparison
The pressure peaks captured by the collapse detector are grid dependent as long as the spatial resolution is not sufficient to capture the smallest relevant physical length scales [3]. Therefore, the output of the collapse detector can be used for a relative comparison of different cases. Mihatsch et al [13] showed that a scaling procedure can be applied to compensate resolution dependencies but in order to perform this scaling, experimental data is required. Due to the lack of experimental validation possibilities for the cases investigated in this study, a scaling of the data obtained by the collapse detector is not yet possible. Therefore, the pressure magnitudes reported in this paper, when obtained with the collapse detector, are a measure for relative comparison between cases and setups but not necessarily absolute values.

Pressure initialization process
For this work, the initial pressure distribution was simplified in the following way. The pressure inside the bubbles is the saturation pressure, while the pressure within the whole domain filled with liquid was initialized with the value of 100 bar. This setup enables a straight forward reproducibility for other research groups who might be interested in reproducing the cases presented in this paper. However, from a physical point of view, it would be more realistic to assume that the Laplacian of the initial pressure field is zero but without featuring a jump in pressure at the surfaces of the bubbles. Such initialization is described by Schmidt et al [3] and will be used for upcoming investigations as well.

Summary and conclusions
This study shows the effects of the cloud interaction parameter and of the stand-off distance on the collapse of a cloud of vapour bubbles in a liquid ambient. Maximum loads at the surface are recorded by application of a numerical force/pressure transducer. Collapse focussing is analysed by application of a procedure which detects isolated collapses (collapse detection algorithm by Mihatsch at al [13]). Two series of numerical simulations with in total seven cases were performed and analysed.
The major findings are: 1) The cloud interaction parameter proposed by Brennen [4] is a suitable measure for bubble interactions and shock focussing. From this investigation it can be concluded that an interaction parameter of O(1) or smaller results in weak interaction without shock/collapse focussing. Once the interaction parameter is O(10) or higher, strong interaction and pronounced focussing are to be expected. For partial cavities [17] a large interaction parameter can be estimated by visual inspection of clouds observed in experiments. Thus, collapse focussing might be a very relevant physical process in cavitating flows occurring in hydraulic machinery. Increasing the interaction parameter β leads to an enhancement of geometrical focusing and of the collapse pressure of individual bubbles. Additionally, the intensification of vapor production close to the wall during rebound is observed and, consequently, an increase of maximum loads at the wall is detected.
2) The stand-off distance is well suited to quantify the intensity of the primary cloud collapse with respect to its erosive potential at an adjacent material surface. However, this investigation shows that secondary and even third collapses subsequent to rebound processes can cause intense loads at material surfaces, too. Therefore, primary collapses and collapses of rebounded vapor structures are both potentially relevant for erosion. The underlying physical mechanism is purely inertia controlled. The primary collapse induces a velocity field which targets towards the wall and rebounding vapour patters are advected with this velocity field. Secondary and tertiary collapses occur closer and closer to the wall and the resulting load at the wall increases as well.
Further studies of this work will include nonhomogeneous initial bubble radii in order to assess effects of micro-scale details aside from macro-scale quantities such as the interaction parameter. Further numerical developments will be the inclusions of viscous effects and non-condensable gas content. The investigation of collapse processes using non-Newtonian liquids is currently under development.

Acknowledgment
This study is part of a project that has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement № 642536. .