# A Cavitation Erosion Model for Ductile Materials

Nabila A. Berchiche, Jean-Pierre Franc, Jean-Marie Michel

► To cite this version:

Nabila A. Berchiche, Jean-Pierre Franc, Jean-Marie Michel. A Cavitation Erosion Model for Ductile Materials. Journal of Fluids Engineering, 2002, 124 (3), pp.601-606. 10.1115/1.1486474. hal-00212305

HAL Id: hal-00212305

https://hal.science/hal-00212305v1

Submitted on 21 Jan 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A Cavitation Erosion Model for Ductile Materials

An analytical model is proposed for the prediction of cavitation erosion of ductile materials. It is based upon a physical analysis of the work-hardening process due to the successive bubble collapses. The material is characterized by its classical stress-strain relationship and its metallurgical behavior is analyzed from microhardness measurements on cross sections of eroded samples. The flow aggressiveness is determined from pitting tests, using the material properties to go back to the impact loads. The histogram of impact loads is applied numerically a large number of times on the material surface and the evolution of the mass loss with the exposure time is computed. The approach is supported by experimental tests.

N. Berchiche

Department of Naval Architecture and Ocean Engineering, SE-412 96 Göteborg, Sweden
e-mail: bena@na.chalmers.se

J. P. Franc

e-mail: Jean-Pierre.Franc@hmg.inpg.fr

J. M. Michel

e-mail: Jean-Marie.Michel@hmg.inpg.fr

Laboratoire des Ecoulements Géophysiques et Industriels,

BP 53, 38041 GRENBOLE Cedex 9, France

## 1 Introduction

A lot of work has been done in the field of cavitation erosion. Roughly speaking, the techniques of prediction of cavitation erosion can be classified into three main categories:

* empirical correlations with material properties or with electrochemical or noise measurements

* simulation techniques using special test devices to reproduce a given aggressiveness in an accelerated way

* analytical methods.

Detailed information on the two first categories of techniques can be found in review papers such as Hammitt [1] Karimi and Martin [2], Franc and Michel [3]. The present paper is especially devoted to analytical techniques whose objective is to predict cavitation erosion without model tests or at least with a limited request to experiments as firstly imagined by Kato et al. [4]. Such techniques are still in development and represent a real challenge to research workers in cavitation erosion for the next years. They require extensive research efforts and, therefore, it is not surprising that the first attempts, including the present work, do not lead to models fully operational and satisfactory. However, such a new research direction is worthwhile being explored, in our opinion, as it might open a new field of techniques of prediction of cavitation erosion.

The present paper is a contribution to this subject. It presents a model of prediction of the erosion damage applicable to ductile materials only. Other limitations of the model will be pointed out along the presentation. The originality of this work lies in the fact that the proposed model is fully predictive and involves no parameters to be adjusted on the basis of experimental data. It is based upon the original work of Karimi and Leo [5]. Contrary to Karimi and Leo in which each pit is treated as a whole, the present approach is much more local. A regular mesh, with a characteristic size of the order of a few micrometers only, is defined on the material surface, so that each pit is described by a significant number of cells. This technique allows to compute, with a space resolution which can be refined as much as required, the strain field on the material surface and inside, from which the erosion rate is deduced. The main input is the surface distribution of impact loads which is determined from pitting tests.

The characteristic of a ductile material exposed to cavitation is to be progressively hardened by the successive collapses. The

work-hardening process is here characterized by the thickness of the hardened layers together with the shape of the strain profile inside the material.

The first step of the proposed model consists in the quantification of the hydrodynamic aggressiveness of the cavitating flow. This is done from classical pitting tests, each pit being characterized by its diameter and its depth. The surface distribution of the impact load responsible for each pit is deduced from this couple of data, using the strain profile and the stress-strain relationship of the material. The flow aggressiveness is finally characterized by a distribution of impact loads.

In a second step, this distribution derived from short duration tests is numerically applied a large number of times on the material surface. The present model computes the mass-loss as a function of the exposure time. To support this approach, a few pitting and mass loss tests were conducted on an experimental device which produces cavitation erosion from the collapse of a cavitating vortex (Dominguez-Cortazar et al. [6] Filali and Michel [7], Filali et al. [8]).

## 2 Presentation of the Model

The principle of the model is presented in the simplified case of a perfectly reproducible impact (Fig. 1). Let us consider a given point of the material surface on which a stress $\sigma_1$ is applied due to a bubble collapse. The material is initially supposed to be virgin, i.e., the strain is zero everywhere, on the surface and inside the material. If $\sigma_1$ is lower than the elastic limit $\sigma_e$, the material is supposed to return to its original state after unloading. Hence the impact loads below the elastic limit have no effect. In particular, fatigue mechanisms are not taken into account. Thus, the present model is applicable to sufficiently aggressive flows, which present a substantial number of impacts beyond the elastic limit.

After the first impact, the strain on the surface of the material has become $\varepsilon_1$, which is deduced from the stress-strain relationship of the material (Eq. (6)). The distribution of strain inside the material is supposed to be given, for $x \le l$, by the following empirical relation (Fig. 1):

$$ \varepsilon(x) = \varepsilon_s \left(1 - \frac{x}{l}\right)^{\theta} \quad (1) $$

where $\varepsilon_s$ is the surface strain at the point of impact, $l$ the depth of the hardened layer, $\theta$ the shape factor of the strain profile and $\varepsilon(x)$ the strain at the distance $x$ from the surface. After the first impact, we have $\varepsilon_s = \varepsilon_1$ and $l=l_1$. The energy absorbed by the material is the shaded area.

A second collapse of exactly the same amplitude $\sigma_1$ is supposed to occur at exactly the same point. The surface strain will be increased up to a certain value $\varepsilon_2$, which is determined from the
<figure><img src="image_2.png" /><figcaption>Fig. 1 Principle of the model. The stress-strain relationship, together with the strain profiles inside the material are presented as a function of the exposure time.</figcaption></figure>

conservation of the impact energy. We suppose that the same energy (indicated by the shaded area on Fig. 1) is absorbed by the material, but this time, the material is no longer virgin. The conservation of the impact energy consists in writing that the area under the $\sigma(\varepsilon)$ curve between points ① and ② is the same as the one between points ① and ③. This condition determines the surface strain $\varepsilon_2$. The strain profile is still given by Eq. (1) with $\varepsilon_s = \varepsilon_2$ and a greater depth $l=l_2$ of the hardened layer which will be determined later, from Eq. (4).

Suppose that a third, still identical, impact is applied. In the particular case of Fig. 1, the conservation of energy leads to a surface strain beyond rupture denoted by $R$. In this case, the $\sigma(\varepsilon)$ curve is extrapolated to get the virtual surface strain $\varepsilon_3$. The strain profile is still given by Eq. (1). As the material cannot withstand a strain greater than its rupture strain, it is supposed that the thickness $a_3$ is removed, which corresponds to the domain in which $\varepsilon$ is greater than $\varepsilon_r$. It can easily be shown that the thickness of the eroded layer is given by:

$$a=L\left[\left(\frac{\varepsilon_s}{\varepsilon_r}\right)^{1/\theta}-1\right] \quad (2)$$

$\varepsilon_s$ is the virtual surface strain (greater than $\varepsilon_r$) and $L$ the maximum thickness of the hardened layer, corresponding to the rupture strain $\varepsilon_r$.

If a fourth identical impact is applied, the surface strain is increased from $\varepsilon_r$ to $\varepsilon_4$. The virtual strain $\varepsilon_4$ is still determined by the conservation of energy, the area below the $\sigma(\varepsilon)$ curve between

points R and ④ being equal to the original impact energy (shaded area). The thickness of the eroded layer is still given by Eq. (2).

As soon as the surface strain has reached the rupture strain $\varepsilon_r$, hardening is maximum. The strain profile inside the material remains unchanged and given by:

$$\varepsilon(x)=\varepsilon_r\left(1-\frac{x}{L}\right)^{\theta} \quad (3)$$

$L$ appears as the maximum thickness of the hardened layer. For partial hardening, leading to a surface strain $\varepsilon_s<\varepsilon_r$, it can easily be shown that the thickness $l$ of the hardened layer is smaller than $L$ and given by:

$$l=L\left(\frac{\varepsilon_s}{\varepsilon_r}\right)^{1/\theta} \quad (4)$$

This equation results from the assumption that the strain profile for partial hardening (Eq. (1)) corresponds to a truncated part of the complete profile described by Eq. (3).

Above, we have examined in detail the simplified case of a perfectly reproducible impact. The principle of the method remains applicable to the more general case of variable loading occurring in the real process of cavitation erosion. The main difference is that the energy absorbed by the material does not remain constant and has to be evaluated for each impact. However, the method is still based upon the principle of energy conservation.

Another difference with respect to the above simplified presentation lies in the evaluation of the energy. In the computation, the energy is not limited to the one absorbed by the surface of the material (as we could believe from Fig. 1), but it corresponds to the total energy actually absorbed by all the hardened layers inside the material. The energy absorbed by the material per unit surface area between an initial virgin state and a state characterized by the strain profile given by Eq. (1) can be calculated as follows:

$$W=\int_{0}^{l_1}\left[\int_{0}^{\varepsilon}\sigma d\varepsilon\right]dx=\frac{\sigma_e\varepsilon_1l_1}{\theta+1}+\frac{K\varepsilon_1^{n+1}l_1}{(n+1)(n\theta+\theta+1)} \quad (5)$$

The quantities $\sigma_e$, $K$, and $n$ are characteristics of the stress-strain relationship and are defined in Eq. (6).

The model is purely one-dimensional. It is supposed that there is no interaction between two neighboring points situated on the material surface or at the same distance from it. The limitations of this assumption are not yet fully understood. In the following computations, a regular surface mesh of 440×440 points is defined on the material surface. The distance between two consecutive points was chosen equal to 5 µm. This value appeared to be a good compromise between the accuracy of the computation and the CPU time. In particular, the mesh size must be small enough to allow a good description of the smallest pits. In the present case, the pits with a diameter smaller than 20µm where not considered. Hence, the smallest pits are defined by a mesh of about 5×5 points.

In conclusion, the present model consists in computing, at each time step, the distribution of strain limited to the material surface, from which all other data can be deduced, including mass-loss and strain field inside the material.

## 3 Material Characteristics

Two series of tests are used to characterize the material. The first one is the classical tensile test which allows to determine the stress-strain relationship. For ductile materials, it is correctly represented by a Ludwig type equation:

$$\sigma=\sigma_e+K\varepsilon^n \quad (6)$$

Because of the high value of the rupture strain for ductile materials, the elastic part of the curve can be considered as almost vertical and the elastic energy can be neglected. In the case of stain-
less steel 316L considered here, we obtained: $\sigma_e = 400$ MPa, $\sigma_r = 1020$ MPa, $n = 0.5$ and $K = 900$ MPa. The rupture strain is $\epsilon_r = 47\%$.

It has to be emphasized that these data were obtained from classical quasi-steady tests, with a very small strain rate of the order of $10^{-4} s^{-1}$. They are used here without any modification for the analysis of cavitation erosion which is known to be characterized by an incomparable higher strain rate of the order of $10^4-10^5 s^{-1}$. The influence of the strain rate is difficult to take into account. A possible approach could consist in artificially increasing the elastic and rupture limits of the material, but this was not done in the present work due to the lack of data.

The two main metallurgical parameters introduced in the model, the maximum depth of the hardened layer $L$ and the shape factor of the strain profile $\theta$ are determined from micro-hardness measurements on cross-sections of an eroded target. We obtained the following values $L = 200 \ \mu m$ and $\theta = 5.0$.

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2 Photograph of the impacted zone on Stainless Steel 316 after a pitting test of 30 shots on the Cavermod. (a) Nomarski interferometric technique; (b) Mirau interferometric technique ($100 \ \mu m$ corresponds to 0.8 cm and 0.6 cm, respectively).</figcaption></figure>

## 4 Pitting Tests

Two pitting tests have been carried out on stainless steel 316L after 30 shots of the experimental device. The number of shots was selected to get a large enough number of pits without significant overlapping. Figure 2 presents two photographs of the same eroded surface with two different observation techniques. Figure 2(b) is obtained on a metallurgical microscope using a Mirau interferometric technique (Belahadji et al. [9]). The main advantage

<figure><img src="image_5.png" /><figcaption>Fig. 3 Histograms of pit number (a) and deformed volume (b) corresponding to the pitting test presented in Fig. 2. The deformed volume is defined as the volume of the pits below the original surface.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4 Pitting test reproduced by the model. The difference in height between two consecutive curves is .46 $\mu m$</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Fig. 5 Step-by-step description of the entire predicting process</figcaption></figure>

of this technique is to allow the estimation of the depth of each pit as the distance between two black fringes corresponds to half a wavelength (0.273 µm).

A large pit is observed in the center, surrounded by smaller pits distributed randomly on the impacted zone. The diameter of this zone is approximately 1000 µm. It depends on the size of the cavitating vortex and the resistance of the material (Filali and Michel [7]).

The analysis of a pitting test consists in determining, for each pit, the coordinates of its center, its maximum depth and its diameter ($2r_e$). The deformed volume is here estimated assuming a simplified conical shape for each pit. This assumption, which consists in supposing that the fringes in Fig. 2(b) are circular and equidistant, proved to be a reasonable approximation. More accurate techniques for the determination of the complete 3D-shape of the pits have been developed (see e.g., Belahadji et al. [9]) but were not available for this work.

Figure 3(a) shows the histogram of pits number versus the pit diameter, resulting from the analysis of the photograph of Fig. 2. Although the number of small pits is very large, they have a relatively small contribution to the deformed volume, which results mainly from the larger pit as shown on Fig. 3(b).

The aggressiveness of the cavitating flow in terms of applied stresses is deduced from the analysis of the pits produced on the material surface during the early stage of erosion. Let us consider a pit of maximum depth $h_{max}$. By integration of the strain profile (Eq. (1)), we obtain the following relation between pit depth and surface strain $\epsilon_s$:

$$h_{max} = \int_{0}^{l} \epsilon(x) \cdot dx = \frac{l\epsilon_s}{\theta + 1} \quad (7)$$

The measurement of pit depth allows to determine the surface strain, and as a result, the original stress $\sigma_{max}$ by the use of the stress-strain relationship (Eq. (6)).

Once the maximum load at the center of the pit is known, the radial distribution is determined by assuming that it follows a gaussian law:

$$\sigma = \sigma_{max} \left[ \frac{\sigma_{max}}{\sigma_e} \right]^{-r^2/r_e^2} \quad (8)$$

where $r_e$ is the measured pit radius. This assumption would not be necessary in case of a complete 3D measurement of the pit shape. The Mirau interferometric technique is considered to give a good

estimate of the size of the plastic zone. This equation takes into account that the stress is equal to the elastic limit $\sigma_e$ at the limit $r=r_e$ of the plastic zone.

By considering all the pits which were identified on Fig. 2(b) it is possible to go back to the distribution of stresses. The reproduction of this distribution by the numerical model allows to reconstruct the image of the surface after the pitting test. Figure 4 presents the results of the “numerical” pitting test. Each pit in Figs. 2(a) and 2(b) can be identified in Fig. 4. The main difference is the perfectly circular shape of each pit in the model.

## 5 Computation of the Erosion Rate and Comparison With Experiment

Once the distribution of impact loads resulting from a reference pitting test is determined, it is applied randomly over the exposed area a large number of times until mass loss occurs. For each pit, only the coordinates of its center are chosen randomly, whereas the impact load and the pit diameter are kept unchanged. A step-by-step description of the entire prediction process is given in Fig. 5. Results of the prediction are presented in Figs. 6, 7, and 8. Let us notice that the number of shots is here equivalent to a classical exposure time as encountered in the case of continuous cavitation.

<figure><img src="image_8.png" /><figcaption>Fig. 6 Calculated and measured erosion rates versus the number of shots (stainless steel 316L). The two calculated erosion rates are based upon two different pitting tests used to characterize the flow aggressiveness.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7 Computed evolution of the strain on the surface of the material (stainless steel 316L). The blue color corresponds to zero strain, whereas the red color corresponds to the rupture strain in %.</figcaption></figure>

Figure 6 presents the two calculated erosion rates obtained from two different pitting tests. Three main stages are predicted in the erosion process: incubation, acceleration and steady state erosion. The mass-loss fluctuations are important during the period of acceleration and vanish in the steady state. This effect is typical of a random distribution of impacts loads. In the case of a repetitive single impact, these fluctuations do not exist (Berchiche [10]). They are due to the fact that the same impact can lead to a very different mass loss according to the degree of hardening. Mass loss is minimum if the impact falls on a virgin surface and maximum if hardening is completed.

Although the pitting tests were carried out under the same experimental conditions, we observe a variation of 30% between the two calculated erosion rates (Fig. 6). Hence, the accuracy of the long-term prediction depends strongly upon the pitting test from which the flow aggressiveness is characterized. In order to limit the sensitivity of the prediction to the pitting test, we suggest to use several pitting tests for the determination of the cavitating flow aggressiveness. This observation is not surprising in so far as pit size and pit load were kept constant throughout the modeling procedure. It could be envisaged to use a more complicated model, by considering statistical laws for the distribution of size

<figure><img src="image_10.png" /><figcaption>Fig. 8 Computed evolution of the strain field on a cross section of the material and of the shape of the eroded surface</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9 Cross-section of an eroded sample after 14400 shots (SS316L). The marks for microhardness measurements are hardly visible.</figcaption></figure>

and amplitude of the impact loads. However, such an approach would require a further analysis to determine precisely these laws and know if they can be considered as independent ones or not.

The evolution of the surface strain as a function of the number of shots is given in Fig. 7. During the acceleration period, the fraction of the surface which is fully hardened (in red) progressively increases. When the whole exposed area is hardened, the erosion rate becomes constant: it is the steady state period.

Figure 8 presents the predicted evolution of strain on a cross section, together with the evolution of the shape of the material surface. During the incubation period, the pit depth results only of plastic deformation. Once the rupture strain is reached on the surface, mass loss occurs.

In order to validate the proposed model, mass loss tests have been carried out on stainless steel 316L in the same experimental conditions than those which were chosen for the determination of hydrodynamic solicitations (Fig. 9). The sample was weighted after each series of 900 shots. The calculated and experimental erosion rates are compared in Fig. 6.

The order of magnitude of the predicted erosion rate in the final steady stage of erosion appears to be in reasonable agreement with the experiment, whereas the duration of the incubation period is significantly underestimated. Concerning the total mass loss after 14400 shots, the experimental value is 2.2 mg whereas the two cumulative mass losses predicted by the model using the two different pitting tests are 2.4 and 3.4 mg. Once more, the order of magnitude of the computed mass loss is consistent with experiments. It is clear that further comparisons to experiments are required for the improvement and the validation of the model. However, the present results are considered as satisfactory in so far as the whole model is fully predictive. It does not involve any adjustable parameter. All the data required for the prediction are determined in a unique way from the material parameters which, in their turn, are determined from classical mechanical or metallurgical tests.

Finally, let us mention that the sensitivity of the prediction to the metallurgical and mechanical parameters of the material was studied. A material characterized by a large value of the shape factor and a small thickness of the hardened layer is more resistant to cavitation erosion. In addition, the duration of the incubation period is independent of the thickness of the hardened layer, but decreases when the shape factor increases.

## 6 Perspectives

Although the present model is fully predictive, we must be aware that several assumptions or shortcuts were necessary to complete the modelling. Among the most critical ones, we can mention the influence of the strain rate which was ignored. In addition, the material was characterized from tensile tests whereas

the actual solicitation in cavitation erosion is a compression. The consequences of the one-dimensional nature of the model are also difficult to estimate.

In the future, it would be interesting to link this kind of model to a classical computation of the cavitating flow. Such a tool would allow a prediction of cavitation erosion based only on the flow geometry, the operating hydrodynamic conditions (pressure and flow velocity) and the mechanical and metallurgical properties of the material. The main steps would be the following.

1. The cavitating flow is computed using a Navier-Stokes solver completed by a cavitation model. Bubble models (see e.g., Kubota et al. [11]) are probably the most suitable for a further prediction of cavitation erosion.

2. The hydrodynamic aggressiveness is deduced from the former computation. In the cavitation model of Kubota et al. for example, the evolution of the bubble radius is determined from the resolution of a Rayleigh-Plesset equation. As a consequence, for each bubble, the interface velocity during the phase of collapse is computed. It can be considered as the key parameter for the estimation of the aggressiveness of the collapsing bubble. Although this procedure has still to be cleared up, it seems physically realistic.

3. The last step consists in modelling the material response and computing the erosion rate using a model as the one presented in this work.

## Nomenclature

a = thickness of the eroded layer

K = constant in the stress-strain relationship (Eq. (6))

l, L = depth of hardened layers

n = exponent in the stress-strain relationship (Eq. (6))

r = radius

x = distance from the surface

ε = strain

σ = stress

θ = metallurgical shape factor (Eq. (1))

## Subscripts

e = elastic

r = rupture

s = surface

## References

[1] Hammitt, F. G., 1979, "Cavitation Erosion: The State of The Art And Predicting Capability," Appl. Mech. Rev., **32**(6), pp. 665-675.

[2] Karimi, A., and Martin, J. L., 1986, "Cavitation Erosion of Materials," Int. Met. Rev., **31**(1), pp. 1-26.

[3] Franc, J. P., and Michel, J. M., 1997, "Cavitation Erosion Research in France: The State of The Art," Journal of Marine Science and Technology, **2**, pp. 233-244.

[4] Kato, H., Konno, A., Maeda, M., and Yamaguchi, H., 1995, "Possibility of Quantitative Prediction of Cavitation Erosion Without Model Test," ASME J. Fluids Eng., **118**, pp. 582-588.

[5] Karimi, A., and Leo, W. R., 1987, "Phenomenological Model for Cavitation Rate Computation," Mater. Sci. Eng., **95**, pp. 1-14.

[6] Dominguez-Cortazar, M., Franc, J. P., and Michel, J. M., 1997, "The Erosive Axial Collapse of a Cavitating Vortex: An Experimental Study," ASME J. Fluids Eng., **119**(3), pp. 686-691.

[7] Filali, E. G., and Michel, J. M., 1999, "The Cavermod Device: Hydrodynamic Aspects And Erosion Tests," ASME J. Fluids Eng., **121**, pp. 305-311.

[8] Filali, E. G., Michel, J. M., Hattori, S., and Fujikawa, S., 1999, "The Cavermod Device: Force Measurements," ASME J. Fluids Eng., **121**, pp. 312-317.

[9] Belahadji, B., Franc, J. P., and Michel, J. M., 1991, "A Statistical Analysis of Cavitation Erosion Pits," ASME J. Fluids Eng., **113**, pp. 700-706.

[10] Berchiche, N., 2000, "Erosion De Cavitation d'un Metal Ductile: Etude Experimentale Et Modélisation," PhD thesis Institut National Polytechnique de Grenoble (in French).

[11] Kubota, A., Kato, H., and Yamaguchi, H., 1992, "A New Modelling of Cavitating Flows: A Numerical Study of Unsteady Cavitation on a Hydrofoil Section," J. Fluid Mech., **240**, pp. 59-96.