On the accuracy of predicting cavitation impact loads on marine propellers

Abstract Predicting the cavitation impact loads on a propeller surface using numerical tools is becoming essential, as the demand for more efficient designs, stretched to the limit, is increasing. One of the possible design limits is governed by cavitation erosion. The accuracy of estimating such loads, using a URANS approach, has been investigated. We follow the energy balance approach by (Schenke and van Terwisga, 2019), (Schenke et al., 2019), where we take account of the focusing of the potential energy into the collapse center before it is radiated as shock wave energy in the domain. In complex flows, satisfying the total energy balance, when reconstructing the radiated energy, has always been an issue in the past. Therefore, in this study, we investigate different considerations for the vapor reduction rate, in order to minimize the numerical errors, when estimating the local surface impact power. We show that when the vapor volume reduction rate is estimated using the mass transfer source term, then all the energy is conserved and the total energy balance is satisfied. The model is verified on a single cavitating bubble collapse, and it is further validated on a model propeller test case. The obtained surface impact distribution agrees well with the experimental paint test results, illustrating the potential for practical use of our fully conservative method to predict cavitation implosion loads on propeller blades.


Introduction
Cavitation is the formation of vapor bubbles of a flowing liquid in regions where the liquid accelerates such that the local pressure of the liquid drops below its vapor pressure. The collapse of these cavities in pressure recovery regions, often leads to vibration and damage of mechanical components, for instance bearings, fuel injectors, impellers, pumps, propellers and rudders. Propeller cavitation might have a big impact on the whole operation of a ship and its environment. Hull vibrations, thrust breakdown, propeller erosion, and underwater radiated noise are the main threats. Cavitation erosion when experienced, normally leads to significant repair and maintenance costs, and even safety issues. Thus, there is a clear need for more efficient designs. Stringent regulations on the fuel efficiency of ships and their propulsion units, make things even more urgent, and designs should be stretched to the limits. For a propeller design this means taking away the margins against unwanted cavitation phenomena as much as possible in favor of higher propulsive efficiency. Most often the performance of propulsion systems is sub-optimal because countermeasures are needed to prevent erosion. Hence, the challenge for modern designs is to find the right balance between efficiency and the allowable level of cavitation. Evaluation of different designs requires objective numerical tools * Correspondence to: Wärtsilä Netherlands BV, Drunen 5151DM, The Netherlands. E-mail address: themis.melissaris@wartsila.com (T. Melissaris). and methods, capable of providing accurate and reliable prediction of cavitation implosion loads, to secure the selection of the best design.
The history of cavitation research dates back to the mid-eighteenth century, when Euler discussed the possibility of a phenomenon, which influences the performance of a water wheel. The first to introduce the phenomenon as it is known today, was Reynolds in 1850's, by discussing the effect it had on the performance of a screw propelled steamer. However, the first observer of cavitation was Parsons, whose ship, Turbinia, suffered from severe thrust breakdown due to cavitation, in 1894. Since then, a lot of knowledge has been gained on cavitation dynamics on propellers, however the underlying physics behind the mechanisms of cavitation erosion due to the implosion of cavitating structures in the vicinity of the propeller surface, remain yet unclear. And even though cavitation dynamics and cavitation behavior on marine propellers are being extensively investigated, for more than a century experimentally [1], and over 20 years numerically, cavitation erosion prediction is still a major challenge.
The assessment of cavitation erosion risk on marine propellers from numerical simulations had not been studied thoroughly, until recently. The lack of physical understanding on how the cavitation implosion loads correlate with the material properties, increases the complexity T. Melissaris   small scales (from mm down to μm and from μs to ns [2]). These are cardinal reasons why research was first focused mainly on the hydrodynamic aspects of cavitation erosion, and primarily on the dynamics of large scale cavitating structures. In one of the first attempts to estimate the cavitation erosion risk numerically on marine propellers, Hasuike et al. [3] investigated the risk of cavitation erosion in four differently loaded propellers operating in a wake flow. They used the erosion indices proposed by Nohmi et al. [4] to estimate the aggressiveness of the cavitation impact loads. However, it is not very clear how these indices are derived and they also seem quite empirical. Ponkratov and Caldas [5] and Ponkratov [6], tried to predict the cavitation erosion risk on a ship scale rudder and propeller, respectively. To estimate the cavitation erosion aggressiveness, they used several numerical functions, developed by the Lloyd's Register Technical Investigation Department (LR TID). Unfortunately, the formulation of these functions was not reported. Usta et al. [7] and Usta and Korkut [8] estimated the erosion aggressiveness using different indicators, as found in the literature. All the indicators are based on the potential energy hypothesis [9,10], which states that the initial potential energy of a cavitating structure is proportional to its volume and the pressure difference driving the collapse. One of the indicators they used was the Erosion Intensity Function by Li et al. [11]. However, in all the indicators, an artificial threshold needs to be defined for the erosive cavitation impact loads. All of the aforementioned attempts to estimate the cavitation erosion risk on the propeller blades are more or less based on the potential energy hypothesis. None of the attempts was able to ensure energy conservation, while they all ignore the spatial and temporal focusing of the potential energy, which actually takes place during a cavity collapse.
Another philosophy applies the microjet erosion model [12,13], as elaborated by Peters et al. [14], considering only collapses of single bubbles near the blade surface. In a previous work [15], we have discussed in detail the importance of the water hammer originating from a microjet impingement, and the pressure wave originated from the collapse of a cavity, with respect to erosion. We hypothesized that the impact load due to a collective cloud collapse, may be much more powerful than the one from a microjet, formed on a single bubble collapsing close to the surface. A more detailed comparison, can be found in Joshi et al. [16], who simulated a single bubble collapse close to a wall, using SPH. They focused on the relation between the shock at the moment of the jet impact, and the shock wave from the eventual collapse of the bubble. The authors also modeled the material plastic deformation. They concluded that, although the water hammer can produce twice the maximum plastic deformation compared to a shock wave from the collapse, the volume of material that is plastically deformed, is miniscule. Basically, a collapse shock wave can plastify almost 800 times larger volume, and hence, the erosion rate will be higher. Considering that both the collapse shock and the microjet impact are eventually fed by the initial potential energy content, and the fact that reliable prediction of a microjet formation and its water hammer would require extremely high resolution in time and space, an energy balance consideration, based on the initial potential energy contained in cavitating structures, is believed to be more successful on a macroscopic scale.
In the present study, an energy conservative method is used to predict the cavitation erosion aggressiveness on a surface. This method is based on the potential energy hypothesis, but it allows for the spacial and temporal focusing of the potential energy during the collapse of cavitating structures. First, the initial potential energy is converted into kinetic energy in the surrounding liquid, and focused in space before the conversion to shock wave energy, and eventually to local surface impact power takes place, at the final stage to the collapse. The model is first applied to a single cavitation bubble, collapsing under ambient pressure of 1 bar, in order to investigate the source of different numerical errors made on the estimation of the radiated power during the collapse. Then, the flow around the KCD-193 model propeller is investigated and the erosion risk on the blades is assessed. The propeller performance under cavitating conditions is compared with the experimental observations and measurements, whereas the identified high erosion risk areas are compared with paint test results.
Previous studies [15,17] have shown, that the overall energy balance can be satisfied if and only if the numerical errors, involved in the reconstruction of the radiated energy, are minimized. In simple cases, where condensation can be separated from evaporation (for instance a single bubble or a bubbly cloud collapse initially at rest), those numerical errors can easily be eliminated. However, in complex flows where it is impossible to isolate a cavity collapse (flow over a hydrofoil, propeller blade etc.) the errors made on the reconstruction of the radiated energy cannot easily be avoided, resulting in eventual violation of the energy balance. Thus, it is imminent to investigate different possible approaches to reconstruct the radiated energy in such a way that all the initial potential energy is conserved.

Governing equations
The Unsteady Reynolds Averaged Navier-Stokes (URANS) equations for momentum and mass continuity to be solved, are given by where is the velocity tensor, is the fluid density, the pressure, the external force per unit mass and the viscous part of the stress tensor. The finite volume method is employed to discretize the continuous governing equations, and a segregated approach is adopted, where the flow equations are solved in a SIMPLE like manner to achieve pressurevelocity coupling. A second-order implicit method, and a second-order upwind scheme are used for the time marching and the convective terms respectively. The SST k-turbulence model is employed with a low-+ wall treatment similarly to previous studies [15,18,19]. The empirical reduction of turbulence dissipative terms in the two-phase regions has been applied, by modifying the turbulent eddy viscosity [20] = ( ) where is the vapor density, the liquid density and the mixture density. For the constant the recommended value = 10 has been used.

Cavitation modeling
In Eqs. (1) and (2), the fluid density and the turbulent eddy viscosity are given by the linear mixture relations respectively, where 0 < < 1 is the vapor fraction. A homogeneous multiphase mixture model is employed to achieve phase transition, and to track the interfaces between the two phases. The pure liquid ( = 0) and vapor ( = 1) phases are modeled as incompressible. The mixture regime is also incompressible, however, compressible behavior can be mimicked during phase transition.
An additional conservation equation that describes the transport of vapor volume fraction is solved In Eq. (6), represents the mass transfer source term. In order to account for bubble growth and collapse, a cavitation model should be introduced. The Schnerr-Sauer [21] cavitation model is used in this study, based on a simplified Rayleigh-Plesset equation, which neglects the influence of bubble growth acceleration, as well as viscous and surface tension effects. Within the control volume, the vapor phase is assumed to be present in the form of bubbles. Each bubble has the same radius . The number density of seeds is defined as 0 , which corresponds to the number of bubbles per unit volume [22]. Therefore the mass transfer source term becomes where the bubble radius can be expressed as: The seed density 0 and the seed diameter 0 = 2 0 are user specified, and the latter provides a lower limit for the bubble size. Finally the bubble radius change rate is estimated using the inertia controlled growth model where is the saturation pressure, is the local pressure around the bubble and is the fluid density.
, and , are scaling factors to adjust the source term magnitude for the condensation and the evaporation process, respectively. They work similarly to the condensation and evaporation coefficients found in other cavitation models [23][24][25].
In mass transfer models, compressibility is mimicked only in regions subjected to phase transition. Frikha et al. [26] and Morgut and Nobile [27] have shown that such models are at least able to correctly reflect the inertia driven kinematics of cavitating flows. Additionally, Koukouvinis and Gavaises [28] and Schenke and van Terwisga [29] have pointed out that the equilibrium assumption for a barotropic flow can theoretically be mimicked by the mass transfer model, if the finite transfer rate tended to infinity. This transfer rate is controlled by the mass transfer model coefficients, which in this case are the scaling factors , and , . These factors should be set as such, that the transfer rate is forced towards infinity, and thus a sharp transition is achieved from liquid to vapor phase and vice versa. Given that, the mass transfer source term always provide enough capacity to establish the equilibrium flow condition, where the time scale of phase transition is not important within the advective time scale of the flow. The density-pressure trajectory then remains close to vapor pressure during phase transition and the behavior of more realistic thermodynamic models is mimicked. Thus, in strong inertia driven flows, mass transfer models, and fully compressible models, can give very similar results as far as the inertial dynamics of cavitating flows are concerned. However, a fully compressible flow model still behaves differently in many aspects, such as wave propagation, and compressibility of pure phases and mixture.

Cavitation erosion modeling
The increased demand for the prediction of cavitation erosion has paved the way for the development of computational tools that can give a numerical estimation of high erosion risk areas. Model testing of the propeller cavitation behavior in a depressurized towing tank is nowadays the most typical way a propeller designer can get an assessment of the erosion risk on the propeller blades. However, assessing the cavitation intensity through optical observation requires a high degree of experience, and thus remains rather subjective. Besides, it does not provide more detailed information than high fidelity numerical simulations. Recent work shows that cavitation erosion risk assessment with numerical methods has a great potential, and it is expected to become integrated into the design process in the near future [15,17,18,30]. The presented method to estimate the cavitation erosion risk, has been developed within the European project ''CaFE'', as proposed by Schenke et al. [30].

Energy balance
In previous studies [15,17,19], we assumed that the potential energy within the vapor structures is instantaneously radiated in the domain, during the condensation process. However, it has been shown that the initial potential energy is partitioned into different forms of energy during the collapse, while the total energy is conserved [30][31][32]. When the vapor structure enters a pressure recovery region, then the potential energy is converted gradually into kinetic energy till the final moment of the collapse, where the initial potential energy is fully converted into kinetic energy [30]. At the final stage of the collapse, Tinguely et al. [32] showed that the initial potential bubble energy ,0 is eventually partitioned into shock wave energy , dissipative thermal energy and rebound energy , such that The dissipative thermal energy has been shown to be negligible, as thermal processes typically absorb negligible energy fractions [33]. The rebound energy depends strongly on the initial gas content and the pressure driving the collapse. As we assume that the cavitating structures are completely filled with vapor and there is no non-condensable gas in the flow, the rebound energy depends only on the driving pressure, and thus it becomes relevant only for low ambient pressures, significantly lower than 1 bar [32]. Therefore we can assume that the initial potential energy in the structure is first feeding into kinetic energy until it has fully collapsed and the kinetic energy is completely converted into shock wave energy .

Potential energy in collapsing cavitating structures
Following the notion by Hammit [9], Vogel and Laterborn [10] proposed that the potential energy of a cavity is equal to the work done by the driving pressure difference − on its vapor volume where is the ambient pressure driving the collapse and the vapor pressure. The instantaneous change of volume specific potential energy then, can be estimated in every cell at location from the material derivative of ( , ), written aṡ ℎ ∶= + ⋅ ∇. = Considering the change of potential energy in time in the domain, the potential energy is going to convert into kinetic energy. However, each term of Eq. (12) will contribute in a different way. It is crucial, at this point, to distinguish two different sources of kinetic energy. The first term on the r.h.s. of Eq. (12) includes a change in cavity volume. This will result in a collapse induced energy flux, and part of the initial potential energy in the bubble will convert into collapse induced kinetic energy, distributed to the liquid around the cavity interface. That is the kinetic energy which is responsible for the relative motion of the cavity interface with respect to the cavity (collapse) center, with velocity , proportional to the pressure difference between the ambient driving pressure and the vapor pressure . The second term on the r.h.s. of Eq. (12) is activated only when the driving pressure is time dependent in the Lagrangian frame of reference, and it represents a change in pressure. Following the cavity in time, it can experience a change in driving pressure in two ways. One is when the cavity is subjected to a steady state driving pressure gradient. In this case, the change of driving pressure can only be seen by following the cavity in space, hence in the Lagrangian reference frame. Due to the pressure gradient, part of the initial potential will convert into inertial kinetic energy. This change of potential energy can only contribute to the acceleration or deceleration of the cavity as a rigid body (or deformable body if we assume that the shape of the cavity can change, but without changing its volume). The cavity interface is then forced to move along with the center of the cavity with velocity . The second way that the cavity can experience a change in driving pressure is if the driving pressure somehow changes in the entire control volume in time. This effect can be seen in the Eulerian reference frame as well. As the pressure is uniform everywhere, and the driving pressure gradient is zero, there will be no change in the inertial kinetic energy. However, there will be a change in the potential energy of the cavity. This effect appears in the first term on the r.h.s. of Eq. (12) by expanding the driving pressure in time using a first order Taylor series approximation in the Lagrangian reference frame, as shown by Schenke and van Terwisga [17]. Fig. 1 illustrates the contribution of each term on the r.h.s. of Eq. (12) by depicting three conditions with changing potential energy. In the first example, Fig. 1a, a cavitation bubble is collapsing under steady driving pressure. Therefore the second term on the r.h.s. of Eq. (12) is zero. In this case, the collapse center does not move during the collapse, and the liquid around the cavity interface cannot accelerate without volume change. Thus, the bubble will start shrinking with collapse induced velocity , and all the potential energy contained in the bubble will convert into collapse induced kinetic energy , . In the second example, Fig. 1b, the cavitation bubble is subjected to a negative pressure gradient in the -direction. We assume that the volume of the bubble remains constant during this process. Then, the first term on the r.h.s. of Eq. (12) is zero. Consequently, all the potential energy initially contained in the bubble will convert into ''inertial'' kinetic energy , , which will accelerate the bubble towards the positive x direction. In a Lagrangian frame of reference (following the bubble motion), the bubble will experience a time dependent driving pressure. In the Eulerian reference frame this is translated as a partial time derivative term , and an advective term ⋅ ∇ (see second term on the r.h.s. of Eq. (12)). In this specific case, the driving pressure field is constant in time but not in space. This means that the partial derivative of pressure is zero, and the pressure gradient ∇ is the only responsible for the bubble motion with advective velocity . If the bubble volume is free to change its volume, the first term on the r.h.s. of Eq. (12) is not zero anymore. Then, the bubble will start shrinking, and eventually collapse, as the liquid pressure is always higher than the vapor pressure . However, the bubble will collapse at a location different than the initial location of the bubble 0 . In the last example, Fig. 1c, the flow around a hydrofoil is investigated. We assume that a cavitation bubble is shed from the sheet cavity on the suction side. The bubble reaches its maximum volume (and hence potential energy) the time instant = 0 . As the bubble is advected due to the main flow towards the trailing edge, it experiences an adverse pressure gradient. Since the pressure is higher than the vapor pressure , the bubble will start to collapse. At the same time, due to the positive pressure gradient, the fluid force is opposed to the bubble motion, resulting in a deceleration of the bubble. At a time instant = 0 + , part of the initial potential energy has been converted into collapse induced kinetic energy , , responsible for the bubble volume change, while another part has been converted into ''inertial'' kinetic energy , , responsible for the deceleration of the bubble. At the final stage of the collapse, the ''inertial'' kinetic energy will be minimum, while the collapse induced kinetic energy will be maximum, which will eventually feed into acoustic energy. Thus, only the first term on the r.h.s. of Eq. (12) is contributing to the radiated energy from the collapse of the bubble, and consequently to the surface accumulated energy.
Then, the volume specific potential energy reduction rate in every cell is given bẏ where ( ) − denotes specific volume change due to condensation.
In Eq. (13), the unknowns are the material derivative of , and the driving pressure . The material derivative of is proportional to the velocity divergence ∇ ⋅ and the mass transfer source term according to Eq. (6). Furthermore, the void fraction is defined as = − − , and from the local mass conservation + ∇ ⋅ ( ) = 0 we deduce [34]: Thus, we end up with the following three formulations for the material derivative of : All the considerations above, should theoretically give the same volume change rate. Nevertheless, each formulation introduces a different numerical error. For instance, the advective term ⋅ ∇ in Eq. (15)(a), and the velocity divergence term ∇ ⋅ in Eq. (15)(b), are discretized as one term in the transport equation of the vapor fraction, and thus, they cannot be computed separately from the solver. Nonetheless, a combined term is computed, that includes the contribution of both terms. Consequently, each one of them needs to be reconstructed. The velocity divergence term ∇ ⋅ can easily be reconstructed directly from the face fluxes. However, the reconstruction of this term introduces a non-negligible numerical error. One possible reason is that the face fluxes of are the result of an interpolation from the cell centers due to the collocated grid arrangement.
On the other hand, the advective term ⋅ ∇ involves the reconstruction of ∇ . As we use a homogeneous multi-phase method to model the different phases in the flow, the cells at the interface consist of a homogeneous mixture of liquid and vapor. In cases where a sharp interface is pursued (e.g. a single bubble or a cloud of separated bubbles) the derivatives of with respect to space are not defined, due to the sudden jump in vapor fraction from 0 to 1, and therefore, this term cannot be reconstructed. In case of a cloudy mixture (with unresolved sub-grid bubbles) the interface is more diffused and the advective term ⋅ ∇ can be computed. However, the accuracy of reconstructed gradients on unstructured meshes depends strongly on grid quality and grid effects (e.g. mesh stretching, curvature, skewness), and therefore a (sometimes significant) numerical error is introduced.
Furthermore, the mass transfer source term, on the r.h.s. of Eq. (15)(c), is obtained from the instantaneous pressure and the vapor fraction . Thereby, it is a direct outcome of the solution at each time step. In that case, the introduced error is equal to the iterative error during each time step. The source of the aforementioned errors and their magnitude are investigated in Section 4.
Finally, the volume specific potential energy reduction rate can be estimated, accordingly, by the following three considerationṡ where only the volume change due to condensation is considered in each case, and it is indicated by the subscript . In Eq. (16)(a), the advective term has been omitted, since either it cannot be reconstructed, when a sharp interface is achieved, or its error cannot be estimated, and thus we introduce an additional modeling error. For simplicity, from now on, we will refer to Eqs. (16)(a), (c) and (c), as ''partial derivative'', ''divergence'' and ''source term'', respectively.
The unsteady term of Eq. (15)(a), , represents the rate of change of in time, and it can be accurately predicted assuming a sufficient temporal discretization. In Eq. (14), moving the advective term to the right hand side gives As long as only condensation takes place (for instance a single cavitation bubble collapsing under ambient pressure higher than its vapor pressure, while far from the bubble the liquid is assumed at rest), the advective term ⋅ ∇ is zero, and thus, the partial derivative of will be equal to the material derivative of . Similarly, in cases where we know a priori that the advective term is negligible compared to the velocity divergence term (see r.h.s of Eq. (17)), then the rate of vapor reduction can be approximated by the partial time derivative of .

Effective driving pressure
The second unknown in Eq. (13) is the pressure driving the collapse . The determination of this quantity is not straightforward for complex flows, and it introduces the largest uncertainty in the prediction of the cavitation impact loads on a surface. For an isolated single cavitation bubble, collapsing in an infinite liquid volume, without the effect of gravity, the pressure effectively driving the collapse is just the pressure at infinity, which will result in a spherical collapse. Nevertheless, when the bubble is collapsing close to a wall (or close to another bubble), the driving pressure across its interface will vary due to the interaction with the wall (or the other bubble), as there is no space for the pressure to recover to ambient pressure. Similarly, in a flow around a lifting body, a hydrofoil or a propeller for instance, the interaction between the cavities formed on the blade or the hydrofoil surface, and their interaction with the wall, should be taken into account, as the pressure recovery gradients along the surface are considered important for the cavitation dynamics.
Unfortunately, the determination of the driving pressure in such complex flows is not an easy task. The local instantaneous pressure could not be used to estimate the driving pressure, because it would imply that the pressure at the cavity interface is higher than the vapor pressure, in order to obtain a non-zero driving pressure difference − , at locations where energy is radiated. However, looking at the density-pressure trajectory, it remains very close to vapor pressure during phase change, thus the corresponding pressure difference would physically be nearly zero, in case the local cell pressure is used as a measure for the driving pressure. Therefore, a different way to estimate the driving pressure should be found.
In previous studies [15,17,18], the time-averaged pressure field at time instant , , computed from the instantaneous pressure field in cavitating conditions , was assumed to be the ambient pressure field driving the cavity collapses. In this way, a rough estimate of the conditions that collapsing cavities experience on statistical average can be achieved. However, this implies that the time-averaged pressure field is sufficiently converged (almost constant in time), which always requires several cycles. In the present study, we investigate the influence of several more instantaneous pressure field definitions, computed by averaging the instantaneous pressure field over a moving time window of size Computing the moving average requires to store all the values within the chosen time window of the pressure signal in a buffer at each time step, for every cell of the computational domain. Depending on the size of the window, the amount of data can easily exceed the available random-access memory (RAM) limit [35]. Therefore, we apply the method by Welford [36] to approximate the moving average of the pressure field at time instant , , over a sliding window at each computational cell where is the number of time steps within the sliding window, and the time step size.

Radiated energy and surface impact power
Let us consider the potential energy released instantaneously at each time step and each location where condensation takes place. Then, the instantaneous radiated power is given directly by the reduction rate of the volume specific potential energy: where the minus (−) sign is a consequence of the energy conservation. However, in order to follow the energy consideration, as discussed in Section 3), the potential energy should somehow be stored as kinetic energy around the cavities, focused towards the collapse center, and eventually converted into radiated (acoustic) energy, at the final stage of the collapse, giving a more physical representation of the collapsing process of cavitating structures.
To account for the conversion of the potential energy into kinetic energy, and the focusing of the kinetic energy into the collapse center, a novel approach is employed, as introduced by Schenke et al. [30]. In this ''focusing'' model the reduction of the potential energy due to condensation is absorbed by an accumulated kinetic energy field ( , ), until a criterion for the conversion of the kinetic energy into acoustic energy is met. This process is described by an additional transport equation: and by applying the product rule, Eq. (21) becomes where is the collapse induced velocity, and is the consequent kinetic energy induced by the volume reduction of the cavitating structures in the flow. Practically, this collapse induced kinetic energy is induced from locations of negative velocity divergence only, or in general, only from locations where condensation takes place. The sum of the two terms on the left hand side of Eq. (22), i.e, the unsteady term and the advection, is the material derivative of that gives the rate of change of following a fluid particle. In particular, the term ⋅ ∇ represents the conservative advective transport of the collapse induced kinetic energy, and it is responsible for its spatial distribution around the cavitating structure as the kinetic energy is focusing towards the collapse center. However, the exact spatial distribution of ( , ) is unknown, and consequently, a modeling assumption is required for the conservative transport of the accumulated kinetic energy ⋅ ∇ .
Furthermore, the two source terms on the right hand side, represent the kinetic energy source (or creation of kinetic energy), and the radiation source (or reduction of kinetic energy), respectively. The kinetic energy source is directly connected to the reduction of potential energy. When the volume specific potential energy is reducing at any location , as a result of condensation, then the kinetic energy source term is activated. On the other hand, the kinetic energy source term remains zero, when an increase of potential energy, associated with cavity growth, takes place. Consequently, the kinetic energy source term (∇ ⋅ ) can be modeled by the volume specific potential energy reduction rate, discussed in Section 3.2.
The radiation source term,̇( ), is activated only when a certain criterion is met. To physically represent a cavity collapse, most of the energy, initially contained as potential energy in the structure, should be released at the final stage of the collapse. To chose a suitable criterion, we should take a close look to what happens as we approach the final stage of the collapse. The pressure in the mixture regime cannot be much higher than the vapor pressure, and therefore, a high amplitude pressure wave can only form in the liquid phase. Besides, the absorbed kinetic energy should accumulate on the low pressure side, and propagate inwards, towards the collapse center. Considering these assumptions, we derive the following criterion for : According to Eq. (23), the parameter is defined such that = 1 at locations where the vapor has completely disappeared and the local pressure is higher than the pressure at infinity. Considering also the fact that the value of is of interest only at locations where the potential energy has been converted into collapse induced kinetic energy, then we can isolate all the locations where a collapse of a vapor structure has taken place. Consequently, a volume cell, where there is some collapse induced kinetic energy focused, is a candidate to radiate energy, if and only if its vapor volume has fully condensed and if a pressure wave is emitted. Then Eq. (22) can be replaced by the following modeling assumption: The term ( ) is a model for the unknown conservative advective transport of the kinetic energy ⋅ ∇ in Eq. (22). The terṁ, ( ) represents the source term for the reduction of potential energy, and consequently the generation of kinetic energy (∇ ⋅ ) in Eq. (22). Hence, ( ) −̇, ( ) serves for the kinetic energy flux. The term is a model for the radiation source terṁ( ) in Eq. (22). The energy radiation is modeled as a discrete event, because all the focused energy is released within the time interval , which means that the energy content of the radiated wave is conserved, while the exact wave shape across the wave front is not captured. Finally, the term (1 − ) in Eq. (24) has been added to make sure that at locations where volume specific energy is radiated, there is no production of kinetic energy, and therefore the overall energy balance is satisfied.
The advective transport term ( ) in Eq. (24) is given by [30] ( ) =̇, − P (∇ ) where is assumed constant in space and is defined such that and P (∇ ) is the normalized projection of ∇ on the local flow velocity vector where division by zero is prevented by adding a small number to the denominator ( = 10 −15 ), such that it does not affect the accuracy of the final result. In Eq. (25), the first term on the right hand side is associated with the production of collapse induced kinetic energy, and it is proportional to the reduction rate of potential energy. The second term is related to the transport of the kinetic energy, and it is based on the assumption that the flow around the interface of the collapsing cavity is directed towards the center, thus aligned with ∇ , considering that is stored at the cavity interface. To get the collapse induced kinetic energy at each time step , we explicitly forward the solution in time using a first order Taylor expansion, which gives and the volume specific radiated power is given bẏ Knowing the radiated volume specific power at each cell in the domain, and assuming that each point source emits its converted potential energy as a radial wave with an infinitely large propagation speed, the instantaneous surface specific impact poweṙ( , ), at some surface location is given by Schenke and van Terwisga [17] aṡ where is the surface normal vector at the impact location. Finally, cases where the cavity collapse rate goes from negative to zero before the cavity is fully collapsed, and the condensation process stops, or it is even reversed and evaporation takes place, are not considered in this model. In such cases, the kinetic energy accumulated at the cavity interface should be dissipated or even converted back to potential energy, which cannot be handle by the model as described in this section. However, such cases are assumed to occur quite rarely when the cavities undergo high pressure differences.

Single cavitation bubble collapse
A verification study is conducted on a single cavitation bubble collapse, considering an inviscid flow. This study serves for verification of the cavitation and erosion model. First, an idealized spherical bubble is simulated in an infinite liquid. The sensitivity of the collapse time to the condensation scaling factor , , the time step size , and the grid size has been investigated. The obtained converged solution is then compared to the analytical solution of the Rayleigh-Plesset equation [37]. Then, the focusing model is applied to verify that the potential energy is stored at the interface as kinetic energy, until the final stage of the collapse, when it is radiated to the domain as acoustic energy. Subsequently, we apply an inflow velocity at the one boundary, parallel to the -axis, so that the flow around the bubble is advected downstream during the collapse. This configuration serves for identification of the numerical error sources on the prediction of the total radiated energy . As mentioned in Section 3 we consider three alternative methods to estimate the change in potential energy during condensation. The volume reduction rate can be estimated by using the partial derivative of the vapor fraction, the velocity divergence, or the cavitation mass transfer source term. In the case where the flow is initially at rest, the total volume change can be accurately predicted by the volume integration of the partial derivative of the vapor fraction ∫ ( ) , as the advective term ⋅ ∇ is zero then. However, it is interesting to see how this term behaves when there is large advection in the flow, and how it compares with the other two considerations for the prediction of the volume change.
Finally, we let the bubble collapse close to an infinite flat surface to verify that half of the energy is accumulated on the surface and the total energy is conserved [17].

Sensitivity study
The collapse of an isolated single cavitation bubble with a radius of = 3.84 mm has been simulated. Six different meshes are used for the sensitivity studies. Table 1 shows the number of cells per diameter, the total amount of cells and the volume off-set from the theoretical bubble volume. For all grids, apart from the coarsest grid (Grid 1), the volume off-set is less than 1%. Fig. 2 shows the domain size and the grid refinement around the bubble in Grid 5. The domain boundaries are located 250 radii away from the bubble center.  A pressure outlet boundary condition has been applied to all boundaries with ∞ = 1 bar, while the vapor pressure inside the bubble is equal to = 2340 Pa (see Fig. 3). The pressure field in the rest of the domain is initialized such that it satisfies the Laplace equation ∇ ⋅ ∇ outside the bubble surface. An initialization algorithm has been used to compute the Laplacian pressure field and the vapor volume fraction at each cell for each grid. The initial vapor fraction field is such, that all the cells entirely located inside the bubble are assigned a vapor fraction of = 1, while those entirely located outside the bubble are assigned a vapor fraction of = 0. To determine the vapor fraction of the cells at the interface, a sample algorithm has been used. The algorithm divides every cell that is cut by the bubble interface in 100 3 subcells and determines the amount of the subcells that are entirely inside the interface, assigning a percentage of vapor to each cut cell. The initial vapor fraction field for the second coarsest (Grid 2) and second finest (Grid 5) mesh is depicted in Fig. 4, where the better representation of the bubble interface (red line) can be seen as we refine the mesh around the bubble. The bubble interface is shown as an iso-line of = 0.5. The sensitivity of the collapse time to the condensation scaling factor , (see Eq. (9)), and the temporal and spatial resolution has been investigated. Fig. 5 shows the evolution of the dimensionless vapor volume for different values of the condensation scaling factor , , and a fixed time step size = 1 × 10 −6 , in grid 3. For very small values, a significant delay of the collapse time is observed, while for , ≥ 1 the solution becomes independent of the coefficient value. However, this behavior can only be achieved if the temporal and spatial resolution are sufficient. Fig. 6 depicts the evolution of the dimensionless vapor volume in time, for a constant condensation scaling factor , = 10, and a systematic variation of the time step size in grid 3. Surprisingly, the collapse time is not very sensitive to the temporal resolution, as already with a time step size of = 1×10 −5 , a time step size independent solution is achieved. In other words, we need no more than 40 time steps during the collapse, in order to obtain a time step size independent solution. The sensitivity of the grid size is depicted in Fig. 7 for 6 different grid sizes. It shows the evolution of the dimensionless vapor volume for , = 10 and = 1 × 10 −6 , and how they compare to the Rayleigh-Plesset equation. As we move to a grid density higher than grid 2, the collapse time does not show high sensitivity. Looking more thoroughly at the predicted collapse time, grid 5 gives a fully converged and grid size independent solution, which deviates less than 1% from the collapse time obtained from the Rayleigh-Plesset equation.

Energy consideration
As discussed in Section 3, in this study we assume that the initial potential energy ,0 , contained in the vapor structures, is continuously feeding into kinetic energy , during the condensation process, until eventually it is fully converted into radiated energy , at the final stage of the collapse. In this Section, we compare the energy focusing approach, with the original consideration that energy is continuously released in the domain, and the potential energy is instantaneously converted into radiated energy (see Eq. (20)). Fig. 8 shows the evolution of the radiated energy, normalized by the initial potential energy, in time for the non-focusing and the focusing approach. In the non-focusing approach, the instantaneous potential power is computed using the three different considerations, as explained in Section 3.2. The volume integrated rate of vapor reduction is predicted (see Eq. (16)) using the partial derivative (red), the velocity divergence (green), and the cavitation mass transfer source term (blue dotted). While the total energy in the domain is conserved, when the partial derivative and the source term are considered for the prediction of the instantaneous power, with the velocity divergence, part of the energy is being lost and only 63% of the total energy is finally being released. The source of this error is the divergence term ∇⋅ . The velocity divergence cannot be computed directly from the solver, as has been demonstrated in a previous study [15], and as a result it needs to be reconstructed. The reconstruction of the velocity divergence introduces a numerical error, which results in the error shown here. In this particular test case, this error could be eliminated by correcting the velocity divergence field by a constant factor c, determined at each time step as reported by Schenke and van Terwisga [17], as the vapor volume change is subjected only to the condensation process during the cavity collapse, and there is no advection in the flow. However, in most of the cases the condensation process cannot be estimated from the total volume change and this correction is not possible.
The focusing approach has only been applied in the case where the mass transfer source term is used to predict the vapor volume reduction rate. It is demonstrated that the potential energy, initially contained in the bubble, is successfully converted into kinetic energy during the condensation process, until we approach the final stage of the collapse where all the energy is radiated radially to the domain and the total energy is conserved. We could also use the partial derivative, however, this would hold for that particular case only, due to the fact that the flow is at rest and no advection takes place. The partial derivative term includes an advective contribution, which in this case is negligible.   Nevertheless, in case we would have a non-zero advective velocity, then the predicted total energy would not be conserved. This can be shown by applying an inflow velocity in positive x direction, as shown in Fig. 9. Now, the computed partial derivative of the vapor volume includes an advective contribution which cannot be neglected. The error made in the energy conservation is getting more dominant with increasing the inflow velocity, as illustrated in Fig. 10. The total energy released in the domain after the bubble has collapsed is overpredicted, while with the material derivative consideration is still underpredicted. However, in the latter case, the magnitude of the error is diminishing as the flow velocity increases. Finally, when the volume integrated rate of vapor Fig. 7. Evolution of the dimensionless vapor volume, normalized by the initial vapor volume, over time for = 10, = 1 × 10 −6 , and a systematic variation of the mesh density. Fig. 8. Total radiated energy, normalized by the initial potential energy, with the non-focusing and the energy focusing model. For the non-focusing model, the volume integrated rate of vapor reduction is predicted using the partial derivative (red), the velocity divergence (green), and the mass transfer source term (blue). For the energy focusing model, the volume integrated rate of vapor reduction is predicted using the mass transfer source term only (black). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 9. Boundary conditions when a non-zero inflow velocity is applied to the domain. A velocity inlet boundary condition is assigned to the left boundary, and a flow advection is forced in positive x direction. reduction is predicted from the mass transfer source term, all the energy is fully conserved, regardless of the magnitude of the flow velocity.

Collapse near an infinite flat surface
The collapse of the same spherical cavitation bubble close to an infinite flat wall has been simulated (see Fig. 11). The bottom boundary has been placed close to the bubble at a distance = 5 mm ( = = 1.3) from the bubble center, and a slip wall boundary condition has been assigned. In the previous section, 4.1.2, we showed that all the initial potential energy is released and converted into shock wave energy during the collapse, such that the total energy is conserved. It can be argued that half of this energy should be distributed to the lower surface of infinite dimensions. Fig. 12 shows the total accumulated surface energy over time, normalized by the initial potential energy, obtained from the non-focusing cavitation intensity model, and the energy focusing approach. The surface energy is estimated at twelve time instants during the collapse. Similar to Fig. 8, the total energy is conserved when the volume integrated rate of vapor reduction is predicted from the partial derivative or the mass transfer source term, as half of the initial potential energy is transmitted to the surface. When the volume change is predicted from the velocity divergence term, part of the energy is lost, due to numerical error in the reconstruction of the velocity divergence, however the final energy, distributed on the flat surface, is still half of that in the case of the bubble collapse in an infinite liquid.
The accumulated specific energy on the surface can be estimated for the two different approaches, which should give an indication of the cavitation erosion intensity of the impact. Fig. 13 depicts the distribution of the accumulated surface specific energy, after the bubble has collapsed. In both cases, there is a distinct axisymmetric footprint with a clear peak value in the center. However, the magnitude of the maximum specific energy in the middle, is higher in the focusing approach. This can be explained by the fact that in this approach the potential energy is continuously focused into the collapse center, therefore right before the collapse, all the energy is concentrated at a single point, and is released all at once, giving a footprint with high intensity towards the middle. While with the non-focusing approach, the energy is continuously being released in the domain, such that the energy remained in the cavity right before the collapse is less than in the focusing approach. This results in an accumulated surface energy, which smears out throughout the collapse, and the energy is continuously distributed over a larger area. Fig. 10. Total radiated energy, normalized by the initial potential energy, with the non-focusing and the energy focusing approach for a systematic variation of the inlet flow velocity. For the non-focusing model, the volume integrated rate of vapor reduction is predicted using the partial derivative (red), the velocity divergence (green), and the mass transfer source term (blue), while for the energy focusing model the volume integrated rate of vapor reduction is predicted only using the mass transfer source term. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 12. Total accumulated surface energy over time, normalized by the initial potential energy, with and without focusing. For the non-focusing model, the volume integrated rate of vapor reduction is predicted using the partial derivative (red-circle), the velocity divergence (green-square), and the mass transfer source term (bluetriangle). The latter is also used to predict the accumulated surface energy with the energy focusing model (black-diamond). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Case description
The cavitation erosion intensity is estimated on the King's College-D (KCD)-193 model propeller. This propeller has been tested at the Emerson Cavitation Tunnel of Newcastle University. Experimental paint test has been conducted on the propeller blades [38]. A 2D wake screen was used, located 0.4572 m (1.5 ) from the propeller center, to create a non-uniform inflow to the propeller, and test the propeller performance under more realistic conditions. The propeller geometry, and the velocity and turbulence intensity measurements at the 2D wake screen plane, were provided by the University of Strathclyde, Glasgow [7]. Table 2 shows the propeller characteristics. The propeller was tested in atmospheric condition (see Table 3).  the experimental set-up by Mantzaris et al. [38]. The location of the inlet boundary is located 1.5 diameters from the propeller center, where the non-uniform inflow velocity is specified only in the -direction. The size of the rotating domain (light blue) is chosen such, that the interfaces are not located in areas, where high pressure gradients occur. The propeller rigid motion is simulated using a sliding mesh approach. In Fig. 15 the velocity distribution in the inlet is shown, and in a plane just in front of the propeller hubcap. The flow is slightly accelerated by the propeller motion, before entering the rotating domain, however, the character of the non-uniform overall flow is maintained. A pressure boundary condition, with a constant ambient pressure, is applied 15 diameters downstream, while the top, bottom and the side boundaries are set to slip wall condition, as they are assumed to not interfere with the propeller motion. The Reynolds number at 70% of the radius is of the order of 10 6 (see Table 3), which is higher than the critical Reynolds numbers for transition from laminar to turbulent flow (typically 10 5 < < 10 6 ). Besides, the non-uniform inflow reduces even more the risk of regions with laminar flow, and the flow over the whole blade can be considered, more or less, fully turbulent.
For the grid generation, trimmed hexahedral cells are used with local refinements and prism layers along the wall to resolve the viscous sub-layer. First a coarse grid, with a refinement at the propeller leading and trailing edge, is generated to roughly estimate the amount of vapor on the propeller blades. Then, an additional refinement is applied on one of the blades to assess the cavitation intensity on the surface, in order to follow the guidelines for the spatial resolution, as proposed in Section 4. First, the size of the cavitating structures needs to be estimated, so that the required grid resolution and time step size can be determined. In Fig. 16, a 2D section of the developed cavity is T. Melissaris et al.  shown, approximately at the location on the blade, where the maximum amount of paint was removed during the paint test. The width of the cavity in the -direction is estimated to be around 1.2 cm. The corresponding Rayleigh collapse time is found to be ≃ 4.27 × 10 −4 s. In Section 4 we showed that approximately 18 cells per diameter, and about 40 time steps, are required to get a grid and time step size independent solution for the single bubble collapse test case. For a cavity with diameter = 1.2 cm, this translates to a minimum cell size of ≃ 0.667 mm, and a time step size of = 1.28 × 10 −5 or a rotation rate of 0.115 deg per time step. Eventually, the minimum cell size of the current domain is selected, ≃ 0.45 mm, which corresponds to approximately 27 cells over the cavity width, and a time step size of = 1.11 × 10 −5 , or a rotation rate of 0.1 deg per time step. Then, the refinement is applied at locations where the maximum value of over one rotation is not zero, and only on one blade to reduce the computational cost. This results in a refinement region that covers the whole cavity in the vicinity of the blade surface, at any time instant, during a full rotation. Fig. 17 shows the mesh refinement and the cavity size when the refined blade is located at the top position (0 • ).
Three finer grids have been generated to assess the temporal and spatial discretization error. In order to keep the different grids as geometrically similar as possible, special attention has been paid to the   generation of the prism layers close to the blade surface (for details see Crepier [39], Melissaris et al. [15]). Eca et al. [40] have shown that for the SST k-turbulence model at high Reynolds number, wall + equal to 1 is not sufficient to obtain low numerical uncertainties, and + values as low as 0.1 (or at least below 0.5) are required to obtain numerical uncertainties similar to other turbulence models. Thus, average wall + values at the refined blade are kept lower than 0.25 for all grids. Table 4 shows, for each grid, the number of computational cells in the whole domain, as well as on the refined blade, and the number of prism layers, and the average wall + values on the refined blade. Finally, Fig. 18 shows the wall + distribution on the refined blade, for the coarsest (Grid 1) and the second finest grid (Grid 3). Clearly, the + values drop as we go to finer grids and more prism layers, even below 0.1 for the finest grid (Grid 4, see Table 4).
Finally, the condensation and evaporation coefficients, , and , , are selected based on the results from the single bubble collapse test case (Section 4). Coefficient values equal or larger than one gave solutions, which were independent of these coefficient values. However, in the propeller test case, we observed that values higher than one introduce some numerical instabilities to the solver, and thus the default value of one was chosen for both coefficients for all cases.

Propeller performance and numerical uncertainty
The propeller performance is predicted and compared with experimental measurements as reported in Usta et al. [7]. Propeller thrust and Table 5 Propeller performance for the coarsest grid (Grid 1) and four different time-step sizes, and comparison with experimental measurements.  torque are computed for the refined blade, and then multiplied by the number of blades. The numerical uncertainty is assessed following the procedure by Eca and Hoekstra [41], using the Validation & Verification tool, developed by MARIN. 1 In Fig. 19 the cavitation development on the propeller blades is compared with the experimental observation. Because of the non-uniform inflow, the vapor volume varies at different blade positions, and unfortunately, due to the limited documentation of the vapor volume development during the cavitation test, it is hard to compare between the experiment and the simulated cavity. However, qualitatively, the cavitation extent seems comparable. First the temporal discretization error is estimated. Table 5 shows the propeller thrust and torque for the coarsest grid (Grid 1), and different time step sizes, and how they compare with the experiment. The convergence of the propeller thrust and torque with respect to the rotation rate, is depicted in Figs. 20 and 21 respectively. The corresponding numerical uncertainty is indicated by an interval that contains the exact solution with 95% coverage. Propeller thrust and torque are not very sensitive to the time step size. The deviation between the computed thrust and torque, and the experimental values is less than 2% even for the highest time step size, while the numerical uncertainty is estimated less than 0.5% for all time step sizes.
The grid discretization error is also assessed for each grid. Table 6 compares the propeller performance for each grid, and Figs. 22 and 23 1 http://www.refresco.org/verification-validation/utilitiesvv-tools/.  Fig. 22. Convergence of the propeller thrust with the refinement ratio ℎ ∕ℎ 1 . Impression of the numerical uncertainty estimate, and comparison with the experimental measurement. depict the convergence of propeller thrust and torque with respect to the refinement ratio ℎ ∕ℎ 1 , and the corresponding numerical estimates.
Here, ℎ is the typical cell size of the ℎ grid, and ℎ 1 the typical cell size of the finest grid. Propeller thrust is slightly more sensitive to the grid density. The observed order of convergence is equal to the theoretical one ( = 2), and the numerical uncertainty is below 3% for all grids, while it drops below 1% only for the finest grid, probably because the average wall + is below 0.1 only for Grid 4. On the other hand, propeller torque is not that sensitive to the grid size. The observed order of convergence is close to the theoretical one ( = 1.61), and the numerical uncertainty is very low for each grid (about 1% for Grid 1 and 0.3% for Grid 4).

Erosion risk assessment
The cavitation intensity on the propeller blades is estimated by determining the local surface impact poweṙ( , ), according to Eq. (30). First, an investigation on the effect of the pressure field effectively driving the collapses has been conducted. The accumulated volume specific potential energy in the domain is estimated using the mass transfer source term and different driving pressure fields. The driving pressure fields are obtained by averaging in time the local cell pressure over different time windows. Fig. 24, demonstrates the influence of the different driving pressure fields, as computed for different sliding windows, on the total potential energy. The total accumulated potential energy obtained from the instantaneous pressure field, is compared with the one obtained from pressure fields averaged over a time window of 1 deg, 10 deg, 30 deg, 90 deg, 360 deg, and eventually, the pressure field obtained from a time averaging over the total sample time. We observe, that when the chosen time window is smaller than one full propeller rotation, the predicted energy is quite sensitive to the window size. As the time window decreases, the predicted energy is approaching what we obtain when the instantaneous pressure is used as driving pressure. Thereby, high amplitude pressure peaks are expected by the periodic cavity collapses and thus, we cannot account for pressure recovery gradients, as typically present along lifting bodies. As the time window is getting larger, so does the total accumulated volume specific potential energy. Eventually, when a time window of one whole revolution is used, the predicted energy is almost identical to the one obtained with a time window of the total sample time. Thus, a time averaged pressure field over one rotation, or one shedding cycle, is considered sufficient for an estimate of the effective driving pressure field.
Afterwards, the distributions of the accumulated surface specific energy on the refined blade, obtained by the non-focusing and the energy focusing model, after eight propeller revolutions, using the three different considerations for the vapor volume reduction rate are compared. Fig. 25 depicts the distributions of the non-focusing model. It is obvious that there is a huge discrepancy between the distribution obtained with the partial derivative and the one obtained with the source term. This is explained by the fact that the partial derivative of in Eq. (15)(a) includes an advective term. As the flow accelerates in higher radii, the advective term contributes more and more to the radiated energy, and as a consequence, it deviates from the radiated energy obtained with the source term. The higher the velocities occurring around the blade, the larger the error. Similarly, the distribution obtained with the divergence, still shows a discrepancy from the one obtained with the source term, although quite smaller. The reason for that is that the numerical error involved in the reconstruction of the radiated energy using Eq. (15)(b) is getting smaller as the flow velocities get larger, and the inertia in the flow increases.
To further analyze the errors made in the reconstruction of the radiated energy in a more quantitative manner, the total accumulated surface energy per rotation is obtained by the surface integral of the accumulated surface specific energy at the end of each rotation and it is depicted in Fig. 26. The total accumulated energy obtained by the partial derivative of is almost an order of magnitude smaller than the other two, while the accumulated energy obtained by the divergence term is about 20% smaller than the one obtained by the mass transfer source term..  Now, comparing the impact distribution between the non-focusing model, where the potential energy is radiated instantaneously, Fig. 25, and the energy focusing model, Fig. 27, the energy appears to be accumulated in a much larger area, with the non-focusing model. Any dynamics of the sheet cavity, result in energy transfer to the surface, thus leaving a footprint, which extends almost over the whole blade along the cavity trailing edge. On the other hand, the surface impacts are much more scattered, and shifted towards the trailing edge, when the energy focusing approach is applied. In addition to that, these localized events are expected to be more aggressive, as big amounts of energy are concentrated into very small areas of the blade surface.
Looking at the total accumulated energy per rotation as obtained by the energy focusing model for the three different energy reconstruction methods (see Fig. 28), one can notice that the error magnitudes are different than the ones with the non-focusing model. The accumulated energy obtained by the partial derivative is higher than before, while the one obtained by the velocity divergence is smaller, which shows that the errors involved in these two energy reconstruction methods are highly dependent on the flow characteristics. As shown in the single bubble collapse test case, when the inertia in the main flow increases, the error involved in the reconstruction of the radiated energy using the partial derivative of and the velocity divergence are getting larger and smaller respectively. The total inertia in the flow depends on the flow velocity and the collapse induced velocity. While the collapse induced velocity is the highest at the final stage of a collapse, the flow velocity is almost zero as the collapsing cavities are decelerated by the adverse pressure gradient. Thus, the total inertia at locations where we approach the final stage of a collapse is lower than locations with  higher flow velocity but low collapse induced velocity. Since with the energy focusing model, the radiated energy is reconstructed only at locations of the final stage of the collapse the errors involved in the reconstruction of the radiated energy using the partial derivative will be smaller and the ones using the velocity divergence will be larger. All in all, we see a similar behavior as in the single bubble collapse in Section 4, and the radiated energy reconstruction errors obtained from the partial derivative of and the velocity divergence are strongly dependent on the inertia of the flow.
Figs. 29 and 30 show the surface integrated accumulated energy per rotation for eight propeller revolutions for the two grids, using the non-focusing and the energy focusing model respectively. The radiated energy on the blade is rather insensitive to the grid density, for both models. Similar energy content is predicted on the blade surface during each revolution for both grids. These results are well in line with the findings in Section 4, and the resolution of the coarsest grid is already sufficient to obtain a grid independent solution for the amount of accumulated energy on the blade. However, the total surface accumulated energy obtained from the energy focusing model is about an order of magnitude less than the one obtained from the non-focusing model. When no potential energy focusing is considered, the energy is radiated instantaneously, when there is a reduction of vapor volume. On the other hand, with the energy focusing model, energy is radiated only from structures that actually reach a final collapse stage, resulting in less accumulated energy on the surface.
Furthermore, some small variations are observed between the propeller revolutions. These variations are partly related to cloud cavitation and system instabilities. One can show from the localized Euler equations at the cavity closure that the stagnation point is highly unstable, even for a globally steady cavity flow [2]. The instability affects the re-entrant jet and the whole closure region, leading to irregular break-up patterns. It appears that the re-entrant jet and cloud cavitation instability preferably occurs in short cavities with pronounced adverse pressure gradients, whereas system instability mostly effects long cavities with weak adverse pressure gradients. Long cavities can be very sensitive to variations in the incoming flow, which is the case for a propeller operating in a non-uniform inflow. These 3-D aspects of cloud shedding are referred to as being intrinsic instabilities or often referred to as self-excited instabilities [2], and it has been shown that they are essentially inertia controlled, and developed naturally, separated from instabilities due to turbulence [42]. In a RANS modeling, most of those instabilities are restrained due to the increase of eddy viscosity at the mixture interface. However, employing the correction to suppress the eddy viscosity at the interface, we allow for partial cavity shedding and therefore we introduce intrinsic instabilities. The other part of these variations at each propeller revolution is related to modeling and iterative convergence errors. Therefore, more than one propeller rotations are necessary in order to obtain a reliable estimate for the average accumulated surface energy. Still, these discrepancies are rather small, and already after the first revolution we get a good idea of the amount of the accumulated energy on the blade.
Comparing the accumulated energy distribution on the blade obtained by the non-focusing model for Grid 1 and Grid 3 after 8 propeller rotations (Fig. 31), they look quite similar. On the other hand, when we compare the distributions obtained by the energy focusing model, some differences are rather noticeable. While for the coarse grid (Grid 1) most of the energy is radiated towards the blade trailing edge, for the fine grid (Grid 3) some part of the initial potential energy is distributed closer to the leading edge. This is because small scale properties are typically dependent on the chosen resolution in unsteady solution of T. Melissaris et al. cavitating flows [42], and thus with the finer grid we resolve more vapor structures in the domain than in the coarse grid. To generate the coarse grid, we estimated the width of the sheet cavity close to the trailing edge, at the area where the maximum paint removal was observed during the paint test. Therefore, the grid resolution in areas closer to the leading edge, where the width of the sheet cavity is thinner, is perhaps not sufficient to resolve any possible collapse, but only bigger structures that are shed towards the trailing edge. However, the finer grid probably provides this resolution to resolve some of the smaller structures closer to the leading edge. On top of that, due to the smaller cell size in the finer grid, the energy at the final stage of the collapse is focused on a much smaller area, and therefore the radiated energy is more localized to smaller surface areas. That explains why for the coarse grid the high energy density areas are broader than for the fine grid. However, the total energy on the blade for both grids is similar, since large scale properties such as shedding frequencies and characteristic void fraction distributions, regarding the basic shape of the sheet and cloud cavities, seem to be less dependent on the grid density.
Furthermore, we compare the distributions of the accumulated surface specific energy on the refined blade, obtained by the energy focusing model for the coarse and the fine grid (Grid 1 and Grid 3 respectively), with the one obtained from the paint test (see Fig. 32). Although the impacted areas are not as extensive as in the test, this is reasonable since the erosion pattern on the blade was obtained after half an hour of testing, and thousands of propeller revolutions, while we assessed the impact distribution after only 0.32 s and 8 revolutions. Due to the intrinsic instabilities of the partial cavity dynamics and the modeling/numerical instabilities we introduce, the location of the final collapse of shed cavities of smaller scale, will vary between different propeller revolutions and, thus much more time is needed to get a fully converged impact distribution on the blade, especially for the finest grid, where more small scale structures are resolved. However, in both grids, high energy density is obtained in blade areas, where paint was removed during the test, and the main impact areas are captured. Finally, it is important to mention that areas with high energy density are not necessarily areas with the highest erosion risk. As we have indicated in previous studies on hydrofoils [15,17], the impact power is of much higher importance. Surface areas with similar high energy content can be observed either by repetitive low amplitude events, or by extreme events of lower frequency. While in both cases the energy content might be comparable, in an extreme event large amount of energy is released, and distributed to the surface, during a small period of time, which might render such an event more aggressive. Such analysis of the rapidness of the cavity collapses could be done using aggressiveness indicators (see Schenke and van Terwisga [17], Melissaris et al. [15], Schenke et al. [35]), however that is beyond the scope of this paper.

Conclusion
A novel methodology, developed within the European project ''CaFE'', by Schenke et al. [30], has been used in this study to predict cavitation implosion loads on a surface. This model takes into consideration the time accurate energy balance during a cavity collapse. While the cavity starts to condense, the initial potential energy is continuously feeding into kinetic energy, around the cavity interface. At the final stage of collapse, all the potential energy has been converted into kinetic energy, and has been focused into the collapse center. That is the point, when all the energy is released in the domain as shock wave energy.
The model has been verified on a single cavitating bubble collapsing in an infinite fluid and close to a wall, where an investigation on the numerical errors involved in the reconstruction of the radiated energy has been conducted. The vapor volume reduction rate is predicted using three different considerations. One involves the partial derivative of the vapor fraction , the second involves the velocity divergence, and the last one involves the mass transfer model source term. The model is further applied to the KCD 193 model propeller test case, that was tested in the Newcastle cavitation tunnel behind a 2D wake screen. From these studies, we draw the following conclusions: • When the vapor volume reduction rate is estimated using the cavitation model mass transfer source term, then the numerical error, involved in the reconstruction of the radiated energy, is minimal. All the energy is then conserved and the total energy balance is fully satisfied. The error (modeling or numerical) in the other two considerations depends strongly on the inertia of the flow. They could provide a good estimate of the surface accumulated energy in specific cases, however one should be aware of their sensitivity to the flow conditions. • The temporal and spatial discretization errors have been assessed for the propeller performance. Propeller thrust and torque are rather insensitive to the time step size and low uncertainty estimates are obtained (below 0.5%). Propeller thrust is more sensitive to the grid density than torque. However, the uncertainty of the propeller thrust is below 3% for all grids, and it drops below 1% for the finest grid. The uncertainty of the propeller torque is rather low for all grids (about 1% or lower). • The moving time average of the instantaneous pressure field over a time window equal to one shedding cycle (or one revolution) or higher, is considered sufficient for an estimate of the effective driving pressure field. It has been shown that when the sliding window is equal to one shedding cycle or higher, the difference on the total accumulated potential energy in the domain is negligible. Nevertheless, as long as the driving pressure distribution is not exactly known, some uncertainties regarding the potential energy content remain. • The total accumulated surface energy on the refined blade was found to be insensitive to the grid density for the non-focusing and the energy focusing model. Comparable amount of energy was obtained for the coarse and the fine grid (Grid 1 and Grid 3) during each propeller revolution. However, the total surface accumulated energy obtained from the energy focusing model is about an order of magnitude less than the one obtained from the non-focusing model. In the latter, each volume variation results in radiated energy, leading to over-prediction of the surface accumulated energy. • When the potential energy focusing is considered some discrepancies are observed on the impact distribution between the coarse and the fine grid, which are related to the grid resolution. The finer mesh resolves smaller scale structures, and provides sufficient resolution for thinner parts of the sheet cavity, while the energy at the final stage of the collapse is focused to a much smaller area. Therefore, the events are more localized, and the energy scattering is increased, so does the time to get a fully converged impact distribution on the blade. Yet, the impact distribution obtained by the energy focusing model agrees well with the high erosion risk areas, as indicated from the paint test, for both grids, and high energy density is observed in the areas of maximum paint removal. The basic cavity dynamics (shape of sheet and cloud cavities) are well captured even with the coarser grid, while, in the finer grid some part of the energy is distributed closer to the leading edge, however it does not necessarily indicate that the erosion risk in those areas is high, as the rapidness/aggressiveness of the collapsing events has not been taken into account. • Finally, the capabilities of our fully energy conservative model to predict the implosion loads of large scale cavitating structures on propeller blades illustrates its potential for engineering applications. Nonetheless, a good estimate of the size of the collapsing cavitating structures to be resolved is a prerequisite for this model, in order to apply the required grid and temporal resolution, as proposed in this study.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.