New method for predicting the incipient cavitation index by means of single-phase computational fluid dynamics model

Due to its serious consequences, cavitation in pipeline systems is one of the main concerns of engineers. In this context, the incipient cavitation index defined in many international standards such as those from International Society of Automation and International Electrotechnical Commission is a very important parameter as it accounts for the onset of cavitation. However, the standards define only the experimental determination of the incipient cavitation index, which is furthermore difficult to perform due to the considerable technical and economical burden of the experimental tests, above all in case of large-size systems. In this work, we propose a new method for predicting the incipient cavitation index by means of computational fluid dynamics, avoiding the need of making any experiment. The method is based on the generalized pressure criterion and requires only one Reynolds-averaged Navier–Stokes simulation of single-phase incompressible flow to provide estimates of the incipient cavitation index. The method is applied to predict the incipient cavitation index of multi-hole orifices with different geometrical characteristics, in terms of equivalent diameter ratio (0.40–0.70), relative thickness (0.73, 1.00), and number (13–52) and disposition of the holes. The experimental data revealed the reliability of the method. Its applicability was also confirmed for more complex geometries, and an application example regarding a control valve is briefly illustrated at the end of this article.


Introduction
Cavitation in pipeline systems is one of the main concerns of engineers as it can result in serious damages of the plants. It is well known that cavitation is more likely to occur in case of a high variation in the flow velocity, as in correspondence of either moving parts (vanes of turbines, pump impellers) or non-moving parts (pipeline fittings, valves). Among the non-moving elements, multi-hole orifices ( Figure 1) are widely used to reduce flow nonuniformities or, placed side by side a valve, to attenuate noise and vibrations and reduce the cavitation itself.
Cavitation is a two-stage process associated with the flow of liquids, which is based on the formation and collapse of vapor cavities in low-pressure regions, causing noise and vibrations of the plant. According to the international standards, different cavitation regimes are identified. Following ISA-RP75. 1 they are: (1) Regime I: absence of cavitation; (2) Regime II: incipient cavitation, the onset of cavitation, when only small vapor bubbles are formed in the flow stream; this condition is detected by measuring increased sound pressure level (SPL) or vibration; (3) Regime III: constant cavitation, involving a sufficiently large volume of vapor to produce a uniform and constant level of cavitation; and (4) Regime IV: maximum vibration, that is, the level of cavitation associated with peak vibration measurements. The cavitation level that causes damages to the components of a plant is difficult to determine a priori, and it is usually indicated by the manufacturers on the basis of experience. Clearly, the more conservative choice consists in restricting all operations to a cavitation-free regime even if typically the incipient cavitation is of mild level and therefore may be acceptable in the design phase. 1,2 Cavitation is often studied by means of specific indexes, such as the different cavitation numbers proposed in the literature, [2][3][4] or the net positive suction pressure common in the context of pumps and turbines. 5 In this work, we make reference to the cavitation index s reported in ISA-RP75.23-1995 1 and Tullis 2 : where P U and P D are the pressure upstream and downstream the device, respectively, and P V is the vapor pressure. The sections U and D are placed 2 pipe diameters and 6 pipe diameters from the device, respectively (see Figure 1, in which a multi-hole orifice is the device being tested). Different threshold values of s are defined: the incipient cavitation index s i , at the beginning of Regime II; the constant cavitation index s C , between Regimes II and III; and the maximum vibration cavitation index s mv between Regimes III and IV. The incipient cavitation index-subject of this article-is particularly important as it identifies the real detectable onset of cavitation and it is a conservative design limit. According to ISA-RP75.23-1995 1 -and also IEC 60534-8-2 6 -the incipient cavitation index s i can be estimated experimentally by drawing in a semi-log plot the SPL versus the cavitation index s and identifying the value of s which corresponds to a sudden increase in the slope of the data ( Figure 2). Actually, the same procedure can be applied considering the pipe wall acceleration instead of the SPL.
A plot of the flow rate Q versus the square root of the pressure drop across the device ffiffiffiffiffiffi ffi DP p ( Figure 3) allows identifying another level of cavitation, the choking cavitation, for which the flow rate is limited to some maximum value (horizontal line in Figure 3), depending on both the geometry of the device and the pressure upstream of it. The onset of choking cavitation is associated with a specific cavitation index (s ch ) which cannot be determined from Figure 2. Conversely, it is hard to identify in Figure 3 the cavitation regimes previously defined. The no-cavitation regime is evidenced by the straight line through the origin; the inception of cavitation-that is, the incipient cavitation index s i -corresponds to the starting point of the curve linking the no-cavitation regime with the choked flow regime, in most cases not distinguishable on the basis of the experimental results.
Many researchers investigated experimentally the cavitation characteristics of different devices; among them, Tullis, 2,7 Kim et al., 8 Testud et al., 9 Maynes et al., 10 and Malavasi et al. 11 focused on single-and multi-hole orifices. However, the already mentioned procedure reported in the ISA-RP75.23-1995 1 and IEC 60534-8-2 6 for estimating the incipient cavitation index by means of SPL or pipe wall acceleration measurements is not free from difficulties mainly arising from the fact that the change in slope in the plot of Figure 2 is not always clearly identifiable, even in the presence of a large number of sampling points. In addition, the economical burden of the tests may be considerable especially for large-size systems. Computational fluid dynamics (CFD) has been explored as a possible alternative approach to deal with cavitation. Some authors 12-15 made use of two-phase models to describe the process of rising and implosion of the vapor bubbles under cavitating flow conditions. These models appear able to reproduce the position of the cavitation pockets as well as the flow field of the cavitating fluid very accurately, but compared to a single-phase Reynolds-averaged Navier-Stokes (RANS) model they introduce numerical difficulties and are more time-consuming. Above all, the estimation of the cavitation limits (which are the really significant parameters in most engineering applications) is not straightforward when using these two-phase models; focusing on the incipient cavitation index, several flow conditions may need to be simulated in order to determine in the end the one which corresponds to the onset of cavitation.
In this work, a new method for estimating the incipient cavitation index using CFD is proposed. The method relies on the generalized pressure criterion that we previously developed to predict the presence of cavitating conditions in single-hole orifices. 16 The proposed procedure requires just one RANS simulation of singlephase incompressible flow to estimate the incipient cavitation index. This feature makes the new method attractive and very useful for engineers, especially with regard to its application to complex geometries and large-size systems.
This article essentially concerns the application of the method to multi-hole orifices which, despite their apparent simplicity, are devices of considerable interest for the applications. The good predictive capacity of the proposed method is tested with respect to our own experimental data, paying particular attention in quantifying the uncertainties of both measurements and computations. Table 1 reports the geometrical characteristics of the orifices considered, which differ in terms of diameter d h , number n h and distribution of the circular holes, and thickness of the plate t ( Figure 1). The internal diameter of the pipe D is equal to 77.9 mm. The values resulted in different combinations of the following characteristic dimensionless parameters: equivalent diameter ratio b = ffiffiffiffi ffi n h p Á d h =D and relative thickness t=d h .
The remainder of the article is divided into four sections, followed by the conclusion. Section ''The proposed method'' explains the proposed method for estimating the incipient cavitation index. Section ''Experimental setup'' describes the experimental design and methods used to collect the data, together with the quantification of measurements uncertainty. Section ''Numerical model and consistency of the numerical solution'' focuses on the numerical model with reference to the key settings (domain configuration, mesh resolution, differencing schemes, boundary conditions, turbulence modeling). Section ''Results'' proves the capacity of the proposed method to correctly estimate the incipient cavitation index in multi-hole orifices by comparison to our own experimental data. The end of the same section is dedicated to showing the reliability of the method in case of more complex devices and briefly reports an application example regarding a control valve.

The proposed method
As already mentioned, the proposed method relies on the generalized pressure criterion that we developed in previous works. At first, we will briefly illustrate the fundamentals of such criterion (the reader may refer to Malavasi and Messa 16 for more details). Afterward, we will explain the procedure for estimating the incipient cavitation index, which constitutes the novelty of this article.
The starting point of the generalized pressure criterion is that when running a RANS simulation of  incompressible flow, the computed pressure field P Ã is defined up to a constant; in other words, it is possible to add or subtract the same amount of pressure to every cell without affecting the other flow variables. The computed pressure field is uniquely determined by imposing a fixed pressure in one cell, typically at the outlet section; by doing so, the computed pressure field P Ã is physically admissible only if P Ã M .P V , P Ã M being the minimum computed pressure in the domain and P V the vapor pressure. This characteristic implies that when simulating the incompressible flow through a device, the relationship between the flow rate Q Ã and the square root of pressure drop ffiffiffiffiffiffiffiffi ffi DP Ã p is a straight line through the origin that extends up to ffiffiffiffiffiffiffiffi ffi DP Ã p ! + '; therefore, the deviation from such straight line depicted in Figure 3, caused by cavitation, cannot be reproduced by a single-phase incompressible RANS model. The basic idea of the generalized pressure criterion is to consider the difference between P Ã M and P V as an indicator of the presence or absence of cavitation: if P Ã M .P V , the system is not subjected to cavitation; if P Ã M \P V , the system is subjected to cavitation. In the latter case, the modulus of the difference P Ã M À P V is assumed as an indicator of the level of cavitation. In Malavasi and Messa, 16 we noticed that with reference to some works of the research group of Joseph and colleagues, 17-19 the equality P Ã M = P V is not suitable to identify the physical onset of cavitation, that is, the formation of the very first vapor bubble; therefore, we limited the applicability of the method to cases in which the difference P Ã M À P V is ''sufficiently'' far from 0. However, according to many technical standards (e.g. ISA-RP75.23-1995 1 and IEC 60534-8-2 6 ), the conditions under which cavitation begins-expressed by the incipient cavitation number s i -are deduced by the first detectable effect of cavitation such as vibrations and noise rather than from detection of the first vapor bubble. Such consideration, together with the technical outline taken by this article, led us to the possibility of evaluating the incipient cavitation index from the condition P Ã M = P V . The numerical estimation of the incipient cavitation index requires adding a constant value to the computed pressure field P Ã obtained from the simulations; as per the considerations already reported, such transformation is allowed when solving the single-phase incompressible RANS equations. Hereafter, we will denoteP as the pressure field obtained from P Ã after adding a certain constant pressure in each cell, which will be from time to time specified. For the estimation of s i , the proper transformation is so that the minimum value ofP is equal to the vapor pressure P V , and the numerically estimated incipient cavitation index s Ã i is It is worth noticing that although it is defined making reference to the translated fieldP, s Ã i can be obtained directly from P Ã , that is, from the output of the simulations, avoiding the need of actually evaluat-ingP. In fact, substitution of equation (2) in equation (3) yields In section ''Results,'' we will show that the values of s Ã i estimated from equation (2) do not depend significantly on the flow rate; therefore, one RANS simulation of single-phase incompressible flow is sufficient for predicting the incipient cavitation index.

Experimental setup
The experiments have been performed in a pilot plant located at Pibiviesse S.r.l., Nerviano, Italy. The rig, sketched in Figure 4, consists of 10$ and 12$ steel pipes, supplied by a pump that guarantees pressures up to 10 bar at the reference section upstream the plate. Control valves placed upstream and downstream the test area allowed setting the proper fluid-dynamic conditions for each experimental test. Pressure drop has been measured with a series of absolute and differential pressure transducers in reference sections located, as reported above, 2D upstream and 6D downstream the plate according to the ISA-RP75. 23-1995. 1 Several other measurement points were placed up to 8D downstream the device to verify the proper occurrence of pressure recovery. The accuracy in the pressure measurements was 44, 155, or 260 Pa depending on the transducer used. The flow rate has been measured by an electromagnetic flowmeter of 10$, placed upstream the test section. As indicated by the manufacturer of the flowmeter on the basis of laboratory tests, the uncertainty on the pipe bulk-mean velocity V p was taken as 0:002V p . During the tests, the water temperature has been measured in order to monitor the values of density, viscosity, and vapor pressure. Either the SPL or the pipe wall acceleration was considered to characterize the cavitation regime. In the former case, we employed a sound meter level with sensitivity of 50 mV/Pa. In the latter one, we employed a single PCB Piezotronics 352A60 accelerometer with sensitivity of 1.02 mV/(m/s 2 ), range equal to 6 4906 m/s 2 and broadband resolution equal to 0.02 m/s 2 .
The tests have been performed maintaining constant pressure at the upstream reference section, and decreasing the downstream pressure, in order to increase the Reynolds number and display the cavitation phenomena. The vapor pressure was taken as 20.98 bar, the temperature being close to 20°. The determination of the uncertainties of the experimental data was performed with respect to the International Organization of Standardization-GUM; 20 the maximum relative error on the cavitation index (equation (1)) was found to be about 0.8%. A more detailed description of the test rig and the experimental procedure is reported elsewhere. 21

Numerical model and consistency of the numerical solution
The PHOENICS 2009 CFD code was used to simulate the flow. PHOENICS uses the finite-volume method. The calculations were performed following the ellipticstaggered formulation in which the scalar variables are evaluated at the cell centers and the velocity components at the cell faces. Central differencing is employed for the diffusion terms, while the convection terms are discretized using the hybrid differencing scheme of Spalding, 22 as in previous studies concerning similar flows. 23 A parallel version of the SIMPLEST algorithm of Spalding 24 was used to solve the finite-volume equations. The code was run until the sum of the residuals over the whole solution domain was less than 0.1% of the variables, based on the total inflow. Convergence was reached within 9000 iterations for all the multi-hole orifice simulations. Additional checks revealed that the minimum computed pressure P Ã M fluctuated less than 0.01% between successive iteration cycles. Figure 5 shows the configuration of the computational domain. Due to the geometrical symmetry of the multi-hole orifices simulated, and because the steady RANS was being solved, only a quarter of the pipe section was simulated. At the inlet, a 1/7th power law velocity profile was set, while the distribution of the turbulent kinetic energy and of its dissipation rate was derived imposing a turbulent intensity of 5% and an inlet mixing length equal to 7% of the hydraulic diameter. At the outlet, the normal gradients of all variables and the value of the relative pressure were set to 0. As already remarked, with a single-phase incompressible RANS solver, the static pressure at the outlet section must be specified only to obtain a unique pressure   field, which, otherwise, would be defined up to a constant. This property allowed us running all the simulations with a zero relative pressure at the outlet, and then, if necessary, add or subtract an arbitrary pressure to the whole computed pressure field. The inlet boundary was located 7D upstream the plate; such distance was verified to be sufficient to achieve fully developed flow conditions, confirming the results of Morrison et al. 25 20D were modeled downstream the device to provide accurate estimation of the region in which the pressure field is affected by the plate. Specific tests were performed to guarantee that such distance is sufficient for all the case study considered.
No-slip conditions were imposed on the walls of both the pipe and the plate. In the near-wall region, the non-equilibrium wall function of Launder and Spalding 26 was set, as per the indications of previous workers who dealt with similar kind of flows. 23,27 In all cases, it was checked that the application of such condition was consistent with the non-dimensional distance of the first grid points from the walls.
To close the system of equations, use is made of the k-e renormalization group (RNG) turbulence model which, taking care of the effect of rapid strain in complex flows, allows predicting the gross flow behavior in the recirculation region downstream the plate. The choice was made on the basis of previous literature studies. Koronaki et al. 28 report that, for separated flows, among the two-equation turbulence models, the k-e RNG is the one which performs best when coupled with non-equilibrium wall functions.
A block-structured mesh of about 4.1 million cells was used to discretize the domain. In the section normal to the flow, the grid consisted of 87 3 87 uniformly distributed cells. Along the flow direction, 255 and 258 cells were set upstream and downstream the plate, respectively, while the number of cells in the inner region was varied between 28 and 35 according to the thickness of the plate. A grid sensitivity study was made to choose the optimum discretization for the investigations. The tests were performed for the case of Plate C in Table 1 and Q = 15.1 L/s. The flow simulation was repeated employing five different grids in which the number of cells was increased from 0.5 to 9 million, and the effect of the number of cells upon the minimum compute pressure P Ã M was investigated. The results, reported in Figure 6, indicate that the third mesh used, consisting of about 4 million cells, allows obtaining a grid-independent solution in terms of the target parameter P Ã M , in the sense that further increase in the number of cells does not result in appreciable variation in P Ã M itself. Moreover, the solution is strongly dependent upon the disposition of the cells. In particular, the value of P Ã M appears considerably affected by the axial dimension of the cells within the holes, which had to be set below 0.0027D to procure a grid-independent solution.
As a preliminary step, the good predictive capacity of the numerical model was checked with reference to the Euler number Eu = (P U À P D )=2rV 2 ; good agreement between computations and measurements was observed for all the devices in Table 1. However, the results of this validation are not reported in this article, which is focused on the cavitation index.

Results
Before using the procedure described in section ''The proposed method'' to estimate the incipient cavitation index of the four orifices whose geometrical characteristics are reported in Table 1, we applied the generalized pressure criterion to check its capability to predict whether, under certain flow conditions, multi-hole orifices are subjected to cavitation or not, thereby extending the validation study which, in our earlier article, 16 was limited to single-hole orifices. However, we no longer made reference to the plot of the flow rate Q versus the square root of pressure drop ffiffiffiffiffiffi ffi DP p (Figure 3) but, rather, we considered the plot of SPL versus cavitation index (Figure 2). The analysis is performed for the multi-hole orifices referred to as A and B in Table 1. Details of the flow conditions simulated are reported in Table 2. Since the cavitation index s (equation (1)) depends on the plant pressure P U , for each device we considered a series of experimental data characterized by the same value of P U , namely, 2.85 and 2.17 bar for Plates A and B, respectively. As briefly remarked in section ''The proposed method,'' the generalized pressure criterion in its original formulation 16 estimates the presence or absence of cavitation from the difference between the minimum pressure computed in the Figure 6. Grid independence study for the case of Plate C in Table 1 and Q equal to 15.1 L/s. Effect of mesh resolution on the minimum computed pressure P Ã M . The vertical line represents the discretization considered in this work.
domain and the vapor pressure. However, in order to validate this procedure with respect to experimental data, we cannot refer directly to P Ã -the output of the simulations, in which a zero static pressure is arbitrarily imposed at the outlet section-but we must instead make reference to a transformed pressure fieldP evaluated from P Ã in such a way that the upstream (plant) pressureP U equals that of the experiments P U . Therefore, for the comparison between the numerical and experimental cavitation indexes, the proper transformation isP It has already been remarked that the addition of the constant value P U À P Ã U to the whole computed pressure field is allowed when solving the incompressible RANS. Therefore, the indicator of cavitation isP M À P V and the corresponding cavitation index is calculated ass =P U À P V =P U À P D . After substitution of equation (5), the above-reported quantities becomẽ The differenceP M À P V is reported in Table 2 for all the flow conditions considered. According to the generalized pressure criterion, a negative sign ofP M À P V means that the system is subjected to cavitation; therefore, the generalized pressure criterion predicts that the flow conditions corresponding to points A 3 , A 4 , B 2 , B 3 , and B 4 are subjected to cavitation while those corresponding to points A 1 , A 2 , and B 1 are not. Table 2. Application of the generalized pressure criterion in its original formulation 16 to Plates A and B in Table 1 and different flow rates.

Device
Point Q(L=s)P U = P U (bar)P D (bar)P M (bar)P M À P V (bar) Cavitation?s( À ) The vapor pressure is P V = À 0:98 bar, as estimated during the experiments. Here, the transformation is as of equation (5), and therefore, the differenceP M À P V and the cavitation numbers are given by equations (6) and (7), respectively. Figure 7. Measured values of SPL as a function of the cavitation index s for plates A and B reported in Table 1. The vertical lines represent the estimates of the cavitation indexs deduced from the numerical simulation of the flow conditions reported in Table 2. P U is the pressure at the upstream reference section during the experimental tests. Figure 7 shows the trend of SPL as a function of the cavitation index for the experimental data considered for comparison. Vertical lines representing the numerically derived valuess for the flow conditions of Table 2 are also displayed. Following the ISA-RP75.23-1995 1 standard, since the presence of cavitation is accompanied by a sudden increase in SPL, the transition between Regimes I and II in Figure 2 corresponds to a variation in the slope of the measured data of Figure 7. It appears thus rather clear that, as correctly predicted by the generalized pressure criterion, the flow conditions A 3 , A 4 , B 2 , B 3 , and B 4 are subjected to cavitations, while A 1 and A 2 are not. The identification of the cavitation regime of the flow condition B 1 is rather less evident; nevertheless, the prediction of the criterion is compatible with the experimental evidence.
The new method for predicting the incipient cavitation index is now checked by comparison to the experimental data regarding all the four plates in Table 1. Since the incipient cavitation indexs i -the abscissa of the point which corresponds to a change of the slope in Figure 2-does not depend significantly on the plant pressure P U , 11 several series of experimental data characterized by different values of P U were considered together to allow a more accurate estimation ofs i . As far as the CFD simulations are concerned, three different flow rates were simulated. For each run, the incipient cavitation index was calculated directly from the computed pressure field P Ã by means of equation (4); the values obtained are reported in Table 3.
For all plates, the estimated incipient cavitation index s Ã i seems to be only marginally affected by the flow rate Q, with deviations always below 8%. This feature provides to the method proposed in section ''The proposed method'' a simplicity which makes it particularly attractive for the applications. As already mentioned in section ''Introduction,'' equation (4) allows obtaining an estimation of the incipient cavitation index, for a given geometry of the device, just running one single-phase RANS simulation. No particular conditions on the flow rate and reference pressure are needed. Figure 8 reports the measured SPLs versus the cavitation index together with the predictions of incipient cavitation index obtained by averaging the values reported in Table 3, corresponding to different flow rates. The predictions of equation (4) show agreement with our experimental evidence. We emphasize that whilst a large number of experimental points are required (and may not be even sufficient) to estimate s i with proper accuracy by ISA-RP75.23-1995 1 procedures, only one numerical simulation is sufficient with our method.
The method proposed in this article has been applied to several control valves of different types and geometrical characteristics and proved capable in providing reliable estimates of the incipient cavitation index also for these devices. As an example, we report here the results obtained for the control valve sketched in Figure  9, in which two perforated plates are inserted in the valve body. Particularly, Figure 10 shows the measured value of pipe wall acceleration (normalized with respect to its maximum value) versus the cavitation index s at a valve opening angle of 40°. The data were collected in the test rig described in section ''Experimental setup.'' As in Figure 8, a vertical bar is used to indicate the estimates of the incipient cavitation index obtained by the proposed methodology, that is, by applying equation The incipient cavitation index s Ã i is evaluated from equation (4).
(4). The good agreement between estimates and experiments confirms the reliability of the proposed method also for cavitation in more complex geometries.

Conclusion
In this work, a new method is proposed for predicting the incipient cavitation index-introduced in many technical standards 1,6 as an indicator of the onset of cavitation-by means of CFD simulations. The method is based on the generalized pressure criterion that we developed and reported in an earlier article 16 in order to estimate whether, under certain flow conditions, a  Table 1, together with the estimate of the incipient cavitation index s Ã i obtained from equation (4). The values of s Ã i are the average of those reported in Table 3. P U is the pressure at the upstream reference section during the experimental tests. Some of the experimental data of Plates A and B have been previously reported in Figure 7.  single-hole orifice is subjected to cavitation or not. The strength of the new method relies in its simplicity, as it allows estimating the incipient cavitation index from just one single-phase incompressible RANS simulation without performing any experimental test. This feature makes the proposed method very attractive for engineers, especially in view of its application to complex geometries that are very hard to investigate experimentally, above all in case of large-size systems.
The case of multi-hole orifices is considered in this article, referring to our own experimental data for comparison. After verifying that the sign of the difference between the minimum computed pressure and the vapor pressure is a suitable parameter to provide an estimation of the presence of cavitating conditions also in multi-hole orifices, extending to a more complex flow the results of our previous work, 16 the proposed method is employed for estimating the incipient cavitation index (equation (4)). Good agreement between predicted and measured values of the incipient cavitation index was obtained for multi-hole orifices with different geometrical characteristics as well as for more complex devices including different types of control valves, revealing the reliability of the method and its relevance for engineering purposes.