Numerical Simulation of Cavitation Erosion Aggressiveness Induced by Unsteady Cloud Cavitation

A numerical investigation of the erosion aggressiveness of leading edge unsteady cloud cavitation based on the energy balance approach has been carried out to ascertain the main damaging mechanisms and the influence of the free stream flow velocity. A systematic approach has permitted the determination of the influence of several parameters on the spatial and temporal distribution of the erosion results comprising the selection of the cavitation model and the collapse driving pressure. In particular, the Zwart, Sauer and Kunz cavitation models have been compared as well as the use of instantaneous versus average pressure values. The numerical results have been compared against a series of experimental results obtained from pitting tests on copper and stainless steel specimens. Several cavitation erosion indicators have been defined and their accuracy to predict the experimental observations has been assessed and confirmed when using a material-dependent damaging threshold level. In summary, the use of the average pressure levels during a sufficient number of simulated shedding cycles combined with the Sauer cavitation model are the recommended parameters to achieve reliable results that reproduce the main erosion mechanisms found in cloud cavitation. Moreover, the proposed erosion indicators follow a power law as a function of the free stream flow velocity with exponents ranging from 3 to 5 depending on their definition.


Introduction
Cavitation is an unique phenomenon in the field of hydrodynamics that occurs when the local pressure in a liquid drops below a critical value, usually close to the vapor pressure, and results in the development of various types of vapor structures such as attached cavities, travelling bubbles, vortical cavities and bubble clouds [1]. Cavitation can typically take place in some widely-used hydraulic machines like pumps, turbines and naval propellers. As a matter of fact, cavitation is often associated with some unwanted consequences like machine performance deterioration, cavitation noise, cavitation-induced vibration and cavitation erosion of solid surfaces [2].
Among the problems caused by cavitation, erosion is the one of the most complex ones since it involves the interaction between fluids and structures. Actually, cavitation erosion is caused by the collapse of the cavities. It has been observed that the collapse of a bubble is a condensation process that ends with the compression of the vapor and the subsequent emission of a shock wave creating a pressure pulse with a very strong amplitude. If the collapse takes place close to a solid wall, a high-speed liquid microjet forms, crosses the bubble and impacts the wall resulting in a very high impulsive pressure. If the impulsive forces resulting either from the impact of the microjet or the length, l, the cavitation number, σ, the shedding frequency, f, and the Strouhal number, St. σ and St have been calculated with Equations (1) and (2), respectively, where P in and P v are the inlet and the vapor saturation pressures, respectively. The computational fluid domain corresponding to the cavitation tunnel test section is shown in Figure 1 and a more detailed description of experiment can be found in Escaler et al. [29] and Couty [30].  (1) and (2), respectively, where Pin and Pv are the inlet and the vapor saturation pressures, respectively. The computational fluid domain corresponding to the cavitation tunnel test section is shown in Figure 1 and a more detailed description of experiment can be found in Escaler et al. [29] and Couty [30].    (2) In the experiment, material specimens 30 mm in width made of standard copper and stainless steel polished down to mirror were mounted along the hydrofoil suction side at different chord positions. The specimens were obtained by an accurate electroerosion machining, glued with cyanoacrylate adhesive and removed after the tests by heating to 400 °C. The erosion intensity was assessed with a statistical analysis of the pitting data measured on the specimens during the incubating period. The pitting results were quantified with the mean pitting rate at the specimen location k, τn (k), and the mean volume deformation rate at location k, τv (k), which were respectively defined as: where Nk is the number of cavitation impacts on location k, Lx and Ly are the side lengths of the tested specimen and Tpit is the time duration of the pitting test. Vd (xi, yi) represents the volume of the indentations on the surface at position (xi, yi), which was measured with a 3D profilometer.  In the experiment, material specimens 30 mm in width made of standard copper and stainless steel polished down to mirror were mounted along the hydrofoil suction side at different chord positions. The specimens were obtained by an accurate electroerosion machining, glued with cyanoacrylate adhesive and removed after the tests by heating to 400 • C. The erosion intensity was assessed with a statistical analysis of the pitting data measured on the specimens during the incubating period. The pitting results were quantified with the mean pitting rate at the specimen location k, τ n (k), and the mean volume deformation rate at location k, τ v (k), which were respectively defined as: τ v (k) = Vd(x i , y i ) T pit × L x × L y (4) where N k is the number of cavitation impacts on location k, L x and L y are the side lengths of the tested specimen and T pit is the time duration of the pitting test. Vd (x i , y i ) represents the volume of the indentations on the surface at position (x i , y i ), which was measured with a 3D profilometer. Figure 2 presents the cavitation erosion results obtained for a U inf of 20 m/s for two different materials: copper and stainless steel. Firstly, the influence of the material strength can be seen by comparing the amplitudes of the erosion indicators for copper in Figure 2a with those for stainless steel in Figure 2b. For copper, the maximum values of τ n (k) are of about 0.03 mm -2 s −1 at 40% of the chord and the maximum values of τ v (k) are of about 160 µm 3 mm −2 s −1 at 50% of the chord. Meanwhile, the maximum τ n (k) and τ v (k) values are of about 0.0004 mm −2 s −1 and 1.9 µm 3 mm −2 s −1 , respectively, at 50% of the chord for stainless steel. Therefore, the pitting rate suffered by the copper is almost two orders of magnitude higher than that for stainless steel. Regarding the location of the damages, the most eroded areas for both materials are found approximately in the range from 40 to 50% of the chord, as expected, because they have been submitted to the same cavitation conditions and the cavity closure region was located at 40% of the chord (l/c = 0.4) as indicated in Table 1. For copper, τ n (k) shows the maximum at 40% of the chord but τ v (k) shows the maximum at 50% of the chord. This seems to indicate that the hydrofoil surface suffers a large amount of impacts with small intensity at 40% and a lower number of impacts with a stronger intensity at 50%. A slightly different erosion pattern is found for stainless steel since it presents the maximum τ n (k) in the range from 40% to 50% but the τ v (k) shows a minimum value at 40% and a maximum value only at 50%. Likewise, the pitting results under other operating conditions show a similar trend as for this condition. Therefore, the erosion results presented on Figure 2 will be used as the reference to validate the numerical predictions.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 5 of 23 50% of the chord for stainless steel. Therefore, the pitting rate suffered by the copper is almost two orders of magnitude higher than that for stainless steel. Regarding the location of the damages, the most eroded areas for both materials are found approximately in the range from 40 to 50% of the chord, as expected, because they have been submitted to the same cavitation conditions and the cavity closure region was located at 40% of the chord (l/c = 0.4) as indicated in Table 1. For copper, τn(k) shows the maximum at 40% of the chord but τv(k) shows the maximum at 50% of the chord. This seems to indicate that the hydrofoil surface suffers a large amount of impacts with small intensity at 40% and a lower number of impacts with a stronger intensity at 50%. A slightly different erosion pattern is found for stainless steel since it presents the maximum τn(k) in the range from 40% to 50% but the τv(k) shows a minimum value at 40% and a maximum value only at 50%. Likewise, the pitting results under other operating conditions show a similar trend as for this condition. Therefore, the erosion results presented on Figure 2 will be used as the reference to validate the numerical predictions.
(a) Copper (b) Stainless steel

Numerical Model
The homogeneous mixture assumption has been used to model the two-phase flow. Thus, velocity and pressures are shared by water and vapor phases and the governing Navier-Stokes equations for the mixture are: where u is the velocity, p is the pressure, t is the time, μm is the mixture dynamic viscosity defined by μm=μvαv+ μl (1-αv) and ρm is mixture density calculated by ρm=ρvαv+ ρl (1-αv).
In the present simulation, the Reynolds Averaged Navier-Stokes (RANS) approach was used, in which the instantaneous quantities are decomposed into the mean and the fluctuating components, i.e., the instantaneous velocity ui is equal to: Using Equation (7) to replace the instantaneous quantities of all the variables in Equations (5) and (6), the RANS equations are obtained:

Numerical Model
The homogeneous mixture assumption has been used to model the two-phase flow. Thus, velocity and pressures are shared by water and vapor phases and the governing Navier-Stokes equations for the mixture are: where u is the velocity, p is the pressure, t is the time, µ m is the mixture dynamic viscosity defined by In the present simulation, the Reynolds Averaged Navier-Stokes (RANS) approach was used, in which the instantaneous quantities are decomposed into the mean and the fluctuating components, i.e., the instantaneous velocity u i is equal to: Using Equation (7) to replace the instantaneous quantities of all the variables in Equations (5) and (6), the RANS equations are obtained: Equations (8) and (9) have the same general form as the instantaneous Navier-Stokes equations, with the velocities and other solution variables now representing time-averaged values, and a series of new additional terms such as, ρu i u j , which are known as the Reynolds stresses induced by the turbulence. In order to close this system of equations, there are two approaches: the first one is based on the eddy-viscosity hypothesis, which relates the Reynolds stresses to the mean flow; the second one is based on solving the transport equation for all the components of the Reynolds stress tensor. In our case, the former method has been used, so a turbulent eddy viscosity, µ t , has been introduced to model the effect of the Reynolds stresses, and the final momentum equation to be solved becomes: In addition, a vapor volume fraction transport equation (Equation (11)) is used to consider the mass transfer between vapor and water: where the term . m represents the mass transfer rate between the two phases. For the current simulation, the turbulent viscosity, µ t , in Equation (10) has been calculated with the SST k-ω model according our previous investigation [31], and µ t has been defined as: where a 1 = 0.31, F 2 is a blending function which restricts the limiter to the wall boundary layer, S is an invariant measure of the strain rate and k and ω are the turbulence kinetic energy and frequency, respectively. In addition, f(ρ) is the density correction to reduce the over-predicted turbulent viscosity as proposed by Reboud et al. [32]: where n is the exponential coefficient which should be specified with a value of 10. Besides, the effects of the cavitation models developed by Zwart [33], Sauer [34] and Kunz [35] on the prediction of erosion were investigated. Table 2 lists the mathematical equations describing each of these cavitation models where R B in the Zwart model is the vapor bubble radius with a constant value of 10 -6 m. Meanwhile, in the Sauer model, it is a variable value that is function of the local vapor volume fraction. α nuc is the nucleation site volume fraction with a default value of 0.0005. Then, t ∞ = C/U inf is the mean time scale and C prod and C dest are the empirical coefficients for vaporization and condensation, respectively, which were taken as 50 and 0.01 for the Zwart model, and 100 and 100 for the Kunz model.
The solution strategy adopted for the current simulations was based on our previous work [31] where the influence of the numerical setting was investigated in detail. Consequently, a medium-sized mesh of 29,749 elements with y+ values in the boundary layers ranging from 0.1 to 3 with a mean value of 1.2 was used. A time step of 3·10 −5 s corresponding to a Root Mean Square (RMS) Courant number of 1.6 was set as well as several successive iterations within each physical time step. A very small residual criterion of 10 −8 and a large iterative number were set to march the solution towards convergence in every time step. The pressure-velocity direct coupling method was used to solve the governing equations. The high-resolution scheme was used for the convection terms. The second-order implicit time scheme was used for the transient term. To accelerate convergence, a series of transient simulations were run from previous steady state models. The simulations were carried out in parallel using twelve cores of an Intel ® Core™ i7-8700 K CPU equipped with 32 GB of RAM. Every unsteady simulation was stopped after at least ten cavity shedding cycles were captured, which took around 12 h for each case. In addition, much effort was taken on building the relationship between the cavity structures and their corresponding erosion intensities during the postprocessing of the obtained results. To achieve a very precise prediction of the position of cavity closure region, which is needed to simulate cavitation erosion [24,36], the cavitation number had to be adjusted to match exactly the same cavity length than in the experiment. Therefore, a cavitation number of 1.55 was used in the simulation, which was based on cavitation tunnel inlet pressure, as compared to the cavitation number of 1.58 set in the experiment.

Cavitation Erosion Model
The cavitation erosion model used in the calculations was developed by Fortes-Patella et al. [17,18]. In this model, the pressure waves emitted during the cavitation collapses that reach the solid wall are the main mechanism responsible for erosion damage. The potential energy of a vapor structure in the fluid domain is defined as: where V vap is the volume of the vapor structure and p d is the driving pressure which forces its collapse.
Then an instantaneous potential power, P pot , can be defined with Lagrangian time derivatives as expressed with Equation (15): Because the vapor volume is related to the vapor volume fraction, α v = V vap /V cell , the potential power in each cell or the potential power density, P den , can be written as: Leclerc et al., [36] found that the second term on the right-hand side of Equation (16) is negligible compared to the first term and they assumed that P pot is released instantaneously only when condensation takes place. This implies that only the first term on the right-hand side of Equation (16) contributes to the radiated power and only if the time derivative of α v is positive. Therefore, P den can be simplified to Equation (17): Note that P pot is defined with Lagrangian time derivatives, which can be obtained substituting the sum of the first and second terms in the right-hand side of Equation (11) with the following relation: Besides, the divergence term, ∇ · u, is actually another form of continuity equation (Equation (5)) defined as ∇ ·u = . m 1 ρ v − 1 ρ l [26]. Hence, Equation (18) can be rewritten as: Appl. Sci. 2020, 10, 5184 8 of 22 Assuming that the positive time derivative contributes to erosion and combining Equations (17) and (19), P den can be finally defined for a given cell as: Equation (20) shows that P den is determined by the driving pressure, p d , which can be defined as the pressure in a cell or the averaged pressure according to the references [18,23] and [26,27], respectively, and by the source term, . m, which depends on the cavitation model being used. Therefore, the definition of p d and the selection of the cavitation model should have an influence on the predicted values of P den .
In order to calculate the potential power reaching a certain position j on the surface of the hydrofoil that has been induced by a collapse occurring at a given point source i on the fluid domain at instant t, the method proposed by Sören and van Terwisga [26] has been applied. For that, it has been assumed that P den (i,t) propagates circumferentially with an infinite wave speed in radial direction without energy losses and that it reaches position j at the same time t, as outlined in Figure 3. Note that Ppot is defined with Lagrangian time derivatives, which can be obtained substituting the sum of the first and second terms in the right-hand side of Equation (11) with the following relation: Besides, the divergence term, ∇ • u, is actually another form of continuity equation (Equation (5)) defined as ∇• = ( − ) [26]. Hence, Equation (18) can be rewritten as: Assuming that the positive time derivative contributes to erosion and combining Equations (17) and (19), Pden can be finally defined for a given cell as: Equation (20) shows that Pden is determined by the driving pressure, pd, which can be defined as the pressure in a cell or the averaged pressure according to the references [18,23] and [26,27], respectively, and by the source term, , which depends on the cavitation model being used. Therefore, the definition of pd and the selection of the cavitation model should have an influence on the predicted values of Pden.
In order to calculate the potential power reaching a certain position j on the surface of the hydrofoil that has been induced by a collapse occurring at a given point source i on the fluid domain at instant t, the method proposed by Sören and van Terwisga [26] has been applied. For that, it has been assumed that Pden(i,t) propagates circumferentially with an infinite wave speed in radial direction without energy losses and that it reaches position j at the same time t, as outlined in Figure  3.  Then, the power reaching position j on the surface of the hydrofoil at instant t, P imp (j,t), can be calculated as: where → R is the position vector of the center j of the surface element from the point source i and → n is the normal vector of this surface element. The contribution of all source points to the loaded power on point j at instant t, P load (j,t), can be calculated integrating Equation (21) over the area (surf ) of the whole computational domain as: Finally, to estimate the cavitation erosion risk on the surface element, the total power received at a given element j due to the accumulation of all the collapses occurred during the period of a complete shedding cycle, P tot_load , can be calculated as: where N ref is the number of time steps simulated during a complete shedding cycle and P threshold is the power threshold above which the material is actually eroded by cavitation. Obviously, P threshold could be assumed to be analogous to a particular material resistance property like the yield stress and it must be validated through experimental tests.
More specifically, P tot_load has been calculated for the first layer of elements on the hydrofoil extrados wall along the chord comprising 50 elements with a width of 2 mm each. Consequently, it has been assumed that P tot_load represents the total cavitation erosion power received on the hydrofoil surface during one characteristic shedding cycle.

Results
In the previous section it has been mentioned that the choice of the cavitation model and the definition of p d in Equation (20) might bring some uncertainty to the calculation of P den and that this can lead to different estimates of P load on the hydrofoil wall. Thus, these possible effects have been investigated in detail based on the transient results obtained during a sufficiently long period of time. Figure 4 shows the numerically predicted time evolution of the total vapor volume within the fluid domain over ten shedding cycles obtained with the three different cavitation models. It can be observed that for the Sauer model the total vapor volume fluctuations are repeatable while for the Zwart and the Kunz models their periods and amplitudes are not so constant from cycle to cycle. Nonetheless, the unsteady cavity behavior is actually quite regular for all the models. In addition, the averaged shedding frequencies calculated in the entire period are 140, 139 and 130 Hz for the Zwart, Sauer and Kunz models, respectively. Thus, it is confirmed that all the results are in good agreement with the expected experimental frequency of 132.8 Hz. Consequently, all the models demonstrate a good performance in capturing the cloud cavitation dynamic behavior.

Influence of Driving Pressure Definition
The numerical results from two consecutive shedding cycles have been selected to calculate Pload on the hydrofoil surface with the different cavitation models and the obtained results have been plotted in Figure 5.
The top graphs in Figure 5 have been obtained considering pd as the flow field instantaneous pressure, p(t). Conversely, the bottom graphs have been obtained considering pd as the averaged pressure over ten cycles, ̅ . It can be observed that the space-time distributions of Pload on the suction side are different depending on which value has been used. For example, when taking the results obtained with the Sauer model, a region with high Pload has been predicted from 20% to 65% of the chord using ̅ in the initial time ranges of each cycle from 0.2 to 0.45 T/Tref and from 1.2 to 1.45 T/Tref. It is observed that using p(t) the calculations can only capture a high Pload during a very short instant around 1.45 T/Tref. Moreover, the amplitude of Pload is also different even though the attack occurs in a similar time-space region. Another large amplitude of Pload has been found with the use of both p(t) and ̅ on a location around 40% of the chord from 0.45 to 1.0 T/Tref for the first cycle and from 1.45 to 2.0 T/Tref for the second cycle. However, this amplitude was higher when taking ̅ as the driving pressure. To finish, similar differences have also been found when comparing the results of the other two cavitation models.

Influence of Driving Pressure Definition
The numerical results from two consecutive shedding cycles have been selected to calculate P load on the hydrofoil surface with the different cavitation models and the obtained results have been plotted in Figure 5.
The top graphs in Figure 5 have been obtained considering p d as the flow field instantaneous pressure, p(t). Conversely, the bottom graphs have been obtained considering p d as the averaged pressure over ten cycles, p. It can be observed that the space-time distributions of P load on the suction side are different depending on which value has been used. For example, when taking the results obtained with the Sauer model, a region with high P load has been predicted from 20% to 65% of the chord using p in the initial time ranges of each cycle from 0.2 to 0.45 T/T ref and from 1.2 to 1.45 T/T ref .
It is observed that using p(t) the calculations can only capture a high P load during a very short instant around 1.45 T/T ref . Moreover, the amplitude of P load is also different even though the attack occurs in a similar time-space region. Another large amplitude of P load has been found with the use of both p(t) and p on a location around 40% of the chord from 0.45 to 1.0 T/T ref for the first cycle and from 1.45 to 2.0 T/T ref for the second cycle. However, this amplitude was higher when taking p as the driving pressure. To finish, similar differences have also been found when comparing the results of the other two cavitation models. chord using ̅ in the initial time ranges of each cycle from 0.2 to 0.45 T/Tref and from 1.2 to 1.45 T/Tref. It is observed that using p(t) the calculations can only capture a high Pload during a very short instant around 1.45 T/Tref. Moreover, the amplitude of Pload is also different even though the attack occurs in a similar time-space region. Another large amplitude of Pload has been found with the use of both p(t) and ̅ on a location around 40% of the chord from 0.45 to 1.0 T/Tref for the first cycle and from 1.45 to 2.0 T/Tref for the second cycle. However, this amplitude was higher when taking ̅ as the driving pressure. To finish, similar differences have also been found when comparing the results of the other two cavitation models.  To understand such differences, the obtained . m from vapor to water has been presented in Figure 6 for the second cycle at various relative instants, T/T ref , marked with dotted lines in Figure 5b. Three stages of the cavity shedding process can be roughly identified. Firstly, the formation of the sheet cavity and its detachment occurs from 1.0 to 1.2 at the initial stage of the shedding process. Figure 6a,b at instants 1.0 and 1.1, respectively, show how the main sheet cavity detaches due to the re-entrant jet and then how a cloud cavity forms and begins to be convected downstream. The second stage ranges from 1.2 to 1.45 and corresponds to the collapse of the cloud cavity. In Figure 6c,d, it can be seen at instant 1.3 how the cloud cavity flows downstream and starts to collapse. Afterwards, the final collapse occurs at instant 1.44. Finally, the stage from 1.45 to 2.0 corresponds to the new growth and formation of the attached sheet cavity. As shown in Figure 6e,f, the sheet cavity reaches its maximum length and then a stagnation point at its closure region forms because the flow over the cavity turns towards the surface. This stagnation point creates a high pressure region which drives the upstream re-entrant jet below the attached cavity and triggers the collapse of the small vapor structures detached at the rear part of the sheet. To understand such differences, the obtained from vapor to water has been presented in Figure 6 for the second cycle at various relative instants, T/Tref, marked with dotted lines in Figure 5b. Three stages of the cavity shedding process can be roughly identified. Firstly, the formation of the sheet cavity and its detachment occurs from 1.0 to 1.2 at the initial stage of the shedding process. Figure 6a,b at instants 1.0 and 1.1, respectively, show how the main sheet cavity detaches due to the re-entrant jet and then how a cloud cavity forms and begins to be convected downstream. The second stage ranges from 1.2 to 1.45 and corresponds to the collapse of the cloud cavity. In Figure 6c,d, it can be seen at instant 1.3 how the cloud cavity flows downstream and starts to collapse. Afterwards, the final collapse occurs at instant 1.44. Finally, the stage from 1.45 to 2.0 corresponds to the new growth and formation of the attached sheet cavity. As shown in Figure 6e,f, the sheet cavity reaches its maximum length and then a stagnation point at its closure region forms because the flow over the cavity turns towards the surface. This stagnation point creates a high pressure region which drives the upstream re-entrant jet below the attached cavity and triggers the collapse of the small vapor structures detached at the rear part of the sheet. The results obtained using both p(t) and ̅ have been plotted on the left-and right-hand sides of Figure 7 for comparison. In particular, the corresponding distributions of p(t) and ̅ , the distribution of Pden as well as the evolution of Pload on the hydrofoil surface along chord (red line) have been plotted at different instants. Moreover, the isolines of αv = 0.1 (black lines) have been superimposed to those graphs to show the location of the main vapor cavities. From a general overview, it can be confirmed that the value of pd has a significant influence on the results. The results obtained using both p(t) and p have been plotted on the left-and right-hand sides of Figure 7 for comparison. In particular, the corresponding distributions of p(t) and p, the distribution of P den as well as the evolution of P load on the hydrofoil surface along chord (red line) have been plotted at different instants. Moreover, the isolines of α v = 0.1 (black lines) have been superimposed to those graphs to show the location of the main vapor cavities. From a general overview, it can be confirmed that the value of p d has a significant influence on the results. The results presented in Figure 6a at instant 1.0 show that the highest condensation processes take place at both ends of the attached cavity. In Figure 7, it is also confirmed that the values of p(t) in the condensation regions at this instant are higher than the corresponding values of ̅ . It can be seen that a higher Pden is predicted on the fluid using p(t), and correspondingly it results in a higher Pload on the hydrofoil surface. Likewise, at instant 1.1, the p(t) in the condensation region is also higher than ̅ , leading to a locally higher Pload. However, this is not a reliable result because, at this stage, when the sheet cavity is transformed into a cloud cavity, no strong collapses are expected to occur.
The results presented in Figure 6c at instant 1.3 show condensation regions located on the margin of the cloud cavity because the pressure in its outer region is definitely higher than that within its interior as shown in Figure 7c. Moreover, p(t) is very low and even almost equal to the vapor pressure inside the cloud cavity region marked by the isoline, which implies that the pressure difference between the driving pressure and vapor pressure is close to zero. As a result, a very low aggressiveness on the hydrofoil surface is predicted. On the other hand, the use of ̅ determines a higher Pden and Pload at the condensation region.
The final cloud collapse takes place at instant 1.44, as shown in Figures 6d and 7d, when a local value of p(t) over 1 MPa is predicted inducing the maximum values of Pden and Pload. However, this instantaneous pressure peak is the consequence of the final collapse but it is not the pressure source driving the cavity collapse. Moreover, the appearance of this local pressure peak is considered by Bensow et al., [37] to be a spurious effect of the numerical calculation. In contrast, the Pden and Pload predicted using ̅ are much lower because the ̅ levels are significantly smaller than the p(t) ones.
Finally, the area with the highest erosion power is mainly concentrated on the closure region of the sheet cavity as shown in Figure 7e at instant 1.6. In addition, a small cavity also appears after the final collapse of the cloud cavity due to the vortex rebound [38]. Comparing the results predicted using p(t) or ̅ , the latter gives higher Pden and Pload because its level on the condensation area is higher than for the former. Similarly, the results at instant 1.8 allow us to verify the previous observations. The results presented in Figure 6a at instant 1.0 show that the highest condensation processes take place at both ends of the attached cavity. In Figure 7, it is also confirmed that the values of p(t) in the condensation regions at this instant are higher than the corresponding values of p. It can be seen that a higher P den is predicted on the fluid using p(t), and correspondingly it results in a higher P load on the hydrofoil surface. Likewise, at instant 1.1, the p(t) in the condensation region is also higher than p, leading to a locally higher P load . However, this is not a reliable result because, at this stage, when the sheet cavity is transformed into a cloud cavity, no strong collapses are expected to occur.
The results presented in Figure 6c at instant 1.3 show condensation regions located on the margin of the cloud cavity because the pressure in its outer region is definitely higher than that within its interior as shown in Figure 7c. Moreover, p(t) is very low and even almost equal to the vapor pressure inside the cloud cavity region marked by the isoline, which implies that the pressure difference between the driving pressure and vapor pressure is close to zero. As a result, a very low aggressiveness on the hydrofoil surface is predicted. On the other hand, the use of p determines a higher P den and P load at the condensation region.
The final cloud collapse takes place at instant 1.44, as shown in Figures 6d and 7d, when a local value of p(t) over 1 MPa is predicted inducing the maximum values of P den and P load . However, this instantaneous pressure peak is the consequence of the final collapse but it is not the pressure source driving the cavity collapse. Moreover, the appearance of this local pressure peak is considered by Bensow et al., [37] to be a spurious effect of the numerical calculation. In contrast, the P den and P load predicted using p are much lower because the p levels are significantly smaller than the p(t) ones.
Finally, the area with the highest erosion power is mainly concentrated on the closure region of the sheet cavity as shown in Figure 7e at instant 1.6. In addition, a small cavity also appears after the final collapse of the cloud cavity due to the vortex rebound [38]. Comparing the results predicted using p(t) or p, the latter gives higher P den and P load because its level on the condensation area is higher than for the former. Similarly, the results at instant 1.8 allow us to verify the previous observations. Figure 8 shows the accumulated P tot_load levels along the hydrofoil chord during one shedding cycle from instant 1 to 2 using Equation (23) with P threshold equal to zero. It is observed that P tot_load along the chord is much lower when using p(t) than p, with the exception of a small region near the leading edge when using the Zwart and the Sauer models. Another observation is that the P tot_load levels near the leading edge are even higher than those on the region where the cloud cavity collapses when using p(t). As already discussed before, this is because p(t) gives a higher P den and P load during the transition from sheet to cloud cavity but underestimates P den during the cloud cavity collapse because the pressure difference between p(t) and vapor pressure is close to zero. Therefore, the results given by p(t) are not reasonable because they are in contradiction with the experimental results presented in Figure 2 proving that no significant erosion occurs along the chord in the range from 0 to 20%. Meanwhile, strong erosion is found downstream of the cavity closure region. In summary, it is confirmed that the cavitation erosion prediction results are sensitive to the definition of the driving pressure, and that the results obtained with p seem to be more in agreement with the experiments. Consequently, they will be taken into account in the following sections.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 23 Figure 8 shows the accumulated Ptot_load levels along the hydrofoil chord during one shedding cycle from instant 1 to 2 using Equation (23) with Pthreshold equal to zero. It is observed that Ptot_load along the chord is much lower when using p(t) than ̅ , with the exception of a small region near the leading edge when using the Zwart and the Sauer models. Another observation is that the Ptot_load levels near the leading edge are even higher than those on the region where the cloud cavity collapses when using p(t). As already discussed before, this is because p(t) gives a higher Pden and Pload during the transition from sheet to cloud cavity but underestimates Pden during the cloud cavity collapse because the pressure difference between p(t) and vapor pressure is close to zero. Therefore, the results given by p(t) are not reasonable because they are in contradiction with the experimental results presented in Figure 2 proving that no significant erosion occurs along the chord in the range from 0 to 20%. Meanwhile, strong erosion is found downstream of the cavity closure region. In summary, it is confirmed that the cavitation erosion prediction results are sensitive to the definition of the driving pressure, and that the results obtained with ̅ seem to be more in agreement with the experiments. Consequently, they will be taken into account in the following sections.   Figure 9 demonstrates that the contours of ̅ obtained with the different cavitation models are similar, which reinforces the assumption that the study of the influence of the cavitation model on Pden and Pload can be done using ̅ as driving pressure.     Figure 9 demonstrates that the contours of p obtained with the different cavitation models are similar, which reinforces the assumption that the study of the influence of the cavitation model on P den and P load can be done using p as driving pressure.

Influence of Cavitation Model
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 23 Figure 8 shows the accumulated Ptot_load levels along the hydrofoil chord during one shedding cycle from instant 1 to 2 using Equation (23) with Pthreshold equal to zero. It is observed that Ptot_load along the chord is much lower when using p(t) than ̅ , with the exception of a small region near the leading edge when using the Zwart and the Sauer models. Another observation is that the Ptot_load levels near the leading edge are even higher than those on the region where the cloud cavity collapses when using p(t). As already discussed before, this is because p(t) gives a higher Pden and Pload during the transition from sheet to cloud cavity but underestimates Pden during the cloud cavity collapse because the pressure difference between p(t) and vapor pressure is close to zero. Therefore, the results given by p(t) are not reasonable because they are in contradiction with the experimental results presented in Figure 2 proving that no significant erosion occurs along the chord in the range from 0 to 20%. Meanwhile, strong erosion is found downstream of the cavity closure region. In summary, it is confirmed that the cavitation erosion prediction results are sensitive to the definition of the driving pressure, and that the results obtained with ̅ seem to be more in agreement with the experiments. Consequently, they will be taken into account in the following sections.   Figure 9 demonstrates that the contours of ̅ obtained with the different cavitation models are similar, which reinforces the assumption that the study of the influence of the cavitation model on Pden and Pload can be done using ̅ as driving pressure.        In Figure 10a, the Zwart model predicts a cloud cavity at instant 1.3 with low condensation rate on its outer region and a maximum Pload of about 60 kWm −2 which is the lowest figure compared to the other models. On the contrary, the Sauer and the Kunz models predict a smaller cavity volume and a higher Pload induced by the higher condensation rate on the outer region of the cloud. More specifically, the Kunz model predicts the highest Pload of about 240 kW m −2 at 40% of the chord because a part of the condensation region of the cloud cavity is much closer to the hydrofoil surface than in the rest of results. Similarly, the Zwart model predicts again the smallest condensation rate and Pload at instant 1.4. The other two models predict a relatively higher intensity especially for the Kunz one. At instant 1.5, the cloud cavity still presents a large volume with the Zwart model, meanwhile with the Sauer model the cavity already vanishes at instant 1.44, as shown in Figure 6d, and it reappears with a very small volume due to the collapse rebound. With the Kunz model, the cloud cavity finishes its final collapse at this particular instant. This indicates that the final collapse occurs at different instants and at different chord locations depending on the cavitation model. Furthermore, the collapse process predicted with the Zwart model creates a lower Pload than the rest. Comparing Figures  6d and 10c, it can be seen that the Sauer model predicts a much higher condensation rate than with the Kunz model when the final collapse occurs, which leads to a higher Pden and Pload, as shown in Figure 7d. Figure 11 shows the results at instant 0.5 corresponding to an analogous situation to the one observed at instant 1.5 presented in Figure 10c. By comparing these results, it can be seen that the cavity behavior differs between the two consequent cycles depending on the cavitation model being used. The cavity topologies at these two instants are very similar with the Sauer model, showing that the sheet cavity has reached its maximum length with a high Pload at its closure region and that a very small cloud cavity appears at 70% of the chord inducing a small Pload. Likewise, the cavity behavior and the Pload provided by the Zwart model at these two instants are also similar. However, the results for the Kunz model are different since there is a cloud with a larger volume of vapor that induces the In Figure 10a, the Zwart model predicts a cloud cavity at instant 1.3 with low condensation rate on its outer region and a maximum P load of about 60 kWm −2 which is the lowest figure compared to the other models. On the contrary, the Sauer and the Kunz models predict a smaller cavity volume and a higher P load induced by the higher condensation rate on the outer region of the cloud. More specifically, the Kunz model predicts the highest P load of about 240 kW m −2 at 40% of the chord because a part of the condensation region of the cloud cavity is much closer to the hydrofoil surface than in the rest of results. Similarly, the Zwart model predicts again the smallest condensation rate and P load at instant 1.4. The other two models predict a relatively higher intensity especially for the Kunz one. At instant 1.5, the cloud cavity still presents a large volume with the Zwart model, meanwhile with the Sauer model the cavity already vanishes at instant 1.44, as shown in Figure 6d, and it reappears with a very small volume due to the collapse rebound. With the Kunz model, the cloud cavity finishes its final collapse at this particular instant. This indicates that the final collapse occurs at different instants and at different chord locations depending on the cavitation model. Furthermore, the collapse process predicted with the Zwart model creates a lower P load than the rest. Comparing Figures 6d and 10c, it can be seen that the Sauer model predicts a much higher condensation rate than with the Kunz model when the final collapse occurs, which leads to a higher P den and P load , as shown in Figure 7d. Figure 11 shows the results at instant 0.5 corresponding to an analogous situation to the one observed at instant 1.5 presented in Figure 10c. By comparing these results, it can be seen that the cavity behavior differs between the two consequent cycles depending on the cavitation model being used. The cavity topologies at these two instants are very similar with the Sauer model, showing that the sheet cavity has reached its maximum length with a high P load at its closure region and that a very small cloud cavity appears at 70% of the chord inducing a small P load . Likewise, the cavity behavior and the P load provided by the Zwart model at these two instants are also similar. However, the results for the Kunz model are different since there is a cloud with a larger volume of vapor that induces the maximum P load at 75% of the chord at instant 0.5, which cannot be observed at instant 1.5 in Figure 10c. Similarly, the space-time distributions of P load in Figure 5 are significantly different between the time span from 0 to 1 and from 1 to 2 for the Kunz model, meanwhile they are quite similar for the other two models. Consequently, it can concluded that the Zwart and Sauer models provide similar predictions of the erosion power in terms of erosion aggressiveness and location from cycle to cycle, while the Kunz model provides results that differ from cycle to cycle.   Figure 12 represents compares the values of Ptot_load without setting any threshold value for two consecutive shedding cycles obtained with the three cavitation models. With the Zwart model, a similar shape and amplitude of the Ptot_load distribution is found in both cycles although the second one gives a slightly higher intensity. With the Sauer model, the two cycles show exactly the same result. In contrast, with the Kunz model a significant difference is found in the second half of the chord. More specifically, Ptot_load drops with a faster rate towards a very small amplitude at 80% of the chord during the second cycle, meanwhile Ptot_load decreases more gradually and it is still significant even at the trailing edge during the first cycle. Note that similar differences between two consecutive cycles are also observed along the 10 cycles shown in Figure 5 for all the cavitation models. Another conclusion that can be extracted from Figure 12 is that the cavitation model also influences the erosion distribution along the chord. More specifically, all three models give the maximum values of Ptot_load around 40% of the chord at the attached cavity closure region but they give different results at the cloud cavity collapse region. With the Zwart model, the collapses of the clouds, shown on the left hand side of Figure 10, are less intense and they take place close to the trailing edge presenting a low amplitude of Ptot_load from 50 to 100% of the chord. With the Sauer model,  Figure 12 represents compares the values of P tot_load without setting any threshold value for two consecutive shedding cycles obtained with the three cavitation models. With the Zwart model, a similar shape and amplitude of the P tot_load distribution is found in both cycles although the second one gives a slightly higher intensity. With the Sauer model, the two cycles show exactly the same result. In contrast, with the Kunz model a significant difference is found in the second half of the chord. More specifically, P tot_load drops with a faster rate towards a very small amplitude at 80% of the chord during the second cycle, meanwhile P tot_load decreases more gradually and it is still significant even at the trailing edge during the first cycle. Note that similar differences between two consecutive cycles are also observed along the 10 cycles shown in Figure 5 for all the cavitation models.   Figure 12 represents compares the values of Ptot_load without setting any threshold value for two consecutive shedding cycles obtained with the three cavitation models. With the Zwart model, a similar shape and amplitude of the Ptot_load distribution is found in both cycles although the second one gives a slightly higher intensity. With the Sauer model, the two cycles show exactly the same result. In contrast, with the Kunz model a significant difference is found in the second half of the chord. More specifically, Ptot_load drops with a faster rate towards a very small amplitude at 80% of the chord during the second cycle, meanwhile Ptot_load decreases more gradually and it is still significant even at the trailing edge during the first cycle. Note that similar differences between two consecutive cycles are also observed along the 10 cycles shown in Figure 5 for all the cavitation models. Another conclusion that can be extracted from Figure 12 is that the cavitation model also influences the erosion distribution along the chord. More specifically, all three models give the maximum values of Ptot_load around 40% of the chord at the attached cavity closure region but they give different results at the cloud cavity collapse region. With the Zwart model, the collapses of the clouds, shown on the left hand side of Figure 10, are less intense and they take place close to the trailing edge presenting a low amplitude of Ptot_load from 50 to 100% of the chord. With the Sauer model, Another conclusion that can be extracted from Figure 12 is that the cavitation model also influences the erosion distribution along the chord. More specifically, all three models give the maximum values of P tot_load around 40% of the chord at the attached cavity closure region but they give different results at the cloud cavity collapse region. With the Zwart model, the collapses of the clouds, shown on the left hand side of Figure 10, are less intense and they take place close to the trailing edge presenting a low amplitude of P tot_load from 50 to 100% of the chord. With the Sauer model, the cloud collapses are concentrated from 50 to 60% of the chord and the results show a small peak around 55% of the chord. With the Kunz model, the results are sensitive to the shedding cycle being considered. According to the experimental results presented in Figure 2, the main eroded area was located from 30 to the 70% of the chord and the maximum erosion was found from 40 to 50% of the chord. Based on these results, it can be concluded that the cavitation power distribution predicted by the Sauer model agrees more precisely with the experimental observations.

Influence of Cavitation Model
Based on the previous discussion, it can be stated that different cavitation models give different results of cavitation aggressiveness regarding the condensation . m, which determines the level of P den , and regarding the cavity topology, which determines the efficiency of the transmission from P den to P load on the surface and the size and location of the area subjected to the highest P tot_load . Nonetheless, it can be concluded that the Sauer model giver a more accurate prediction of cavitation aggressiveness than the rest.

The Mechanism of Cavitation Erosion
In this section, the numerical results obtained with the Sauer model have been used to discuss the mechanisms of cavitation erosion based on the comparison with the experimental observations. For that, the second cycle results from instants 1 to 2 have been considered without any loss of generality because the Sauer model has shown a good repeatability for all the simulated shedding cycles. Figure 12b shows that the maximum P tot_load is found both at the closure region of the attached sheet cavity and at the location where the main cloud cavity collapses. Nevertheless, the calculation of P tot_load without taking into account the fact that a P threshold level exists is not an accurate method to estimate the actual erosion risk. Therefore, a P threshold has been set to eliminate the contribution of low intensity collapses not being sufficiently high to cause material damage. For that, the number of effective time steps at which P load is higher than P threshold at any chord location, N eff , has been considered to estimate the average erosive power load, P ave_load , at each effective time step and hydrofoil position with the following expression: If N ref is the number of time steps in one cycle of duration T ref , then the ratio N eff /N ref accounts for the percent time duration during one shedding cycle with P load values higher than P threshold .
In the present study, two P threshold levels with values of 30 kW m −2 and 90 kW m −2 have been considered to investigate the influence of the particular material resistance to cavitation erosion and the corresponding distributions of P tot_load , N eff /N ref and P ave_load along the chord have been plotted in Figure 13. Note here that it can be assumed that the lower level of 30 kW m −2 represents copper and that the level of 90 kW m −2 represents stainless steel since there is no doubt that copper has a lower yield strength than stainless steel.
For P threshold of 30 kW m −2 , the maximum values of P tot_load and N eff /N ref are found at 40% of the chord but at this location P ave_load is relatively small compared to the rest of values farther downstream. This is because small vapor structures are collapsing at the cavity closure region for a long time duration (0.7 T ref ) but they are individually inducing relatively small P load . This prediction is in agreement with the experimental results that show a high pitting rate at this location but with a relatively small deformation volume rate. The second peak value of P tot_load is located at 53% of the chord. In this case, N eff /N ref is around 0.2 while P ave_load shows its maximum value. This indicates that, at the region of the cloud collapses, the erosion intensity is very high but its time duration during one cycle is relatively short. This is again in agreement with the experiment that shows a relatively smaller pitting rate but with the highest deformation volume rate in the same position.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 23 place during a long time at the sheet cavity closure is dampened but it does not affect the results at the cloud collapse region where the Ptot_load is very high and clearly above Pthreshold. Therefore, the setting of Pthreshold is necessary to take into account the material resistance in the numerical prediction of cavitation erosion. Moreover, this is a concept that can also be used to explain the experimental observations found with different types of materials if it is assumed that it reproduces the effects of the material yield stress.
For example, the numerical Ptot_load and the experimental τn and τv have been normalized by their respective maxima along the chord and compared in Figure 14. It can be seen that the predicted results for Pthreshold of 30 kW m −2 compare well with the copper erosion measurements, and the ones for 90kW m −2 compare well with stainless steel especially for the pitting rate observations. In addition, the numerical results show a zone with a low erosion intensity like a groove between the two maximum erosion picks for Ptot_load, which is due to the fact that in the simulations the boundary between the vapor and water phases is very well delimited whereas in reality these boundaries are quite unstable and the cavity wake is much fuzzier due to the flow turbulence [21]. Given the good results obtained with the present simulations, the Sauer model has been further used to investigate the effect of the free stream velocity on the cavitation erosion power. When P threshold increases to 90 kW m −2 , the strongest erosions are also predicted at the sheet cavity closure and at the location where the main clouds collapse. However, it can be seen in Figure 13a,b that P tot_load and N eff /N ref decrease significantly in the sheet cavity closure region but only slightly in the cloud collapse region when compared with the results obtained with 30 kW m −2 . Now the maximum P tot_load is located at around 53% of the chord instead of 40%, which is also in accordance with the experimental results for the stainless steel specimen that showed the highest pitting rate at around the 50%. These numerical predictions confirm the existence of a threshold level and the fact that it permits the correction of the estimation of the cavitation erosion risk based on the cavitation aggressiveness. Increasing P threshold , the contribution to erosion of the small intensity collapses taking place during a long time at the sheet cavity closure is dampened but it does not affect the results at the cloud collapse region where the P tot_load is very high and clearly above P threshold . Therefore, the setting of P threshold is necessary to take into account the material resistance in the numerical prediction of cavitation erosion. Moreover, this is a concept that can also be used to explain the experimental observations found with different types of materials if it is assumed that it reproduces the effects of the material yield stress.
For example, the numerical P tot_load and the experimental τ n and τ v have been normalized by their respective maxima along the chord and compared in Figure 14. It can be seen that the predicted results for P threshold of 30 kW m −2 compare well with the copper erosion measurements, and the ones for 90kW m −2 compare well with stainless steel especially for the pitting rate observations. In addition, the numerical results show a zone with a low erosion intensity like a groove between the two maximum erosion picks for P tot_load , which is due to the fact that in the simulations the boundary between the vapor and water phases is very well delimited whereas in reality these boundaries are quite unstable and the cavity wake is much fuzzier due to the flow turbulence [21].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 18 of 23 place during a long time at the sheet cavity closure is dampened but it does not affect the results at the cloud collapse region where the Ptot_load is very high and clearly above Pthreshold. Therefore, the setting of Pthreshold is necessary to take into account the material resistance in the numerical prediction of cavitation erosion. Moreover, this is a concept that can also be used to explain the experimental observations found with different types of materials if it is assumed that it reproduces the effects of the material yield stress.
For example, the numerical Ptot_load and the experimental τn and τv have been normalized by their respective maxima along the chord and compared in Figure 14. It can be seen that the predicted results for Pthreshold of 30 kW m −2 compare well with the copper erosion measurements, and the ones for 90kW m −2 compare well with stainless steel especially for the pitting rate observations. In addition, the numerical results show a zone with a low erosion intensity like a groove between the two maximum erosion picks for Ptot_load, which is due to the fact that in the simulations the boundary between the vapor and water phases is very well delimited whereas in reality these boundaries are quite unstable and the cavity wake is much fuzzier due to the flow turbulence [21]. Given the good results obtained with the present simulations, the Sauer model has been further used to investigate the effect of the free stream velocity on the cavitation erosion power. Given the good results obtained with the present simulations, the Sauer model has been further used to investigate the effect of the free stream velocity on the cavitation erosion power.

Free Stream Velocity Effects
Four different free stream water velocities, U inf , have been simulated with the same numerical model corresponding to 15, 20, 25 and 30 m/s. The time step has been set inversely proportional to U inf in order to keep the Courant number constant between all the simulations. In addition, the outlet pressure has been adjusted correspondingly in order to get always the same cavity length around the 40% of the chord as in the experiment. The model setup main parameters and the shedding frequencies calculated numerically and measured experimentally have been listed in Table 3. It can be seen that the simulations predict with good accurately the shedding frequencies for any U inf with a maximum percent deviation of 10.3%. Similarly to the results shown in Figure 4b, the time signals of the total vapor volume are repeatable from cycle to cycle for all the cases. Finally, the driving pressure has been based on the averaged instantaneous pressures over ten shedding cycles. To study the relationship between cavitation aggressiveness and flow velocity, the collapse efficiency proposed by Fortes-Patella et al., [39], η collapse , has been used which quantifies the effective energy transfer between the potential power of the vapor volume and the actually erosive power. The value of η collapse is dependent on the initial gas pressure inside the bubble, P g0 , and the environment pressure, P ∞ , but it is weakly dependent on the initial vapor volume so it can be calculated as: η collapse = 0.029 P g0 /P ∞ −0.54 (25) Equation (25) shows that, for a given P g0 , the efficiency is higher the higher is P ∞ which means that for a given cavitation number, the collapses will release more energy for higher U inf . In the present work, P g0 has been taken as 1500 Pa and P ∞ has been calculated with the numerical cavitation number as done in [39]. Therefore, η collapse is a constant coefficient depending on U inf as shown in Table 4. Then, the actual effective erosion power load, P eff_load , can be calculated as: Since the experimentally obtained Strouhal number at different U inf is constant with a value around 0.28 (see Table 1), the shedding frequency will increase linearly with U inf and the collapsing frequency on the hydrofoil will also increase. Therefore, the cavitation erosion aggressiveness per unit time, P agg , can defined as: P agg = P e f f _load /T = P e f f _load · f (27) All the calculated values of the various cavitation erosion indicators defined in this section have been listed in Table 4 and plotted in Figure 15, where several power law relations between them and U inf have been found. Note that P tot_load and P eff_load are the accumulated potential power and the accumulated effective power for one complete shedding cycle, respectively, and P agg is the accumulated effective power per unit time. The maximum P tot_load along the chord, which is found at 40% of the chord for each U inf because of no threshold has been considered, has also been selected to quantify the flow velocity effects.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 20 of 23 All the calculated values of the various cavitation erosion indicators defined in this section have been listed in Table 4 and plotted in Figure 15, where several power law relations between them and Uinf have been found. Note that Ptot_load and Peff_load are the accumulated potential power and the accumulated effective power for one complete shedding cycle, respectively, and Pagg is the accumulated effective power per unit time. The maximum Ptot_load along the chord, which is found at 40% of the chord for each Uinf because of no threshold has been considered, has also been selected to quantify the flow velocity effects. As observed in Figure 15a, Ptot_load increases approximately as Uinf to the third power for one shedding cycle, which is in accordance with the numerical results obtained by other researchers like Carrat et al., [23], Leclercq et al., [24] and Melissaris et al., [27], although they did not take into account the collapse efficiency. As observed in Figure 15b, Peff_load increases approximately as Uinf to the fourth power when considering the efficiency which is in accordance with the numerical results obtained by Fortes-Patella et al. [39]. Finally, because the shedding frequency increases with the flow velocity, Pagg increases approximately as Uinf to the fifth power as shown in Figure 15c but this curve might rise faster, i.e., present a higher exponent, because the shedding frequency has been slightly overestimated for the lowest Uinf and underestimated for the highest Uinf in the present calculations.  As observed in Figure 15a, P tot_load increases approximately as U inf to the third power for one shedding cycle, which is in accordance with the numerical results obtained by other researchers like Carrat et al., [23], Leclercq et al., [24] and Melissaris et al., [27], although they did not take into account the collapse efficiency. As observed in Figure 15b, P eff_load increases approximately as U inf to the fourth power when considering the efficiency which is in accordance with the numerical results obtained by Fortes-Patella et al. [39]. Finally, because the shedding frequency increases with the flow velocity, P agg increases approximately as U inf to the fifth power as shown in Figure 15c but this curve might rise faster, i.e., present a higher exponent, because the shedding frequency has been slightly overestimated for the lowest U inf and underestimated for the highest U inf in the present calculations. In addition, the experimental cavitation numbers are higher for the higher velocities than those in the numerical model, which implies higher environment pressures and higher collapse efficiencies for higher velocities. Consequently, the erosion rate may increase with an exponent higher than 5, which would be in agreement with the experimental investigations by Dular et al., [40,41] who found that the erosion damage followed a power law with n values from 5 to 8.

Conclusions
In the present study, the dynamic behavior and the erosion power of unsteady cloud cavitation on a 2D hydrofoil has been investigated numerically based on an energy balance approach. The influences of the driving pressure and of the cavitation model have been discussed in detail. The numerical prediction of erosion is in agreement with the experimental measurements and with the main erosion mechanisms. Moreover, the influence of flow velocity on erosion power has been quantified. As a result, the following conclusions have been obtained: 1.
The selection of the driving pressure to estimate the power of the cavity collapse has a significant effect on the space-time distribution of the cavitation aggressiveness on the hydrofoil surface. The use of the average pressure gives more similar results to the experiment than the use of the instantaneous pressure. 2.
The cavitation model influences significantly the power loaded on the hydrofoil surface both in terms of magnitude and spatial distribution along the chord. For the cases considered in the present study, the Sauer model performs better than the Kunz and Zwart ones.

3.
Two main erosion mechanisms have been predicted that are in good agreement with experimental observations. One is induced by the high frequency of low-intensity collapses taking place at the closure region of the main sheet cavity attached to the hydrofoil surface. The other one is induced by the low frequency and high intensity collapses of the shed cloud cavities.

4.
Power laws have been obtained that permit the calculation of the erosive cavitation intensity as a function of the flow velocity by taking into account the collapse efficiency and the shedding frequency. More specifically, the effective power load law grows with an exponent of 4, and the erosion aggressiveness per unit time grows with an exponent of 5.