## Abstract

This paper is intended to provide a potential basis for a numerical prediction of cavitation erosion damage. The proposed method can be divided into two steps. The first step consists in determining the loading conditions due to cavitation bubble collapses. It is shown that individual pits observed on highly polished metallic samples exposed to cavitation for a relatively small time can be considered as the signature of bubble collapse. By combining pitting tests with an inverse finite-element modelling (FEM) of the material response to a representative impact load, loading conditions can be derived for each individual bubble collapse in terms of stress amplitude (in gigapascals) and radial extent (in micrometres). This step requires characterizing as accurately as possible the properties of the material exposed to cavitation. This characterization should include the effect of strain rate, which is known to be high in cavitation erosion (typically of the order of several thousands s^{−1}). Nanoindentation techniques as well as compressive tests at high strain rate using, for example, a split Hopkinson pressure bar test system may be used. The second step consists in developing an FEM approach to simulate the material response to the repetitive impact loads determined in step 1. This includes a detailed analysis of the hardening process (isotropic versus kinematic) in order to properly account for fatigue as well as the development of a suitable model of material damage and failure to account for mass loss. Although the whole method is not yet fully operational, promising results are presented that show that such a numerical method might be, in the long term, an alternative to correlative techniques used so far for cavitation erosion prediction.

## 1. Introduction

Hydraulic devices such as hydroturbines, propellers or pumps may be affected by cavitation. Cavitation is the development of vapour bubbles in the flowing liquid. Bubbles appear in regions of low pressure generally associated with high velocity, as stated by the Bernoulli equation, or strong curvature of the streamlines. In a first approximation, cavitation bubbles appear where the local pressure drops below the liquid vapour pressure, which is the critical pressure for cavitation inception. Cavitation bubbles also contain non-condensable gas that diffuses into them from the liquid where it is present as dissolved gas. Cavitation bubbles are entrained by the liquid flow and their vapour content changes back into liquid as they enter higher-pressure regions.

This collapse of vapour structures is generally a violent phenomenon. This is because the counter-pressure inside the bubble, which is close to the vapour pressure, is quite small (2300 Pa at room temperature). As a result, the liquid around the bubble is strongly accelerated and reaches locally high velocities. A well-known consequence of cavitation is erosion of neighbouring solid walls caused by the impact of high-velocity liquid jets. Figure 1 shows a typical example of cavitation erosion damage observed in a gear pump. Several regions of erosion with a different amount of damage are visible in figure 1, from limited damage characterized by an orange peel appearance to more serious damage characterized by significant material removal.

Cavitation damage depends on both liquid flow conditions and material properties. On the liquid side, a key issue is to define the so-called flow aggressiveness or cavitation intensity. An attempt is made in this paper to quantify the cavitation intensity in terms of the loading conditions due to bubble collapse.

On the material side, damage depends upon the mechanical properties of the material. Various correlations are available in the literature for predicting cavitation damage (e.g. [1]). They generally relate the resistance to erosion, often defined as the inverse of the erosion rate (typically in micrometres per hour), to some relevant material property (e.g. hardness). A major difficulty is that such correlations are not universal and depend upon the class of material and type of cavitating flow. It is then hazardous using such correlations to predict cavitation damage for materials or cavitating conditions that deviate from those used for establishing the correlation.

Because of the lack of versatility of correlations, alternative methods are being developed that make greater use of numerical simulation. They are based on a physical analysis of the phenomena involved in cavitation erosion and require modelling as accurately as possible each step of the erosion process for both fluid and material. The objective of this paper is to provide an outline of such an analysis.

The proposed method can be divided into two main steps. The first step consists in determining the impact loads due to cavitation bubble collapse. The idea behind this is that any simulation of cavitation erosion requires firstly a thorough knowledge of the loads applied to the material. The difficulty is that bubble collapse generates very special loading conditions. They are impulsive loadings resulting from shocks of small duration (microseconds), extreme amplitude (gigapascals) and that are very localized in space (micrometres). In this paper, a method will be proposed to determine the spectrum of impact loads from joint pitting tests and finite-element modelling (FEM).

The second step consists in simulating the response of the material to repetitive impact loads. The objective is to predict the behaviour of the material as a function of the exposure time, i.e. for an increasing number of impact loads. The model should be able to account for the various stages of the erosion process. This includes the incubation period where no significant mass loss is observed and the more advanced stages of erosion where material is progressively removed by collapsing bubbles. For such a simulation, the impact loads identified in the first step should be applied repetitively and randomly on the material surface.

This computational step is based on the stress–strain relationship of the material that should include the effect of strain rate, which is quite high in cavitation because of the small duration of the impact loads. It requires an appropriate modelling of the hardening mechanism in order to account for fatigue that is observed during the incubation period. It also requires modelling the damage process, including failure that leads to mass loss.

The approach proposed here is a one-way simulation as the impact loads are supposed to be unchanged by the erosion of the solid wall. In other words, it is assumed that the change in the wall geometry does not significantly affect the cavitating flow and the impact loads. The model is then supposed to reach a steady-state regime of erosion characterized by a constant erosion rate beyond a given exposure time (figure 2). It is not able to predict the longer term behaviour that may be characterized by a more complicated time evolution of mass loss such as an attenuation period [5] which is generally interpreted in terms of a reduction in flow aggressiveness because of a modification of the flow due to the change in wall geometry.

It may be, however, necessary that the impact loads determined in the first step include the damping effect due to the fluid–structure interaction. This is particularly important in the case of soft materials such as polymeric coatings that may undergo a significant deformation during bubble collapse, resulting in a reduction in the amplitude of impact loads in comparison with a perfectly rigid wall. A simple model is briefly described in this paper to account for the fluid–structure interaction.

Various laboratory devices have been developed in order to investigate cavitation erosion such as ultrasonic horns [6,7], cavitating liquid jets [6–10], rotating discs [11] and cavitation tunnels [2,12]. The approach presented in this paper is supported by experimental data (in particular pitting data) obtained in a cavitation erosion tunnel using a radial divergent test section [1–4].

The two steps mentioned above are addressed in a separate section. Section 2 is devoted to the determination of cavitation impact loads using pitting tests combined with an inverse FEM simulation. Section 3 is devoted to an analysis of the response of the material to repetitive impact loads with a special emphasis on the most relevant material properties to be considered for damage prediction. It should be emphasized that the FEM simulations of fatigue and damage presented in §3 are very preliminary results that need further research developments before the proposed method is effective for real applications.

## 2. Cavitation impact loads

In this section, a method is proposed to estimate the loading conditions due to bubble collapse. It is based on pitting tests that have long been recognized as a suitable way to characterize flow aggressiveness [13–15].

Pitting tests are short duration tests conducted on ductile metallic materials. The sample surface is mirror polished prior to exposure to cavitation in order to reveal the tiny plastic deformations or pits due to bubble collapse. Figure 3 presents a typical example of the surface of a stainless steel sample exposed to cavitation for a relatively short time. The exposure time should be adjusted in order to have a sufficient number of pits to enable a statistical analysis, on the one hand, and to limit pit overlapping that would bias the analysis, on the other hand.

The issue addressed here is how to derive quantitative information on the hydrodynamic impact loads responsible for every pit in terms of applied pressure in particular. The idea is to use FEM of single cavitation impacts and adjust the parameters of the applied pressure field in order to reproduce numerically each pit observed experimentally.

In order that the method remains tractable, simple loading conditions have been used in this paper. It is assumed that the applied load due to the collapse of a cavitation bubble has a Gaussian shape given by the following equation:
2.1
where *r* is the radial distance along the wall from the bubble axis. This equation gives the space distribution of the hydrodynamic impact load in terms of amplitude *σ*_{H} and radial extent *d*_{H}. The load given by equation (2.1) is the normal component of the applied stress tensor.

The Gaussian law used here is a simplified model of the actual loading conditions. It is well known that several phenomena occur during the collapse of a bubble near a wall including the development of a microjet and shockwaves [1]. As a result, the actual loading is much more complicated than that described by equation (2.1). In addition, unsteadiness in loading and unloading is ignored in equation (2.1) as no time dependence is included. Then, there is obviously room for improving the present approach. Even though the time variation of loading is ignored, the high strain rate encountered in cavitation is here taken into account to some extent by using, in FEM computations, stress–strain relationships obtained at high strain rate (1000 s^{−1}) as opposed to quasi-steady stress–strain relationships.

When the amplitude of the impact load is large enough compared with the material yield stress, a permanent indentation is formed on the material surface. In the present approach, the resulting pit is characterized by two geometrical parameters, namely its residual maximum depth *h*_{p} and residual diameter *d*_{p} measured at mid-depth (figure 4*b*).

The method proposed in this paper to determine the two hydrodynamic parameters *σ*_{H} and *d*_{H} consists in simulating the response of the material to the Gaussian impact load and computing numerically the geometrical characteristics *h*_{p} and *d*_{p} of the pit using an FEM approach. It is shown below that an inverse computation can be developed in order to determine the two hydrodynamic parameters *σ*_{H} and *d*_{H} from the two measured parameters *h*_{p} and *d*_{p}.

When the Gaussian load is numerically applied onto the surface of an elastic–plastic material, the result after elastic recovery is a permanent deformation as shown in figure 5. A large number of computations were made for values of *σ*_{H} between 1.5 and 3 GPa and values of *d*_{H} between 10 and 240 µm. For each impact load, the depth *h*_{p} and diameter *d*_{p} of the resulting pit have been computed. A summary of the results is presented in figure 6.

Figure 6 shows that a one-to-one correspondence between hydrodynamic load parameters (*σ*_{H}, *d*_{H}) and pit geometrical parameters (*h*_{p}/*d*_{p}, *d*_{p}) exists. The interest of considering the pit shape factor *h*_{p}/*d*_{p} rather than pit depth *h*_{p} is that the isolines of load amplitude *σ*_{H} are almost horizontal. This means that the pit shape factor is strongly correlated to *σ*_{H}. In contrast, the isolines of constant impact load diameter *d*_{H} are roughly vertical, which means that pit diameter depends primarily upon impact load diameter and that the amplitude of impact load has a minor effect on pit diameter.

It is concluded from figure 6 that an inverse method can be developed from a series of FEM computations in order to determine the hydrodynamic characteristics of the impact load for each measured pit. The method is presented in detail in [16]. A semi-infinite axisymmetric model is used. The mesh is parametrized by the hydrodynamic parameter *d*_{H}. Minimum element size is 0.123 µm for *d*_{H} = 20 µm. Infinite elements are introduced in order to meet the zero stress/zero strain condition at infinity. The mesh has been validated by applying the reference Hertzian loading on a purely elastic semi-infinite material for which the analytical solution is available. Computations are made using Abaqus standard. Although the numerical model is quasi-steady, it should be emphasized that the effect of strain rate on the material behaviour is, to some extent, included by considering a stress–strain relationship at high strain rate (see §3). Typical results of an inverse calculation are shown in figures 7 and 8.

Figure 7 shows the evolution of the rate of impact loads as a function of their amplitude *σ*_{H}. Two curves are plotted that correspond to two different ranges in impact diameter *d*_{H}. The vertical axis is the rate of impact loads per unit surface area exposed to cavitation. In addition, it has been divided by the two bandwidths in *σ*_{H} and *d*_{H} to get density functions insensitive to bandwidths. This explains the presence of megapascals and micrometres in the unit.

The data plotted in figure 7 have been obtained from inverse computations using three pitting tests on three different materials: an aluminium alloy Al 7075, a nickel aluminium bronze (NAB) and a duplex stainless steel alloy A2205. As shown in figure 7, each material allows us to investigate a different range in stress amplitude. Soft materials (such as the aluminium alloy) permit the determination of the impact load spectrum for relatively low stress amplitude whereas harder materials (such as stainless steels) permit investigating higher amplitudes. The benefit of using several materials is to broaden the range in stress amplitude. It is clear that the material used for pitting tests must not be harder than the one for which the prediction of long-term erosion is desired. If it is harder, a fraction of the impact loads of small amplitude will not be detected by pitting tests whereas they actually contribute to damage.

It is interesting to observe in figure 7 that the data corresponding to the three materials tend to define a unique curve. This shows that the obtained spectrum is material independent and thus characterizes only the cavitating flow. In other words, the materials are used as sensors that provide a measurement of hydrodynamic parameters. The big advantage of using the material as a sensor in comparison with conventional pressure sensors is that it provides an estimate of the size of the loaded region that is definitely out of reach of pressure sensors. The method obviously requires the material properties to be determined accurately and particularly the stress–strain curve. This point is addressed in the next section.

As the impact loads are here characterized by two parameters, namely the impact stress *σ*_{H} and impact diameter *d*_{H}, it should be interesting to represent the spectrum of rate of impact loads versus both parameters using a three-dimensional picture, as shown in figure 8. The surface thus obtained is considered in this paper as a measure of the cavitation intensity or flow aggressiveness.

Other parameters may be considered in addition to *σ*_{H} and *d*_{H} for a more accurate description of each impact such as the duration of impact. If this were the case, an additional axis could be added in figure 8 and flow aggressiveness would be characterized by a hypersurface. The larger the number of parameters the more accurate the model of flow aggressiveness is.

The next step is to simulate the response of the material to the so-defined flow aggressiveness in order to predict cavitation damage.

## 3. Material response to cavitation impact loads

This section provides a rough outline of a potential approach for a numerical prediction of cavitation damage. It is based on an FEM method that aims at simulating the response of the material to the cavitation intensity defined in §2. An important issue is to determine the material properties in order to account for the phenomena involved in cavitation erosion. This includes the determination of the stress–strain relationship at high strain rate as well as the use of an appropriate hardening model for fatigue and an appropriate failure model for damage.

The FEM simulation of cavitation damage is still the subject of research and the results presented in this section are very preliminary results that need to be confirmed and improved before a consolidated approach for cavitation damage prediction is actually available. In particular, the numerical simulations presented below are not directly related to the response of the material to cavitation impact loads but to repetitive loads generated by a spherical indenter that is supposed to model to some extent a repetitive cavitation impact.

It should be emphasized that the fluid–structure interaction was not considered for determining the impact loads in §2. This is because pitting tests from which impact loads were derived are made using metallic materials that are generally quite stiff. As a result, the hydrodynamic impact loads are only weakly damped by the recoil effect of the wall. This is no longer the case for soft materials that may significantly reduce the impact stress in comparison with rigid walls. In this case, the fluid–structure interaction must be taken into account in order to update the impact loads. This is especially the case for soft coatings such as paints or polymeric coatings.

In a first approach, the importance of the fluid–structure interaction can be estimated using a simple one-dimensional model of the impact of a high-velocity liquid jet on an elastic material. The shock of the liquid on the material surface is supposed to generate two plane shock waves, one propagating forward into the solid and the other one backward in the liquid. The impact pressure is given by [5] 2.2

The numerator is the impact pressure of a liquid jet of velocity *V*_{l} on a perfectly rigid wall. It is derived from the common water hammer formula. The quantity *ρ*c is the acoustic impedance of the liquid (subscript l) or solid (subscript s). The acoustic impedance of the solid is related to Young's modulus *E* by

Equation (2.2) shows that the actual pressure on an elastic material is smaller than that on a rigid wall. According to this model, the damping effect is directly connected to the ratio of the acoustic impedances Figure 9 shows the influence of *α* on the non-dimensional pressure The perfectly rigid wall corresponds to an infinite value of *α*, whereas the compliant wall corresponds to *α* = 0, for which the impact pressure vanishes.

As an example, let us consider a coating made of ultra-high molecular weight polyethylene (UHMWPE) exposed to a cavitating water flow. In this particular case, the mechanical properties are as follows: The ratio of acoustic impedances is then As a result, the impact pressure given by equation (2.2) is only 56% of the impact pressure in the perfectly rigid case.

This simple model of fluid–structure interaction shows that the amplitude of impact loads may be quite significantly reduced in the case of a compliant wall. The same conclusion is obtained when considering a stiffer liquid such as cavitating mercury. In the case of a compliant wall or a stiff liquid, the impact stresses as determined in §2 must be updated to account for the damping effect due to fluid–structure interaction for the correct prediction of cavitation damage.

An essential input of FEM simulations is the stress–strain relationship of the material. For cavitation applications, it must include the effect of strain rate that generally results in an increase in flow stress for ductile materials. Figure 10 presents an example of the influence of strain rate on material flow stress for stainless steel A2205. Stress generally increases with strain rate in a logarithmic way. This behaviour is described by the Johnson–Cook equation:
2.3
where *K* is a constant that can be determined from the slope of the lines in figure 10 and is an arbitrary reference strain rate (e.g. ).

For cavitation impacts, the strain rate may reach high values, larger than several thousands s^{−1}. A typical order of magnitude of can be obtained by assuming that the impact generates a strain of the order of a few per cent (Δ*ɛ* ≈ 10^{−2}), which is a typical order of magnitude of strain for a pit, whereas its duration is of the order of a few microseconds (Δ*t* ≈ 10^{−6} s). This gives (e.g. [1]).

In order to characterize the material behaviour at such a high strain rate, a special experimental device should be used. The cloud of experimental data shown in figure 10 above about 100 s^{−1} were obtained from split Hopkinson pressure bar (SHPB) tests (figure 11*b*). All other points below 100 s^{−1} were obtained from conventional compression tests (figure 11*a*) that allow moderate values of strain rate to be investigated. More details on SHPB tests can be found in [17].

The next step after determining the stress–strain curve is to model the hardening mechanism. A specificity of cavitation erosion is the repetitive loading that induces fatigue. Fatigue is particularly obvious during the incubation period because cavitation bubbles collapse repetitively on the material surface without producing any mass loss. During the incubation period, cumulated plastic strain progressively increases within the material until damage occurs. The occurrence of damage determines the end of the incubation period and the onset of the acceleration period.

The conventional model of isotropic hardening is not appropriate to model the incubation period. This is explained in figure 12*a* where it is assumed, for simplicity, that the material has a bilinear behaviour and the impacts have all the same amplitude *σ*_{H}. After the first impact, the material is in state B and its yield stress has increased from *σ*_{Y} to *σ*_{H}. If a second impact with the same (or smaller) amplitude hits the surface at the same point, the material has a purely elastic behaviour and its representative point moves during loading and unloading along path BC. Thus, the second impact and all subsequent impacts have a purely reversible elastic effect but do not deposit additional energy into the material, thus making this hardening model inappropriate for cavitation erosion modelling.

On the other hand, kinematic hardening allows us to account for fatigue. The loading phase OAB for the first impact is exactly the same as for isotropic hardening. However, during unloading, the material follows path BCD, which comprises not only an elastic part BC but also an additional plastic deformation CD. The stress difference between points B and C is twice the initial yield stress *σ*_{Y} as the yield surface is translated but retains the same size in kinematic hardening; this is in contrast to isotropic hardening, for which it expands.

During the second impact, the material representative point describes the loop DEBCD that has a non-zero thickness. As a result, additional energy corresponding to the area of loop DEBCD is deposited into the material after each new impact. This is confirmed by figure 12*c*, which compares the evolution of work after successive impacts for both hardening models. In the case of isotropic hardening, work after unloading and elastic recovery remains constant whatever the number of impacts. In the case of kinematic hardening, work increases step by step after each impact. This behaviour is consistent with fatigue that develops during the incubation period.

The mechanism of fatigue is illustrated in figure 13 in the case of repetitive impacts that are here modelled in a very simplified manner by quasi-steady repetitive indentation tests using a spherical indenter. The progressive increase in plastic strain in the material is clearly visible. The point is now to model the onset of damage in order to be able to predict the longer term behaviour.

For an elastic–plastic material, damage is characterized by both the degradation of stiffness and the degradation of elasticity (figure 14*a*). A common approach to account for both effects consists in introducing a scalar damage variable *D*. Undamaged material corresponds to *D* = 0, whereas *D* = 1 means complete failure with the occurrence of fracture. FEM codes generally correlate the damage variable *D* with the cumulative plastic strain Several models are available. A simple version is given by a linear correlation between the two as shown in figure 14*b*. Two critical values of the cumulative plastic strain *p* are then introduced. The first critical value *ɛ*_{frac} indicates damage initiation. It corresponds to the point where the damage variable starts increasing. The second critical value *ɛ*_{frac} + *ɛ*_{fail} corresponds to failure, i.e. *D* = 1.

In order to demonstrate, at least qualitatively, the relevance of this model for cavitation damage prediction, the model has been used in an FEM simulation of successive impacts using a spherical nanoindenter to model the impact. Results are presented in figure 15. During the first nine impacts, fatigue takes place in the material and strain progressively increases. This initial stage would correspond to the incubation period in cavitation erosion.

After impact 10, the failure criterion has been reached in one mesh element (in white in figure 15) where cumulative deformation is maximum. Let us note that this element is within the material, below the surface. This is the very beginning of the acceleration period. After one more impact (impact 11), several new elements are damaged, and, after impact 13, a piece of material becomes surrounded by elements that have all lost their load-carrying capacity so that this piece of material becomes detached as a whole.

This simplified simulation shows that, by introducing a suitable hardening mechanism together with a damage model, one should be able to model the fatigue mechanism that takes place during the incubation period as well as the more advanced stage of erosion characterized by material removal and mass loss. In particular, it is expected that the usual mass loss curve that characterizes the kinetics of cavitation erosion damage could be predicted from such a numerical approach.

It should be emphasized that the results presented in this paper were obtained from quasi-steady computations. They need to be revisited in the framework of dynamic computations including a model of strain rate sensitivity such as the Johnson–Cook model. Under such realistic conditions, some strain accumulation may be expected even under isotropic hardening. The reason is that strain rate may change from one impact to another. Strain rate is expected to be maximum for the first impact where the increase of strain Δ*ɛ* is maximum. For the second impact, the material will deform less. The increase of strain will be less and, assuming that the impact duration Δ*t* is the same, strain rate will be less and consequently flow stress will be less. As a result, the material appears softer and can undergo an additional deformation during the second impact. It can be expected that a saturation mechanism will take place after a number of impacts, provided that the material is not damaged.

## 4. Conclusion

Although much remains to be done to apply this approach to cavitation erosion and validate it, the preliminary results presented in this paper tend to prove that prediction of cavitation erosion damage is now within the reach of numerical simulation.

The proposed method requires first estimating the distribution of impact loads due to cavitation bubble collapse. This could be done from pitting tests carried out on ductile materials. It was shown in §2 that an inverse FEM computation could be developed in order to determine both the stress amplitude and the radial extent of the hydrodynamic impact responsible for each pit.

Once the impact load spectrum is determined, the method consists in simulating the response of the material to these impacts of various amplitudes and sizes randomly hitting the wall. This could be achieved by FEM of the material response. By using an appropriate hardening model as well as a damage model, it can be expected that the incubation period as well as the evolution of mass loss as a function of exposure time could be numerically predicted.

The approach requires a precise characterization of material properties that is suited to the physics of cavitation erosion. This includes characterizing the material response at high strain rate. It is also essential to characterize the fatigue mechanism that takes place in the material during the repetitive loading by collapsing bubbles, which requires using a hardening law that accounts for a kinematic component. In addition, an appropriate damage model should be used in order to predict material removal and mass loss. All these models introduce material properties such as critical strains for damage initiation and material failure that need to be determined.

Such an approach is *a priori* applicable to any material provided the mechanical properties are available. In particular, it could be applied to compliant materials, coatings and multi-layered materials of interest for cavitation erosion mitigation. Such a numerical approach could be very helpful for optimization purposes.

The method presented here is not a purely numerical method as it requires pitting tests for estimating the impact loads. It can however be expected that this experimental step could be replaced in the future by a numerical step. The idea is to predict the impact loads from computational fluid dynamics (CFD). Recent results (e.g. [18,19]) are quite promising regarding the numerical prediction of impact load spectra due to the collapse of cavitating structures. As a long-term vision, the CFD could be directly coupled to the FEM in order to better account for the fluid–structure interaction and have a fully numerical approach to cavitation damage.

## Competing interests

We declare we have no competing interests.

## Funding

This research was conducted under a NICOP project funded by the Office of Naval Research.

## Acknowledgements

The authors wish to thank Dr Ki-Han Kim from the Office of Naval Research (ONR) and Dr Woei-Min Lin from the Office of Naval Research Global (ONRG) who supported this work. The authors are very grateful to Prof. Nicolas Ranc, Laboratoire Procédés et Ingénierie en Mécanique et Matériaux (PIMM), Ecole Nationale Supérieure d'Arts et Métiers de Paris, for his help in conducting Hopkinson tests at high strain rate. The authors are also very grateful to Mr Simon Garcin who made the FEM simulations presented in §3 in the framework of an internship.

## Footnotes

One contribution of 13 to a theme issue ‘Amazing (cavitation) bubbles: great potentials and challenges’.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.