Isothermal compressible flow solver for prediction of cavitation erosion

Cavitation erosion can be observed on hydraulic mechanical devices and has long been studied as a challenging research topic. In this paper, a practical method to predict cavitation erosion was developed. To simulate cavitating flows using computational fluid dynamics (CFD), an isothermal compressible cavitating flow solver was developed using the OpenFOAM toolkit. A cavitation erosion coefficient was proposed as a function of the cavitation and Reynolds numbers. From the simulations, the cavitation erosion coefficients were calculated using metamodeling methods. Finally, the developed practical method of cavitation erosion was applied to a marine propeller. As observed in experiments, cavitation erosion was predicted near the propeller tip.


Introduction
Cavitation can be found in lots of hydrodynamic machineries, such as marine propellers, pumps, nozzles, and turbines. Cavitation causes abrasion and erosion on a metal surface, and noise and vibration in systems. Among many harmful outcomes of cavitation, cavitation erosion, which degrades hydrodynamic performance and at the same time results in structural damage, is one of the most critical issues. The prediction and prevention method for cavitation erosion has been an important research topic for many decades and studied by many researchers.
Extensive research with experiments, which include problems from specimen test to single bubble bust, has been undertaken to understand the phenomenon of cavitation erosion, for example by Lush (1983), Grant and Lush (1987), Karimi and Avellan (1986), Okada, Iwai, and Awazu (1989), Richman and McNaughton (1990), Soyama, Kumano, and Saka (2001), Hattori and Nakao (2001), Soyama and Futakawa (2004), Franc (2009), Hattori andKishimoto (2008), Abouel-Kasem and Ahmed (2008), and Keil, Pelz, Cordes, and Ludwig (2011). Besides the experimental studies, many researchers have developed mathematical models to predict cavitation erosion. Kato, Konno, Maeda, and Yamaguchi (1996) proposed a quantitative prediction scenario of the impact force on the surface from cavitation, based on a model for an isolated implosion of a single bubble in a finite space. Fortes Patella, Challier, Reboud, and Archer (2001) proposed the relation between the pressure wave energy, CONTACT Sunho Park spark@kmou.ac.kr created by bubble collapses, and the pit volume, representing the solid damage. They applied the relation to three materials, that is, aluminum, copper, and stainless steel. Berchiche, Franc, and Michel (2002) related the stress-stain and the strain profile for ductile materials and proposed a cavitation erosion prediction model. Dular, Stoffel, and Širok (2006) performed cavitation erosion tests for simple-shaped hydrofoils and developed a theoretical model that could predict the pit extent. Hattori and Kishimoto (2008) predicted the erosion rate with the Vickers hardness. Szkodo (2005) suggested a mathematical model to describe the cavitation erosion resistance of materials based on Weibull's distribution. Van Terwisga (2009) reviewed physical mechanism and cavitation erosion models. Peters, Sagar, Lantermann, and el Moctar (2015) proposed an intensity coefficient to determine the erosion intensity using the critical and impact velocities. Koukouvinis, Gavaises, Li, and Wang (2016) simulated cavitating flow in a diesel injector using large eddy simulation (LES) and correlated cavitation erosion with pressure peaks. Mouvanal, Chatterjee, Bakshi, Burkhardt, and Mohr (2018) suggested a collapse detector algorithm and predicted cavitation erosion in fuel injectors using the unsteady Reynolds-averaged Navier-Stokes (URANS) approach. Computational studies on cavitation erosion have not considered the temperature and time variations.
The objective of the present study was to develop a practical method for the prediction of cavitation erosion.
For computational analyses, a cavitating flow solver that couples velocity, pressure, and phase compositions was developed using the OpenFOAM toolkit based on a finite volume method. The RANS approach was selected for the turbulence closure. The developed cavitating flow solver was validated by applying it to sheet and cloud cavitating flows (Park & Rhee, 2013.
The present study is organized as follows. The development of a practical method for cavitation erosion prediction is presented first, followed by the code development. The cavitation erosion coefficient is introduced and discussed. The developed method is then validated and applied. Finally, conclusions are provided.

Development of practical method
Once cavitation occurs, cavity bubbles grow and suddenly implode. The cavity implosion is followed by a powerful high speed micro-jet and a strong blast wave near the surface, which erodes it. To predict the erosion caused by cavity implosion, two physical quantities, impact and critical velocities, have to be considered. Here, the impact velocity indicates how fast the micro-jet and the blast wave impact the solid surface and the critical velocity is the threshold velocity for a solid surface to be deformed. In this study, the two velocities were defined as follows.

Impact velocity
When a single spherical bubble with incompressible flow in a quiescent fluid collapses near a solid surface, a jet reaches a high speed quite early in the collapse process. The jet velocity can be represented by the fluid density, the non-dimensional distance from the bubble center to the surface, and the difference of pressure between the pressure that would maintain the bubble at equilibrium and the far-field pressure (Brennen, 1995;Plesset & Chapman, 1971).
An unsteady cavity does not consist of a single bubble in a spherical shape, however. A group of bubbles that have different sizes and shapes form an unsteady cavity. So the jet velocity due to the collapse of a single bubble cannot be applied to agglomerated cavitating flow. There could be another approach. While a bubble is collapsing, the difference of pressure between the local and far-field pressure is replaced by the difference of pressure between the maximum and minimum pressure on the surface during a cavity shedding cycle (Chahine, 2014;Peters et al., 2015). The liquid density was replaced by the mixture density. Then the impact velocity can be defined as v impact = C ce (P/ρ m ) Here, the subscript m indicates the 'mixture phase' and C ce is a cavitation erosion coefficient. The cavitation erosion coefficient C ce is a function of reference and vapor pressure, density, dynamic viscosity, velocity, characteristic length and cavity shedding frequency.
To define the cavitation erosion coefficient as a nondimensional property, the Buckingham pi theorem was used. Finally, the coefficient is expressed as a function of the Reynolds, cavitation, and Strouhal numbers. Rather than in a superposition form, C ce is expressed in a multiplication form.
To consider the influence of turbulence, a fluctuating pressure, P , is used in calculating the difference of pressure in the impact velocity equation. The fluctuating pressure is expressed, following Hinze (1975), as Here, k is the turbulence kinetic energy. This expression relates the turbulence kinetic energy with cavitation erosion (Chitrakar, Dahlhaug, & Neopane, 2018). The impact velocity is obtained from the pressure variation, including the mean pressure from the solution of the RANS equations and the fluctuating pressure from the equations for the turbulence closure.

Critical velocity
If the impact velocity of the jet exceeds a critical velocity, cavitation erosion occurs on the solid surface. Therefore, how to define the critical velocity is a key factor in predicting cavitation erosion.
When a liquid mass strikes a solid surface, a normal shock wave propagates against the liquid stream, and behind the normal shock the pressure is increased and the liquid velocity is decreased. It is reasonable to assume that, for a moment, the liquid stream is an isentropic and incompressible flow (Lush, 1983), and the impact velocity may be derived excluding temperature. If the liquid mass is assumed to be infinitely wide, the phenomenon can be treated in one dimension as shown in Figure 1. In a steady flow, the continuity equation can be written as Figure 1 shows the configuration of the global coordinate system. It can also be represented using a relative coordinate system based on the shock wave. Here v indicates the initial velocity of liquid mass, u indicates the velocity behind the shock, c indicates the propagating speed of the shock wave, ρ l indicates the ambient liquid density, and ρ l2 indicates the density behind the shock wave. The momentum conservation equation can be written as Here, P amb indicates the ambient pressure. Both of the above equations can be solved by applying a pressure and density relation, which is where n = 7 and B = 300MPa (Lush, 1983). The ambient pressure can be neglected due to the large value of the liquid and shock speeds. The critical velocity of the liquid can be derived by applying the yield stress of the material, as it causes plastic deformation on a solid surface (Lush, 1983).
If the impact velocity due to the collapse of the cavity exceeds the critical velocity, erosion occurs. Note that the pressure and density variations, which generate the impact velocity, are selected as the key parameters. A much smaller velocity than critical velocity can also cause erosion of the solid surface because of fatigue. However, erosion by fatigue is not considered in this study.

Computational method
The present study focused on hydrodynamic machinery that are operated without a big change in temperature. Note that the temperature change caused by cavitation in hydrodynamic mechanical devices was known to be very small (Ha & Park, 2016;Shin, 2011). Incompressible cavitating flow solvers are widely used for cavitation erosion (Koukouvinis et al., 2016;Mouvanal et al., 2018;Peters et al., 2015). In an isothermal compressible cavitating flow code, the densities of the liquid and vapor phases were assumed to be a function of pressure only, while in a compressible cavitating flow code the densities were assumed to be function of temperature. In this study, an isothermal compressible flow code with cavitation was developed and applied.

Governing equations
The continuity equation, or the mass conservation equation, can be written as The Navier-Stokes equations, or the momentum conservation equations, can be written as where P indicates the static pressure and τ indicates the turbulent stress tensor. For turbulence closure, the ktwo-equation turbulence closure model was used. The same turbulence closure model has been employed for various cavitating flows (Li, Guo, & Li, 2018;Ma, Chen, Yu, & Jiang, 2018;Park & Rhee, 2012, 2013Yang, Teng, & Zhang, 2019).
The density and viscosity for mixture properties are calculated as a function of α, as where α v indicates the vapor volume fraction, and α l indicates the liquid volume fraction. When the cavitation is initiated, α v changes from 0 to 1. The density and viscosity of cavitating flows are computed by volume averaging in each cell. The transport equation of the volume fraction was solved to account for the dynamics of a cavity flow. The cavitation was governed by the kinetics and the thermodynamics of the phase change in the system. The transport equation for cavitating flows from Singhal, Athavale, Li, and Jiang (2002) can be expressed in a general form as follows Here R cond and R evap terms can be expressed as where the condensation coefficient (C c ) and the evaporation coefficient (C e ) are set as 0.01 and 0.01, respectively. v ch is the characteristic velocity and γ is the surface tension (Park & Rhee, 2013).

Pressure correction
In incompressible solver, the momentum conservation equation is discretized.
where the subscripts N and P indicate neighbor and owner cells, respectively. a N and a P are matrix coefficients of a cell. S indicates a source term. For simplicity of the off-diagonal part in the matrix of the momentum equation, a new operator, H( v), is used. The momentum conservation equation with the new operator is expressed as Substituting the expression v P into the incompressible continuity equation, ∇ · v = 0, yields This is the pressure equation form for incompressible flows. The pressure equation used is described in detail in Kissling et al. (2010).
In isothermal compressible flows, the density was not canceled in the conservation equation for mass. The density can be expressed as the pressure by the equation of state for the liquid and vapor phases, respectively. The vapor phase was expressed by the equation of state an ideal gas, and the liquid phase was expressed by the equation of state for an adiabatic or isentropic liquid (Park & Rhee, 2015) where c is the speed of sound.
The equation for conservation of momentum is the same for incompressible and isothermal compressible flows. Substituting the expression for v P into the isothermal compressible continuity equation, This is the pressure equation form for isothermal compressible flows. Here, using the harmonic relation, 1/P is computed: The developed solver for isothermal compressible flow and the solver of Senocak and Shyy (2002) are generally similar. The density variations of the two phases are different. Senocak and Shyy (2002) selected one equation to get the density for each phase: where primed variables represent the correction terms. C is an arbitrary constant. Here, C = 4 is adopted. In the liquid phase, variation of density is not allowed, while variation of density is allowed by the pressure variations in the vapor phase.

Numerical method
The first-order accurate implicit scheme was used for time derivative terms. The selected temporal scheme and carefully selected time step sizes have provided sufficient accuracy for engineering problems. The least-square gradient method was used for the gradient of solutions at cell centers. The total variation diminish (TVD) scheme with the van Leer limiter was applied to the convection terms. The general form of the TVD scheme was expressed as where the van Leer limiter was defined as where Here, P was the cell center of a focusing cell, A was the cell center of the neighbor cell, and f was the face center between P and A cells. The diffusion term was discretized by a central differencing scheme. A Pressure-Implicit with Splitting of Operators (PISO) algorithm was adapted for the velocity and pressure coupling and overall procedure. A Gauss-Seidel algorithm and multi-grid method were used to solve the discretized matrix.

Cavitation erosion coefficient
Cavitating flows in a converging-diverging nozzle and around a hydrofoil were simulated to obtain a cavitation erosion coefficient using the isothermal compressible cavitating flow solver.

Converging-diverging nozzle
In the Cartesian coordinate system adopted, the nozzle height, h, is given by where x is the horizontal channel location and H is the channel height of 0.05 m. The radius of the curvature at x/H = 0.3 was 2.5. A rectangular bump attached to the curvature at x/H = 0 was 1 mm long, 0.5 mm high, as shown in Figure 2. The height from the bump to the tunnel wall was 0.0195 m. More details on the converging-diverging nozzle used are described in Keil et al. (2011) and Kim, Chahine, Franc, and Karimi (2014). The cavity and cavitation erosion observation tests were done in the cavitation facility of the Technical University of Darmstadt (Keil et al., 2011). The test section of this tunnel was 0.5 m in length, 0.1 m in height, and 0.05 m in width. The Reynolds number (Re), based on the channel height and free stream velocity (U ∞ ), was from 2.75 × 10 5 to 3.5 × 10 5 , while the range of the cavitation number was 5.45-6.05. Table 1 lists the test conditions. Cavitation started to form behind the rectangular bump in the low-pressure region and convected downstream. The gap between the rectangular bump and the tunnel wall of the nozzle could be a parameter for cavitation erosion. In the present study, however, the effect of the gap was included in the effects of the Reynolds and cavitation numbers to simplify the cavitation erosion coefficient.
A two-dimensional solution domain was considered according to the test configuration. Figure 3 shows the solution domain extent (−10 ≤ x/H ≤ 20, 0 ≤ y/H ≤ 1). The fixed velocity condition was applied on the left boundary and the fixed reference pressure with the    extrapolated volume fraction and velocity condition was applied on the right boundary. a no-slip condition on the nozzle walls was applied, and a slip-condition on the tunnel walls was applied. Figure 4 shows a structured grid of 14,000 cells with y+ of 30.
To tame the numerical instability from the phase change, the computations of the cavitating flow started from the solutions for the single phase at the same Reynolds number and a large cavitation number. The non-dimensional time step, tU ∞ /H, was 0.00011, which corresponded to a time step of 10 −8 and the simulation time was 6 s. The time advancing was done when the normalized residuals for the solutions had dropped by six orders of magnitude. Figure 5 shows the non-dimensionalized sheet cavity length, l C /H, for various Reynolds and cavitation numbers, and compared with the experimental data (Keil et al., 2011). In the computations with a cavitation number of 6, and Reynolds numbers of 2.35 × 10 5 and 3.5 × 10 5 , the computed length of the cavity was slightly longer than that of the experiment. In the overall computations, the results showed quite a close agreement. Figure 6 shows the pressure difference, which includes the mean and fluctuating pressure differences, on the nozzle surface behind the rectangular bump. The pressure difference is defined as the difference between the maximum and minimum pressure values on the nozzle surface when the cavity is generated and collapsed. A large pressure difference was observed in a strongly fluctuated cavitating flow. The cavity dynamics were seen in a wider range than the sheet cavity length. Figure 7 shows the pressure difference in the frequency domain. The fundamental frequency for four cases was 1.66 Hz. Figure 8 shows the volume fraction contours with cavity shedding at a cavitation number of 5.45 and Reynolds number of 2.75 × 10 5 . The cavity was generated behind the rectangular bump. As the cavity convected downstream, the cavity behind the bump was detached from the agglomerated cavity and the pressure increased significantly at this point. Keil et al. (2011) measured the damage extent, as shown in Figure 9. The pressure difference is related to the damage extent. Using the damage  Figure 9, the pressure difference in the range of x/H can be read from Figure 6, from which the impact velocity can be calculated. By comparing the impact and critical velocities, the cavitation erosion coefficients can be calculated for various Reynolds numbers and cavitation numbers. Table 2 summarizes the cavitation erosion coefficients. The cavitation number was the dominant factor influencing the extent of damage from cavitation erosion. In the experiments (Keil et al., 2011), the observed cavitation erosion extent depended on the cavitation number, not the Reynolds number. Figure 10 shows the hydrofoil with a circular-shaped leading edge and a wedge-shaped trailing edge. The hydrofoil was symmetrical in the span direction and had 0.0016 mm constant thickness. The 0.1079 m chord and the 0.05 m span lengths were determined by the size of the cavitation tunnel. In the Cartesian coordinate system selected here, the chordwise direction was aligned to the positive x-axis, the spanwise direction was aligned to the positive y-axis, and the normal direction of the surface was aligned to the positive z-axis. The observation tests of the cavity and cavitation erosion were done in the cavitation facility of the Technical University of Darmstadt (Dular & Coutier-Delgosha, 2009). The test section was 0.5 m in length, 0.1 m in height, 0.05 m in width. A developed cavitating flow around the hydrofoil with 5°incidence angle was observed. The incidence angle could be a parameter for cavitation erosion. In the present study, however, the effect of the incidence angle was included in the effects of the Reynolds and cavitation numbers to simplify the cavitation erosion coefficient.

Hydrofoil
The cavitation number (σ ) of 2.0 was used, while the Reynolds number (Re), based on the chord length of 0.1079 m and free stream velocity (U ∞ ) of 13.0 m/s, was 1.04 × 10 6 The two-dimensional solution domain (−10 ≤ x/C ≤ 10 and −0.5 ≤ y/C ≤ 0.5) was considered to conforme the test setup (see Figure 11). The fixed velocity condition was applied on the left boundary and the fixed reference pressure with the extrapolated volume fraction and velocity condition was applied on the right boundary. A no-slip condition on the hydrofoil was applied and a slip-condition on the tunnel walls was applied. Figure 12 shows a structured grid of 12,000 cells with y+ of 30. Two hundred and fifty cells were used on the surface and 60 cells were used in the normal direction.
To tame the numerical instability from the phase change, the computations of the cavitating flow started from the solutions of the single phase at the same Reynolds number and a large cavitation number. The cavity was generated around the leading edge and convected downstream. Figure 13 shows the pressure difference ( P), which includes the mean and fluctuating pressure differences, on the hydrofoil surface. Figure 14 shows the volume fraction contours with the cavity. The shedding cavity increased the pressure significantly. The cavitation erosion coefficient could be calculated by the same procedure for the converging-diverging nozzle. The damage extent with Reynolds number of 1.38 × 10 6 and cavitation number of 2 was x/C = 0.78 (Dular & Coutier-Delgosha, 2009), and the cavitation erosion coefficient was calculated to be 28.72.

Determination of cavitation erosion coefficient
The cavitation erosion coefficient, C ce , was expressed as a function of the Reynolds, cavitation, and Strouhal numbers. The authors considered that the Stouhal number should be excluded because high pressure variation finally caused cavitation erosion. Other computational studies on cavitation erosion have not yet considered the  Strouhal number (Koukouvinis et al., 2016;Mouvanal et al., 2018;Peters et al., 2015). In this paper, the Strouhal number was treated as constant.
In the function of the Reynolds number, the characteristic length could be selected as the cavity or geometry  length. The cavity length was more reasonable in smallscale cavity flow. However, in this study, the geometrybased characteristic length was selected to develop a more practical method.
To develop the practical method for the prediction of cavitation erosion, a metamodeling method was considered. Metamodels have been widely used with optimization techniques in fluids and structural engineering applications because they evince excellent performance prediction. Response surface methodology (RSM) and   kriging metamodels were widely used. A RSM metamodel was based on an optimization method that built linear or quadratic local approximations, whereas a kriging metamodel was more flexible than polynomial models and less sensitive than radial basis functions (Meckesheimer, Booker, Barton, & Simpson, 2002). Also, a kriging metamodel showed accuracy in highly nonlinear problems. Thus a kriging metamodel was selected for the present study. Figure 15 show the cavitation erosion coefficient predicted by the kriging metamodel. The predicted cavitation erosion coefficient was available for the cavitation number range of 2-6.05 and Reynolds number range of 2.75 × 10 5 to 13.8 × 10 5 . In Figure 15, the predicted cavitation erosion coefficient included some error beyond the design points due to the small number of design points, that is, the selected cavitation and Reynolds numbers for the converging-diverging nozzle and hydrofoil cases. Also, beyond the cavitation and Reynolds numbers range, the predicted cavitation erosion coefficients included some errors due to extrapolation. The cavitation erosion coefficients for various cavitation and Reynolds numbers are listed in Table 3.

Application
The developed and validated practical method for the cavitation erosion prediction was applied to a marine   propeller. As the object propeller of this study, KP505 was selected; that had been designed for a container carrier (Seo, Lee, Han, & Rhee, 2016). A five-bladed propeller with NACA66 foil section and diameter of 0.25 m was selected. The advance ratio (J), based on a propeller diameter of 0.025 m, revolution rate of 25 rps, and free stream velocity (U ∞ ) of 3.125 m/s, was 0.3. The Reynolds number (Re) was 7.22 × 10 5 , while the cavitation number (σ ) was 1.5. A time step size of 10 −8 was used.
In this study, a moving reference frame was selected to express the revolution of the propeller. Only one blade with a periodic boundary condition was modeled because of the open water condition. The computational domain extent was 12D in length and 3D in radius, where D is the diameter of the propeller. The inlet and outlet boundaries from the propeller were 4D and 8D, respectively.
The fixed velocity condition was applied on the left boundary and a fixed reference pressure with the extrapolated volume fraction and velocity condition was applied on the right boundary. a no-slip condition on the blade surface was applied, and a periodic condition on the side boundary was applied. Figure 16 shows a hybrid grid of 0.4 million cells. One hundred cells in the chordwise and 400 cells in the spanwise directions were used on the blade surface as shown in Figure 17. Tetrahedral cells were used in the subdomain including the blade, while hexahedral cells were used in the other subdomain with a simple domain shape. By means of this hybrid zonal approach, the mesh generation difficulties of complex geometries were excluded. Dense cells near the blade tip and trailing edge were used due to the hybrid mesh approach.
To tame the numerical instability arising from the phase change, the computations of the cavitating flow started from the solutions of the single phase at the same Reynolds number and a large cavitation number. To validate the selected mesh, boundary conditions, and numerical methods, the thrust coefficient (K T = T/ρn 2 D 4 ) and the torque coefficient (K Q = Q/ρn 3 D 5 ) for various advance ratios (J = V/nD) of the open water tests were compared with experimental data (Seo et al., 2016) in Figure 18. Here, T is the thrust produced by the propeller blades, Q is the torque applied to the propeller blades, n is the revolution per second, and D is the propeller diameter. For thrust, there was good agreement with the data, while for torque there was 5% difference from the data. Figure 19 shows the pressure coefficient contours on both sides of the blade in the non-cavitating flow. High and low pressure coefficient peaks were observed at r/R = 0.7 on the pressure and suction blades, respectively. Figure 20 shows the generated cavity on the suction side and the predicted cavitation erosion extent. Using a volume fraction of 0.5, the cavity interface was captured. The cavity was generated at the leading edge of r/R = 0.7 and propagated to the blade tip. The cavitation erosion extent was predicted near the cavity closure at r/R = 0.8. Figure 21 shows a snapshot of the cavitation erosion of a propeller blade (Pfitsch, Gowing, Fry, Donnelly, & Jessup, 2009). The predicted cavitation erosion extent was similar to the measured one.

Conclusions
In this study, a practical method for the prediction of cavitation erosion was developed. When the impact velocity was larger than the critical velocity, it was predicted that cavitation erosion could be observed. To define cavitation erosion, impact and critical velocities were introduced.
To simulate cavitating flows, an in-house computational fluid dynamics (CFD) code, called SNUFOAM-Cavitation, was developed using the OpenFOAM toolkit. Cavitating flow solvers for both incompressible and isothermal flows were developed. To validate the developed solvers, sheet and cloud cavitating flows were simulated and validated against experimental results. For the cavitation in the flow with cavitation erosion, the isothermal compressible flow solver was selected and validated for the cavitating flows. To obtain the cavitation erosion coefficient, the cavitating flow in the nozzle and around the hydrofoil was simulated. The cavitation erosion coefficient was derived by metamodeling and curve-fitting methods.
Finally, the method developed was validated by applying it to cavitation erosion on a propeller blade. The predicted damage extent was nearly the same as that on a damaged full-scale propeller blade. It was confirmed that the proposed method can predict cavitation erosion observed on the hydrodynamic blades of turbines, pumps, and marine propellers.

Disclosure statement
No potential conflict of interest was reported by the authors.