**New Advances of Cavitation Instabilities**

Editor **Florent Ravelet**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Florent Ravelet LIFSE Arts et Metiers Institute of Technology Paris France

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Applied Sciences* (ISSN 2076-3417) (available at: www.mdpi.com/journal/applsci/special issues/ cavitation instabilities).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-1546-5 (Hbk) ISBN 978-3-0365-1545-8 (PDF)**

© 2021 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## *Editorial* **Editorial for Special Issue: New Advances of Cavitation Instabilities**

**Florent Ravelet**

Arts et Métiers Institute of Technology, CNAM, LIFSE, HESAM University, 75013 Paris, France; florent.ravelet@ensam.eu

**Abstract:** This editorial presents the main articles published in the Special Issue: New Advances of Cavitation Instabilities.

**Keywords:** water jet cavitation; cavitation peening; hydrodynamic cavitation; compressible two-phase flow; Kelvin-Helmholtz instability; slamming; vortex rope; bulb turbine; hydrofoil; Francis turbine

Cavitation refers to the formation of vapor cavities in a liquid when the local pressure becomes lower that the saturation pressure [1]. In many hydraulic applications, cavitation is considered as a non-desirable phenomenon, as far as it may cause performance degradation, vibration problems, enhance broad-band noise-emission and eventually trigger erosion.

In this Special Issue, the deliberate use of cavitating jets for peening treatment of materials has been considered by Yoshimura et al. [2] and by Soyama [3]. This very interesting review shows how the aggressive intensity of a cavitating jet varies non-monotonically with the cavitation number in accordance with the variation of mixture sound speed for different regimes and that there exists an optimum.

These findings illustrate the importance of instabilities in cavitating flows. The cavitation is usually associated with flows with very high Reynolds numbers. This formally unsteady flow type is very sensitive to disturbances. In the ordinary case of a 2D profile, it is rightly recognized that for sufficiently low cavitation numbers, periodic or quasi-periodic cavity shedding arises. On more complex problems, like cavitation in inducers, several other periodical behaviors can emerge. Firstly, the local intrinsic flow instabilities will depend on the region that presents hydrodynamic cavitation: tip-leakage vortices, backflow vortices, or blades suction surface. Secondly, interaction between adjacent blades can lead to rotating instability, with cells of various sizes that propagate from blade to blade. Finally, system instabilities are also observed, because of the blockage linked to the volume variation of the pockets, or to a possible positive slope of the inducer characteristics linked to a change in the angle of attack on the blades. A more proper understanding of these instabilities is of crucial interest, especially in the field of turbomachinery, still motivating applied and fundamental research on cavitation instabilities.

In this Special Issue, Pipp et al. [4] and Ravelet et al. [5] present new findings in cloud shedding mechanisms arising in a 2D venturi geometry, respectively at very small scale (micro channel) and at at very low Reynolds numbers, where the development of the Kelvin-Helmholtz instability is evidenced. Viitanen et al. [6] studied numerically the cavitation dynamics on a standard hydrofoil with a compressible approach and a delayeddetached eddy simulation and show evidence of pressure wave fronts associated with a cavity cloud collapse. Pancirolli et al. [7] use a smooth particle hydrodynamics (SPH) model to study the cavity formation during the impact of a 2D wedge on a free surface and show that there is a threshold for the ratio between horizontal and vertical velocities above which a cavity is created, greatly impacting the stability of the body during slamming.

On a more applied point of view, Podnar et al. [8] improved the shape of an hydrofoil in order to reduce the cavitation development in a bulb turbine. Decaix et al. [9] conducted a comprehensive investigation on the importance of the time step tuning to achieve a good

**Citation:** Ravelet, F. Editorial for Special Issue: New Advances of Cavitation Instabilities. *Appl. Sci.* **2021**, *11*, 5313. https://doi.org/10.3390/ app11125313

Received: 1 June 2021 Accepted: 4 June 2021 Published: 8 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

prediction on the cavitating behavior of a Francis Turbine. Finally, Li et al. [10] numerically explored the shape optimization of a globe valve in order to suppress cavitation inside of it.

**Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


## *Article* **Peening Natural Aging of Aluminum Alloy by Ultra-High-Temperature and High-Pressure Cavitation**

**Toshihiko Yoshimura 1,\*, Masayoshi Iwamoto <sup>1</sup> , Takayuki Ogi <sup>1</sup> , Fumihiro Kato <sup>1</sup> , Masataka Ijiri <sup>2</sup> and Shoichi Kikuchi <sup>3</sup>**


#### **Featured Application: Mechanical structural parts such as automobile parts made of aluminum alloy that have long-term stable compressive residual stress and require an extremely hard surface.**

**Abstract:** The peening solution treatment was performed on AC4CH aluminum alloy by ultra-hightemperature and high-pressure cavitation (UTPC) processing, and the peening natural aging was examined. Furthermore, peening artificial aging treatment by low-temperature and low-pressure cavitation (LTPC) was performed, and the time course of peening natural aging and peening artificial aging were compared and investigated. It was found that when the AC4CH alloy is processed for an appropriate time by UTPC processing, compressive residual stress is applied and natural aging occurs. In addition, the UTPC processing conditions for peening natural aging treatment with high compressive residual stress and surface hardness were clarified. After peening artificial aging by LTPC processing, the compressive residual stress decreases slightly over time, but the compression residual stress becomes constant by peening natural aging through UTPC treatment. In contrast, it was found that neither natural nor artificial peening natural aging occurs after processing for a short time.

**Keywords:** multifunction cavitation; water jet cavitation; ultrasonic cavitation; high-temperature high-pressure cavitation; peening natural aging; low-temperature low-pressure cavitation; peening aging

#### **1. Introduction**

The mechanism of high-pressure water jet action on materials has been well understood for decades [1]. Shot peening [2,3] and fine particle peening [4] are techniques for mechanically striking surfaces with particles to impart compressive residual stress. Another processing technology modifies surfaces with an abrasive water jet formed by adding an abrasive to water pressurized to about 300 MPa [5,6].

When water jet cavitation (WJC) is irradiated with ultrasonic waves to increase the temperature and pressure, the cavitation has both mechanical and electrochemical affects. This is called multifunction cavitation [7,8] (MFC, PCT Pub. No. W02016/136656, US Registered Patent [9]). WJC has a bubble radius of about 100 µm, and the impact pressure at the time of bubble collapse is about 1000 MPa [10]. In contrast, in the ultrasonic cavitation (UC) used for cleaning, the bubble nuclei originally existing in the water with a radius of several micrometers are irradiated with ultrasonic waves, and the collapse pressure is small, but the temperature at the time of bubble shrinkage is several thousand kelvins [11–13]. MFC bubbles have a large radius of about 100 µm, a collapse temperature of several

**Citation:** Yoshimura, T.; Iwamoto, M.; Ogi, T.; Kato, F.; Ijiri, M.; Kikuchi, S. Peening Natural Aging of Aluminum Alloy by Ultra-High-Temperature and High-Pressure Cavitation. *Appl. Sci.* **2021**, *11*, 2894. https://doi.org/ 10.3390/app11072894

Academic Editors: Marwan Al-Haik and Theodore Matikas

Received: 6 February 2021 Accepted: 22 March 2021 Published: 24 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

thousand kelvins, and a high collapse pressure of about 1000 MPa. Such bubbles are a feature of WJC and UC. The distribution and unsteady flow of WJC bubble clouds that are changed by sound pressure from incident ultrasonic waves have been clarified analytically and experimentally [14]. Previous studies have developed ultra-high-temperature and highpressure cavitation (UTPC) obtained by attaching a swirl flow nozzle (SFN) used for water jet peening in air to further increase the bubble size, number of bubbles, temperature, and pressure of MFC [15]. In addition, a photon counting head was used to measure the photons (sonoluminescence) emitted from UTPC, MFC, UC, SFN-WJC, and WJC. The UTPC has the highest emission intensity, number of photons, and bubble interior temperature upon collapse. It was clarified that there is a correlation between the sonoluminescence results and the surface modification results of low-alloy steel [9]. Using these high-temperature and high-pressure bubbles, surface modifications of various metals and ceramics have been carried out.

The foremost use of aluminium casting alloys is in the automotive industry (60% of the tonnage world-wide). These are mainly engine components elaborated from alloys with silicon, most of which also contain at least 3% copper. The main applications outside the automotive industry are in mechanical construction, electrical engineering, transport, household electrical appliances and ironmongery [16]. Studies on the processing and natural aging of aluminum alloys have been also reported. The influence of the stirring pin and pressing tool shoulder on the microstructural softening during friction-stir processing and subsequent natural aging behavior was investigated for a 6061-T6 aluminum alloy [17]. Severe plastic deformation (SPD) can be used to generate ultra-fine-grained microstructures and thus to increase the strength of many materials. Unfortunately, highstrength aluminum alloys are generally hard to deform, which imposes severe limits on the feasibility of conventional SPD methods. The low-temperature equal-channel angular pressing method was used to deform an AA7075 alloy, and the influence of natural aging of a solid solution heat treated AA7075 alloy (prior to SPD) on the microstructure and the mechanical properties after equal-channel angular pressing was studied [18].

In this study, the outermost surface of a material was solution-treated using UTPC processing to raise the temperature to at least 500 ◦C higher than the equilibrium phase diagram of the aluminum alloy. Then, the elements added to the aluminum alloy formed a supersaturated solid solution and became distributed uniformly. After that, the possibility of T4 treatment, which is a natural aging process that gradually precipitates a Guinier– Preston (GP) zone at room temperature, was examined from the viewpoints of hardness, structural change, corrosion resistance, surface shape, and changes in compressive residual stress over time.

#### **2. Materials and Methods**

The specimens were Al-Mg-Si alloy (AC4CH), and its chemical composition is shown in Table 1. Rolled-plate tempered material of 30 mm × 30 mm × 5 mm was cut from the rolled state and used as specimens. Therefore, they were in a state where work hardening remained.

**Table 1.** Chemical composition of Al-Mg-Si alloy (AC4CH) (wt%).

UTPC processing aimed at natural aging and low-temperature and low-pressure cavitation (LTPC) processing for comparison were performed on the rolled tempered material. Figure 1 is a schematic diagram of the equipment that performed UTPC and LTPC processing. High-pressure water is injected from a high-pressure pump and generates WJC, which is irradiated with ultrasonic waves from an ultrasonic transducer. The isothermal

expansion and adiabatic compression [19,20] of WJC continues when the sound pressure

**Si Mg Fe Ni Ti Al** 6.80 0.31 0.10 0.01 0.01 Bal.

around WJC remains below the break threshold [21] due to changes in the sound pressure of longitudinal waves. This produces high-temperature and high-pressure MFC. Furthermore, when a swirling nozzle is attached to the jet nozzle, the pressure of the jet center is reduced as the swirling flow is generated, the size of WJC increases, and the number of bubbles of WJC also increases due to the decrease of the cavitation number. When this WJC is irradiated with ultrasonic waves, UTPC is generated [15].

**Figure 1.** Schematic diagram of processing equipment for ultra-high-temperature and high-pressure cavitation (UTPC) and low-temperature and low-pressure cavitation (LTPC).

μ μ μ μ As the study on the cavitation number, there is a report that the formation of cloud cavitation vortex is controlled by the cavitation number in a narrow region similar to the swirling nozzle of this study [22]. According to the calculation of a single bubble using the Rayleigh-Plesset equation, for example, when a sound pressure of 130 kPa due to ultrasonic waves is applied to a WJC with a radius of 100 µm, the bubble expands to the bubble radius of 200 µm and shrinks to 3 µm in a time of 41 µs. The pressure reaches 3700 MPa and the temperature becomes 2.37 <sup>×</sup> <sup>10</sup><sup>5</sup> K inside the bubble. However, in fact, at the final stage of contraction, the increase in the bubble temperature is suppressed by chemical reaction heat and heat conduction accompanying the thermal decomposition of steam. In conventional measurement data of sonoluminescence, it is considered that the maximum bubble temperature is 100,000 K. The dynamics of single-bubble and multi-bubble was experimentally investigated using sonoluminescence [23]. It is known that the pressure and temperature achieved inside a bubble during multi-bubble sonoluminescence (MBSL) are much lower than those during single-bubble sonoluminescence (SBSL). In addition, when MFC or UTPC bubbles reach the minimum volume, they flatten in the direction parallel to the material surface, and the bubble part farthest from the surface becomes a piercing form, forming a high-speed jet (microjet) that pierces the material surface. However, the instability in microjet and two bubble interaction were reported [24]. Further studies are needed on the relationship between interactions in multi-bubble and microjet.

The pressure inside a bubble during collapse is mainly determined by the pump discharge pressure, and the temperature inside the bubble is mainly controlled by the ultrasonic output. The UTPC conditions were pump discharge pressure of 35 MPa, the standoff distance of 65 mm between the sample and the water jet nozzle, ultrasonic output of 800 W, and ultrasonic frequency of 28 kHz, whereas the LTPC conditions were pump

discharge pressure of 20 MPa, ultrasonic output of 100 W, and ultrasonic frequency of 28 kHz. Here, the ultrasonic mode was selected as the dual mode because it had the highest peening effect in previous studies. The LTPC condition is the so-called peening aging [25] condition, in which the highest compressive residual stress more than −170 MPa can be applied to the AC4CH aluminum alloy. The sound pressure of ultrasonic waves controls the temperature of MFC bubbles. In this study, the sound pressure of ultrasonic output 800W that irradiates the water jet cavitation is 15 mV under the UTPC condition. On the other hand, the sound pressure of the ultrasonic output 100 W is 5 mV under the LTPC condition. These are the values of relative sound pressure measured by the sound pressure sensor (Portable Sonic Monitor: HUS-3, HONDA ELECTRONICS CO., LTD.). The processing time was changed to 2 min, 5 min, 10 min, 20 min, and 30 min to increase the outermost surface temperature to the solution heat temperature and clarify the time required for the solution treatment. Furthermore, to investigate the effect of processing for a longer time, UTPC for 60 min was also performed. −

As a method for evaluating the peening natural aging by cavitation, Vickers microhardness measurement (Mitutoyo Corp., HM-210B system), observation and analysis by scanning electron microscope–energy dispersive X-ray spectroscopy (SEM-EDS) performed by a field emission SEM (Hitachi, S-4800), surface roughness measurement, and surface potential (corrosion potential) measurement with a Kelvin force probe microscope (KFM) (Hitachi, AFM5000II) were carried out. In the hardness measurement, the specific locations in the peening region to which the compressive residual stress was imparted were determined, and the change in hardness over time was examined. In addition, 10 points were measured at each processing time under each condition, and the average of 8 points with the maximum and minimum values removed was calculated. In addition, the compressive residual stress applied to the surface was measured by Cr Kα X-ray diffraction (MSF–3 M, Rigaku Corp.) to evaluate the change over time. α

Figure 2 shows the principle of cavitation natural aging. As previously discussed, UTPC (micro-forging) is suitable for surface modification of refractory metals, and LTPC is suitable for low-melting-point metals, such as aluminum alloys. When LTPC is applied to an aluminum alloy, the microjet has few hot spots [11] and the temperature inside the bubbles is low, which is appropriate for artificial aging, so it is possible to impart compressive residual stress and artificial aging simultaneously. However, when a lowmelting-point metal is processed by UTPC with many hot spots in the microjet, the additive elements of Mg and Si become a supersaturated solid solution at room temperature, but the outermost surface is overheated above the solution temperature. Therefore, Mg and Si are mixed uniformly. After the UTPC stops, a GP zone gradually precipitates due to natural aging, and the compressive residual stress formed by the elastic restraint of the surroundings creates new dislocations, which are entangled with the precipitates and hardened by the Orowan mechanism.

#### **3. Results and Discussion**

Figure 3 shows photographs of the LTPC processing surface jetted at a fixed point. When fixed-point injection is performed with using the SFN, an erosion region is formed at the injection center, and a peening region is formed around the transition region of the erosion region. In particular, the distribution of the peening region becomes pronounced after processing of 5 min or more. Figure 4 shows the photographs of UTPC processing surface jetted at a fixed point. Compared to LTPC, the peening area is expanded in the peripheral direction. Although the area of the perforated disc-shaped peening region of LTPC is small, the compressive residual stress is higher than that of UTPC, which has a large peening region, due to the peening artificial aging.

**Figure 3.** Surface photographs after LPTC processing (alloy: AC4CH, jet: 20 MPa, UC: 100 W).

○: ○: **Figure 4.** Surface photographs after UPTC processing (alloy: AC4CH, jet: 35 MPa, UC: 800 W). (#: roughness and hardness measurement area).

Figure 5 shows the changes of surface roughness with processing time in the peening region after UTPC and LTPC processing. Arithmetic mean roughness Ra increases with processing time, but there is no significant difference in surface roughness between UTPC and LTPC up to 20 min. However, from the KFM observation results, the number of peening marks is larger in LTPC than in UTPC, and the surface of LTPC is relatively smooth compared to that of UTPC, but the UTPC surface has many fine irregularities.

**Figure 5.** Change of surface roughness with processing time of LPTC and UTPC (LTPC: jet 20 MPa, UC 100 W; UTPC: jet 35 MPa, UC 800 W) (Evaluation length 5 mm)

The Vickers microhardness of the untreated material was 85.3 HV. Figure 6 shows the time course of hardness after UTPC processing at various times. No increase in hardness was observed after 2 min or 5 min of processing, but an increase in hardness was observed after processing for 10 min or more. The hardness of the material treated for 20 min increased the most. After 20 min of UTPC processing, the hardness reached its maximum value 21 days after processing, and then the hardness decreased and increased again.

**Figure 6.** Relationship between hardness and elapsed time after UPTC processing (WJ: 35 MPa, UC: 800 W) (Average of standard deviation: 2 min 5.9 HV, 5 min 7.5 HV, 10 min 9.9 HV, 20 min 9.5 HV, 30 min 11.2 HV, 60 min 6.9 HV).

The reason the hardness increased after 10 min or more of UTPC processing is as follows. The outermost surface is solution treated by UTPC processing, as shown in Figure 2, and the supersaturated solid solution elements Mg and Si are perfectly mixed in aluminum at the solution temperature. After that, when cavitation stops and the temperature decreases to room temperature, the surface becomes a supersaturated solid solution, and Mg and Si gradually precipitate as GP zones. The surface is plastically deformed by the microjet when cavitation bubbles collapse, and compressive residual stress is imparted. This compressive stress causes dislocations on the surface, and the dislocations are trapped in the precipitates and hardened by the Orowan mechanism. The reason the hardness of material processed for 20 min and 30 min decreases at first and then increases again is thought to be a result of natural aging as the structure changes from GP I phase to GP II phase and θ' phase. It is considered that the compressive residual stress generated causes new dislocations that are related to the pinning effect on the precipitate. θ

In contrast, it can be recognized that after 60 min of processing, unlike the material processed for 20 min, hardening does not occur immediately after processing. The hardness immediately after processing was about 160 HV, which was about the same as that for the LTPC 20-min processing and LTPC 30-min processing shown in Figure 7. However, after 60 days, the 60-min processed surface was hardened to 183 HV.

**Figure 7.** Relationship between hardness and elapsed time after LTPC processing. (WJ: 20 MPa, UC: 100 W) (Average of standard deviation: 2 min 7.7 HV, 5 min 6.6 HV, 10 min 8.5 HV, 20 min 8.4 HV, 30 min HV, 60 min 8.4 HV).

Al-0.62% Mg-0.39% Si alloy was quenched in water after solution treatment and naturally aged. The hardness increased by 12 HV on 1 day after solution treatment and quenching, and an increase by 22 HV on day 11 has been reported [26]. Because the chemical composition of Al-Mg-Si is different between these two cases, a simple comparison cannot be made, but in the UTPC processing of this study, after the solution temperature was reached and the cavitation stopped, specimens were rapidly cooled in water, similar to quenching. However, the decisive difference between natural aging by heat treatment and peening natural aging by UTPC is the presence of compressive residual stress after quenching of peening natural aging on the surface. When a GP zone is formed by natural aging, dislocations occur due to compressive residual stress, and the hardness increase is larger than that of natural aging by heat treatment. In fact, as shown in Figure 6 the hardness of the material processed by UTPC for 20 min increased by about 50 HV from 1 day to 9 days later.

The compressive residual stress after LTPC processing (20 MPa, 100 W, 20 min) was −170.6 MPa, and the compressive residual stress after UTPC processing (35 MPa, 800 W, 20 min) was −133.67 MPa [25]. Figure 7 shows the change in hardness over time after LTPC processing (20 MPa, 100 W, 20 min), which gives the highest compressive residual stress. Unlike the case of UTPC, hardness did not increase over time. This is because the surface did not reach the solution temperature during LTPC processing, and peening natural aging did not occur. Similar to the results for UTPC processing in Figure 6, the material hardness of each LTPC processing time tended to decrease once and increase again. This is the result of a balance between changes in precipitates and the formation of new dislocations due to compressive stress after peening artificial aging, which has the highest compressive residual stress during processing, although it has not reached the solution temperature.

A qualitative SEM-EDS line analysis of the region of interest (ROI (Regions of Interest): the function that creates a window by specifying the X-ray energy range of interest in the spectrum) on the surface after UTPC and LTPC processing was carried out. Here, many-line qualitative analysis at 100,000 magnification was performed on the surface of specimens 1 day and 6 days after UTPC processing. Figure 8 shows the frequency distribution of the ROI analysis intensity of each Al, Mg, and Si concentration. The horizontal axis is the analytical intensity of various elements, and the vertical axis is the frequency distribution of the analytical intensity, showing the concentration distribution of various elements. It can also be recognized that the Mg concentration and Si concentration distribution were broader after 6 days in comparison with that after 1 day. When a GP zone is formed due to fluctuations in the concentrations of solute elements of Mg and Si by spinodal decomposition, the concentration of solvent elements changes to the same extent that the concentration of solute elements fluctuates. However, the fact that the spread of the Mg concentration distributions after 6 days of UTPC and 1 day of LTPC in Figure 8b are large, and the spread of solvent element Al concentration distributions are remarkably large in Figure 8a. This indicates that phase separation occurs due to spinodal decomposition in natural aging, which is consistent with the result of the increase in hardness in Figure 6. Furthermore, according to the result of qualitative analysis of ROI by SEM-EDS line analysis of the surface of a specimen 1 day after LTPC, the surface temperature did not reach the solution temperature. Therefore, phase separation occurred as it does without processing, and the solvent element Al concentration distribution is as wide as that in the specimen 6 days after UTPC.

Figure 9a,b shows the shape image and corrosion potential image of the mirrorpolished surface before processing, respectively, as measured by KFM. The shape image is flat due to mirror polishing, and the potential image has no bias and a uniform potential distribution.

Figure 9c,d shows the KFM measurement results 72 days after UTPC processing of a specimen for 20 min; at this point, the hardness was starting to increase again. Although unevenness after processing is observed in the shape image, the potential image does not have a particularly biased distribution. In addition, the surface was leveled by raising the surface temperature to a level that generated solution treatment, and no remarkable peening marks were observed. Figure 9e,f provides the KFM measurement results after 104 days, when the hardness gradually increased after UTPC (WJ: 35 MPa, UC: 800 W, 60 min) processing. Although large irregularities are observed in the shape image, the potential image does not have a particularly biased distribution. In contrast, Figure 9g,h shows the KFM measurement results 7 days after LTPC processing of a specimen, which applies the highest compressive residual stress. A dent, which is considered to be a peening mark, is observed in the shape image. The surface potential at the peening mark is lower than that in the other regions. Further, the reason the peening marks were observed is that the surface did not reach the solution temperature and surface homogenization did not occur due to the lower processing temperature.

**Figure 8.** Frequency distribution of Mg and Si concentrations obtained from SEM-EDS line analysis 6 days after UTPC processing: (**a**) Al; (**b**) Mg; (**c**) Si. (WJ: 35 MPa, UC: 800 W, 20 min).

Table 2 gives the surface potential results for the various specimens. KFM can measure the work function, which is a measure of the ease with which electrons are emitted from the surface. Because oxidation involves the transfer of electrons, the higher the surface potential, the higher the corrosion resistance. Among aluminum alloys, AC4CH without processing has high corrosion resistance, and its surface potential after mirror polishing is as high as 1200 mV. The surface potential of special steel does not exceed 1000 mV, even after MFC processing [27]. Even if the aluminum alloy is mirror-polished, a passivation film is immediately formed and corrosion resistance is maintained. The surface potential of the specimen subjected to LPTC processing, which imparts the highest compressive residual stress, is almost the same as that of the as-received specimen after mirror polishing. UPTC processing, which causes peening natural aging, produced a decreased surface potential with 20 min of processing and increased potential with 30 min and 60 min of processing. After the outermost surface reaches the solution temperature and is subjected to the solution treatment by UTPC, precipitates such as GP zones are formed on the surface. The disparate surface potential results for the UTPC specimens are probably due to different formation states of these precipitates being created by different processing times (20 min vs. 30 min and 60 min). In the 60-min processed material, no change in surface potential was observed between 41 days and 104 days.

**Table 2.** Surface potential of aluminum alloy (AC4CH) by various processing.


**Figure 9.** Results of shape imaging and corrosion potential measurement of UTPC and LTPC processed material by Kelvin force probe microscope (KFM).

In the UTPC 20-min processing and 30-min processing of the present study, the initial surface roughness was small (Ra: 0.1 µm); therefore, the surface was heated to the solution temperature, and the peening solution process was possible. In contrast, in UTPC 60-min processing, peening solution processing was performed up to 30 min; however, as the surface roughness increases with processing, the microjet impingement is not always vertical, and despite processing under UTPC conditions, UTPC processing does not occur. This is because the processing was performed at low temperature and low pressure. In addition, as can be seen from the shape image of Figure 9e,f due to long-term UTPC processing, high-temperature processing was not possible because the specific surface area increased and the ratio of heat dissipation to water increased rather than the heat input to the alloy interior. Therefore, it is probable that in the latter half of the 60-min processing, the peening artificial aging was performed by the pseudo-LTPC processing from the state where the solution was once processed. Therefore, as shown in Figure 6, hardening by natural aging did not occur after processing.

In UTPC and LTPC processing, the temperature of a metal surface is determined by the balance among heat input to the surface from a large number of cavitation bubbles with hot spots, convective heat dissipation from the surface to water, and heat dissipation due to heat conduction inside the metal. However, it is difficult to measure directly. Therefore, it is carried out by estimating the bubble temperature from the emission intensity (measurement of the number of photons per unit time) at the time of bubble collapse [9] and the temperature estimation is carried out by the structure change based on the equilibrium phase diagram after processing. For example, when UTPC processing is performed on Cr-Mo steel, it has been confirmed that spherical cementite that changes its structure at 500 ◦C or higher is formed just below the surface [28]. Therefore, it is considered that the surface temperature of AC4CH has reached the solution temperature in the UTPC treatment of this study.

Figure 10 shows the values of compressive residual stress applied to the surface immediately after UTPC and LTPC processing and after a long period of time. The compressive residual stress after 7 days of LPTC 20-min processed material is the highest at −171 MPa; however, after 141 days, it drops to the same level as UTPC 20-min processed material. However, the compressive residual stress of the UTPC 20-min processed material after 7 days is −134 MPa, which is not higher than that of the LPTC 20-min processed material, but the compressive residual stress does not decrease much even after 175 days. The compressive residual stress is reduced to −100 MPa or less in the UTPC 30-min processed material and UTPC 60-min processed material that have passed for a long time. The stability of the compressive residual stress of the UTPC peening natural aging material (UTPC 20 min) is superior to that of the LTPC peening artificial aging material (UTPC 20 min). With the artificial peening aging by LTPC, a compressive residual stress of −100 MPa or more can be applied even in a short time of 2 min, and as shown in Figure 7, a hardness of 120 to 140 HV can be obtained. Considering the decrease of compressive residual stress applied by the peening artificial aging for a long time, it can be judged that short-time processing is sufficient for peening artificial aging treatment by LTPC. Regarding the surface hardness, both the peening natural aging material by UTPC and the peening artificial aging material by LTPC reached almost 180 HV and achieved the same level of hardness after 100 days or more with longer-term processing of 20 min or more. From the viewpoint of corrosion resistance, the peening naturally aged material (UTPC 30 min) by UTPC is superior. From the above results, it is necessary to properly use peening natural aging and peening artificial aging depending on the intended use of the AC4CH alloy.

**Figure 10.** Aging variation of compressive residual stress of LTPC and UTPC.

#### **4. Conclusions**

Peening solution treatment was performed on AC4CH by UTPC processing to investigate peening natural aging. Furthermore, peening artificial aging treatment by LTPC was carried out to examine the change over time of peening natural aging and peening artificial aging, and the following findings were obtained:

(1) When the AC4CH alloy is processed for an appropriate time by UTPC processing, compressive residual stress is applied and natural aging occurs.

(2) The UTPC processing conditions of peening natural aging treatment to create a high surface compressive residual stress, surface hardness, and corrosion resistance were clarified.

(3) After artificial aging by LTPC treatment, the compressive residual stress decreases slightly with the passage of time, but the compressive residual stress is stabilized by peening natural aging by UTPC treatment.

(4) Processing for a short time does not cause peening natural aging by UTPC.

**Author Contributions:** All authors contributed to the paper either during the writing or editing phases. All authors made a substantial contribution to this research. Planning of research method and evaluation method, T.Y., M.I. (Masataka Ijiri) and S.K.; cavitation processing, M.I. (Masayoshi Iwamoto), T.Y. and M.I. (Masataka Ijiri); hardness measurement and surface roughness measurement, M.I. (Masayoshi Iwamoto); KFM measurement, T.O. and F.K.; measurement of compression residual stress, S.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was supported in part by JSPS KAKENHI Grant Number 19K04110 (Grantin-Aid for Scientific Research (C)) and by the Light Metal Educational Foundation, Inc.

**Data Availability Statement:** Data is contained within the article. The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Review* **Cavitating Jet: A Review**

#### **Hitoshi Soyama**

Department of Finemechanics, Tohoku University, Sendai 980-8579, Japan; soyama@mm.mech.tohoku.ac.jp; Tel.: +81-22-795-6891; Fax: +81-22-795-3758

Received: 15 September 2020; Accepted: 15 October 2020; Published: 17 October 2020

#### **Featured Application: Cavitation peening, Cleaning, Drilling.**

**Abstract:** When a high-speed water jet is injected into water through a nozzle, cavitation is generated in the nozzle and/or shear layer around the jet. A jet with cavitation is called a "cavitating jet". When the cavitating jet is injected into a surface, cavitation is collapsed, producing impacts. Although cavitation impacts are harmful to hydraulic machinery, impacts produced by cavitating jets are utilized for cleaning, drilling and cavitation peening, which is a mechanical surface treatment to improve the fatigue strength of metallic materials in the same way as shot peening. When a cavitating jet is optimized, the peening intensity of the cavitating jet is larger than that of water jet peening, in which water column impacts are used. In order to optimize the cavitating jet, an understanding of the instabilities of the cavitating jet is required. In the present review, the unsteady behavior of vortex cavitation is visualized, and key parameters such as injection pressure, cavitation number and sound velocity in cavitating flow field are discussed, then the estimation methods of the aggressive intensity of the jet are summarized.

**Keywords:** cavitation; jet; vortex; mechanical surface treatment; cavitation peening

#### **1. Introduction**

Cavitation is a harmful phenomenon for hydraulic machineries such as pumps, as severe impacts are produced at bubble collapse [1,2]. However, cavitation impacts are utilized for mechanical surface treatment in the same way as shot peening, and this is named "cavitation peening" [3,4]. The great advantage of cavitation peening is that shots are not used in the peening process, as cavitation impacts are used in cavitation peening [5]. Thus, the cavitation-peened surface is less rough compared with the shot-peened surface, and the fatigue strength of cavitation peening is better than that of shot-peening [6]. In conventional cavitation peening, cavitation is generated by injecting high-speed water jet into water [3,4], and a submerged water jet with cavitation is called a "cavitating jet". The cavitation peening is utilized for the impacts of cavitation collapses, and it is different from water jet peening, in which water column impacts are used. To use the cavitating jet for peening, it is worth understanding the mechanism of the cavitating jet.

In a cavitating jet, cavitation is produced inside and/or outside a nozzle when sufficient pressure difference is applied across the nozzle. As Monkbadi reviewed, vortex structures in turbulent jets [7], ring vortices (0-mode), a single helical vortex (1st mode) and double-helical vortices (2nd mode) were observed in a cavitating jet [8–10]. In the case of a developed cavitating jet, which has been used for practical applications such as cutting, material testing, drilling and peening [11–17], cavitation clouds are shed periodically [16,18–27]. As it was reported that the lifetime of the cloud is a key factor in the aggressive intensity of the cavitating jet [28], the investigation of cavitation cloud shedding is very important. On the other hand, from the viewpoint of cavitation inception, a turbulent jet with cavitation was investigated [29–32]. In the present review, in order for the cavitating jet to be used

for practical applications, the main subject is the developed cavitating jet, as shown in Figure 1a. In Figure 1, white bubbles are cavitation bubbles, as a used flash lamp was placed at the same side of a camera. As shown in Figure 1a, the cavitation clouds are clearly observed in the cavitating jet in water.

**Figure 1.** Typical aspects of cavitating jet: (**a**) cavitating jet in water; (**b**) cavitating jet in air.

Normally, a cavitating jet is produced by injecting a high-speed water jet into a water-filled chamber, as mentioned above. To apply a cavitating jet for components and/or plants which cannot be put in the chamber, Soyama developed a cavitating jet in air by injecting a high-speed water jet into a low-speed water jet, which was injected into air without the water-filled chamber using a concentric nozzle [33–36]. A typical cavitating jet in air by injecting a high-speed water jet into a low-speed water jet is shown in Figure 1b. Cavitation clouds are also observed in the water column of low-speed water jets of cavitating jets in air, as shown in Figure 1b. Even though the cavitating jet was in air, cloud shedding was observed [35,37,38]. In addition, in the case of the cavitating jet in air at optimum conditions, the wavy pattern of the low-speed water jet was observed, as shown in Figure 1b [34,35,37,38]. The frequency of the wavy pattern is equal to the shedding frequency of the cloud [35]. Even for the cavitating jet in air, unsteady behavior is very important.

In cavitating jets both in water and in air, whereas the cloud shedding frequency is several hundred Hz [18,35,39,40], only a few severe impacts occur per second [34,41–43] when cavitation impacts are measured by special-made PVDF transducers [41,42]. Thus, in order to enhance the aggressive intensity of the cavitating jet in water and air for practical applications of cavitating jets, the unsteady behavior of the cavitating jet should be investigated. In the present review, the normal cavitating jet, i.e., the cavitating jet in water, was mainly discussed, and the cavitating jet in water was simply described as the cavitating jet.

To simulate the cavitating jet numerically, it is important to understand the flow pattern of the cavitating jet. As is well known, numerical simulation of cavitating flow is not yet easy, due to the high Reynolds number and phase change phenomenon. Numerical investigation of three-dimensional cloud cavitation with special emphasis on collapse induced shock dynamics [44], simulation of cloud cavitation on propeller [45], and the hydrofoil [46] were carried out. In the area of numerical simulation of the cavitating jet, considering bubble dynamics, the numerical simulation of vortex cavitation in a three-dimensional submerged transitional jet near inception condition was carried out [47], and bubble growth in the shear layer of the cavitating jet was calculated [48]; residual stresses introduced by cavitating jet considering bubble growth and collapse were also simulated [49]. The cavitating flow in a venturi nozzle was also tried by a large eddy simulation of turbulence–cavitation interactions [50]. From the viewpoint of cloud-shedding of the cavitating jet, numerical analysis of cavitation cloud-shedding in a free submerged water jet was carried out [51,52]. The interaction of cavitation bubbles and materials was also investigated numerically [53–56]. In the case of the aggressive cavitating jet, the clouds shed periodically and form ring vortex cavitation on the impinging surface, and then collapse, producing severe impacts. Thus, the experimental investigation of the flow pattern would be assisted by numerical simulations for more precise investigations of the cavitation impacts produced by the cavitating jet.

In the present review, in order to use the cavitating jet for practical applications, papers that investigated the flow pattern of the cavitating jet experimentally were reviewed. To enhance the aggressive intensity of the cavitating jet for these applications, the key factors of the cavitating jet were also summarized, and the estimation method of the aggressive intensity of the jet was discussed. In the present paper, the aggressive intensity of the cavitating jet means erosion rate measured by weight loss of the target metals and/or the peening intensity measured by the arc height of the metallic plate.

#### **2. Cavitation**

Cavitation is a phase-change phenomenon in which the liquid phase is changed to a gas phase due to a decrease in liquid pressure to vapor pressure by increasing flow velocity [2], and it is called "hydrodynamic cavitation". Cavitation is also generated by ultrasonic vibration, which is named "ultrasonic cavitation". Bubbles are also formed by irradiating a pulse laser into water, in which the bubble behaves as a cavitation bubble [57]; this is called "laser cavitation". A schematic diagram of cavitation is illustrated in Figure 2a. When a cavitation nucleus such as a tiny air bubble is subjected to a low-pressure region, it becomes a cavitation, and it develops and shrinks, then collapses, producing microjet and shock wave during rebound. The microjet and the shock wave produce a severe impact, which can deform metals. After rebound and shrinking, bubbles remain in the water, which is called a residual bubble [58]. In the research area of cavitation, including bubble dynamics, spherical bubbles have mainly been investigated by numerical simulations [59] and experimental studies [57,60,61], and the effects of bubble shape and bubble interactions have also been studied [62].

**Figure 2.** Schematic diagram of development and collapse of cavitation: (**a**) spherical bubble; (**b**) vortex cavitation.

In experimental studies of severe cavitation erosion, it was found that vortex cavitation, as shown in Figure 3 [63], produced severe impacts on fluid machineries such as pumps and valves [1,64,65]. A schematic diagram of vortex cavitation is shown in Figure 2b. Cavitation nuclei accumulate in a vortex, and they become vortex cavitation in a high-speed region, i.e., a low-pressure region. The vortex cavitation develops and shrinks, then collapses. Regarding a model test of vortex cavitation using a rotating device [66], a microjet was observed in the vortex cavitation. Thus, in order to use cavitation impacts for practical applications, the generation of vortex cavitation is very important.

**Figure 3.** Typical aspect of vortex cavitation, which produces severe impact [63].

#### **3. Cavitating Jet**

#### *3.1. Structure of Cavitating Jet*

To show the flow pattern of the cavitating jet, Figure 4 illustrates aspects of the impinging cavitating jet taken with (a) a flash lamp of 1.1 µs and (b) a shutter speed of 1/25 sec, i.e., 40 ms. Thus, Figure 4a reveals an instantaneous aspect of the jet, and Figure 4b shows a kind of time-averaged aspect of the impinging cavitating jet. In the nozzle used for the cavitation jet, the cavitator and the guide pipe were installed to enhance the aggressive intensity of the cavitating jet [39], and standoff distance was defined as the distance from the upstream corner of the nozzle to the specimen surface. As shown in Figure 4a, a cloud cavitation was observed between the nozzle and the impinging surface, and a ring vortex cavitation was observed on the surface. When the cavitating jet was observed by a normal light, as shown in Figure 4b, the cavitating jet seems to be a continuous jet from a nozzle to an impinging surface. However, cloud cavitation sheds from the nozzle to the surface, and the cloud cavitation becomes a ring vortex cavitation.

**Figure 4.** Aspects of cavitating jet (cylindrical nozzle, nozzle diameter *d* = 2 mm, upstream pressure *p*<sup>1</sup> = 30 MPa, downstream pressure *p*<sup>2</sup> = 0.1 MPa, standoff distance *s* = 222 mm): (**a**) observation with flash lamp whose exposure time is 1.1 µs; (**b**) observation with normal light and shutter speed of 40 ms.

To show vortex cavitation in the cavitating jet more precisely, Figure 5a reveals the aspect of the free cavitating jet through a conical nozzle, and Figure 5b shows the aspect of the impinging cavitating jet with a flash lamp [9]. As shown in Figure 5b, the PVDF transducers [41,42] were installed in the impinging surface, which the white arrow shows, and the signal from the PVDF transducer was synchronized to the flash lamp; thus, Figure 5b reveals the aspect of the jet which produces the impact on the impinging surface. As shown in Figure 5a, helical vortex cavitation is observed near the nozzle outlet, and cloud cavitation is observed downstream from the nozzle. Thus, the cloud cavitation results from merging the vortex cavitation. In the case of the impinging cavitating jet, a part of the ring cavitation collapses on the surface producing the impact, which was detected by the PVDF transducer. The ring vortex cavitation produces a severe impact on the impinging surface. In view of the practical applications of cavitation impacts generated by the cavitating jet, the collapse of the ring vortex cavitation is an important phenomenon.

**Figure 5.** Vortex cavitation in cavitating jet: (**a**) helical vortex cavitation in free cavitating jet (conical nozzle, nozzle diameter *d* = 2.1 mm, upstream pressure *p*<sup>1</sup> = 6 MPa, downstream pressure *p*<sup>2</sup> = 0.18 MPa); (**b**) ring vortex cavitation on impinging surface (cylindrical nozzle, nozzle diameter *d* = 2 mm, upstream pressure *p*<sup>1</sup> = 12 MPa, downstream pressure *p*<sup>2</sup> = 0.36 MPa, standoff distance *s* = 40 mm) [9].

(**a**) (**b**)

To reveal periodical shedding of cloud cavitation from the nozzle, as mentioned in the introduction, Figure 6 shows the aspect of the impinging cavitating jet taken by a high-speed video camera. The nozzle shown in Figure 6 was the same nozzle as in Figure 4, and it had the cavitator and the guide pipe. At *t* = 0 ms, the cavitation cloud sheds from the nozzle to the downstream, and the cloud reaches the impinging surface at *t* = 1.5 ms. A part of cloud cavitation becomes a ring vortex cavitation at *t* = 1.75 ms, and it spreads out on the surface and then collapses. Although the water jet is injected continuously into water, the vortex cavitations near nozzle shed downstream coalescing each other and they become a large cloud cavitation, thus the bubble density between clouds is reduced. At *t* = 4.25 ms, the new cavitation cloud, whose shape is similar to that at *t* = 0 ms, sheds from the nozzle. Thus, it is a periodical phenomenon, as previously reported [16,18–27].

According to a previous report [28], as shown in Figure 7, some cloud, i.e., the 1st cloud in Figure 7a, stays near the nozzle, and the 2nd jet core passes through the 1st cloud, as shown in Figure 7b,c; then, the 2nd jet core produces the 2nd cloud downstream of the 1st cloud. The 3rd jet core passes through the 1st and 2nd clouds and produces the 3rd cloud downstream of the 2nd cloud, as shown in Figure 7d. The 5th jet core produces the 5th cloud downstream of the 1st cloud as shown in Figure 7e. Detailed images of the free cavitating jet observed by the high-speed video camera and the schematic cross-sectional diagrams of the progress of a cavitating jet are shown in reference [28]. In Figure 6, the cavitation cloud marked by the yellow arrow stays at nearly the same position. The cloud cavitation staying where it is means that it has a long lifetime. The jet core with a longer lifetime cavitation cloud impinges the surface with high impingement pressure, as a drag of the jet core with the larger cloud is lesser than that of a smaller cloud because of the density of water and bubbles. Thus, the aggressive intensity of the cavitating jet with a longer lifetime cavitation cloud is larger.

Namely, the lifetime of the cloud is a key factor of the cavitating jet, and the longer the cloud lifetime, the more aggressive the intensity of the cavitating jet [28].

Figure 8 shows a typical aspect of a pure aluminum specimen that was exposed to the fixed cavitating jet, revealing the treatment area of the fixed cavitating jet. Plastic deformation pits are observed in a ring region, whose outer diameter is 60 mm and inner diameter is 30 mm. When the nozzle is scanned or the specimen is moved, the treatment area is uniform [67]. In Figure 8, as the nozzle throat diameter was 2 mm, the cavitating jet can treat an area 30 times wider than that of the nozzle throat. While the white region, which was impinged by the jet core, was observed at the jet center, the main treatment area is the ring region. Namely, the jet center is not treated by cavitation impacts. The understanding of the treatment area by the cavitating jet shows that a ring region is very important when using the cavitating jet for the practical applications.

**Figure 6.** Aspects of impinging cavitating jet observed by high-speed video camera (cylindrical nozzle, nozzle diameter *d* = 2 mm, upstream pressure *p*<sup>1</sup> = 30 MPa, downstream pressure *p*<sup>2</sup> = 0.1 MPa, standoff distance *s* = 222 mm).

**Figure 7.** Schematic diagram of the development of cavitation cloud in cavitating jet: (**a**) 1st cloud produced by the 1st jet core; (**b**) 2nd jet core passing through the 1st cloud; (**c**) 2nd cloud produced by the 2nd jet core and 3rd jet core passing through the 1st cloud; (**d**) 3rd cloud produced by the 3rd jet core and 4th jet core passing through the 1st and 2nd clouds; (**e**) 4th cloud produced by the 4th jet core, and 5th cloud produced by 5th jet core passing through the 1st clouds.

**Figure 8.** Typical treatment area by a fixed cavitating jet (pure aluminum, nozzle diameter *d* = 2 mm, upstream pressure *p*<sup>1</sup> = 30 MPa, downstream pressure *p*<sup>2</sup> = 0.1 MPa, standoff distance *s* = 262 mm, exposure time *t* = 1 min) [4].

Regarding the reason for the ring region, Figure 9 illustrates the schematic diagram of the impinging cavitating jet, considering the observations of the cavitating jet by the instantaneous photograph and the high-speed video. As shown in Figure 4b, when the cavitating jet is observed by an instantaneous photograph with normal light, the cavitating region seems to be a continuous jet, as shown in Figure 9a. It was previously thought that swirl cavitation in the shear region around the jet directly hits the impinging surface and that the swirl directly produces a ring treatment area [68–71]. However, this is incorrect, because the aspect of the cavitating jet, as shown in Figure 4b is a kind of time-averaged cavitating region. Considering Figures 4a, 5 and 6, vortex cavitation is initiated in and near the nozzle outlet, and cloud cavitations combine each other; then, the cloud cavitation forms the ring vortex cavitation on the impinging surface. Thus, in order to simulate bubbles in the cavitating jet, the pressure hysteresis of these processes should be considered. For example, if the cavitation in

swirl directly hit the impinging surface, the pressure around cavitation would be the pressure in the shear layer around the jet, and it gradually increases with the shedding of cavitation, as the jet speed decreases from the distance from the nozzle. On the other hand, the pressure of the cloud cavitation, which impinges the surface, increases at the impinging jet center suddenly and decreases in the ring vortex cavitation, then increases again.

**Figure 9.** Schematic diagram of cavitating jet: (**a**) swirl; (**b**) vortex, cloud and ring cavitation.

*σ* To consider the mechanism of the ring erosion occurring on the impinging surface, Figure 10 shows a schematic diagram of the local cavitation number on the surface [72]. The impinging pressure profile on the surface has a maximum at the jet center, *pmax*, and it changes with the injection pressure and cavitation number. The *b* in Figure 10 is a kind of jet width defined by the flow velocity [72]. The cavitation number of the cavitating jet, σ, is defined by the injection pressure, i.e., the upstream pressure of the nozzle *p*1, the downstream pressure of the nozzle *p*<sup>2</sup> and the vapor pressure of water *p<sup>v</sup>* as follows [2], and is simplified as Equation (1) and as *p*<sup>1</sup> >> *p*<sup>2</sup> >> *pv*.

$$
\sigma = \frac{p\_2 - p\_v}{p\_1 - p\_2} \approx \frac{p\_2}{p\_1} \tag{1}
$$

*σ −* ଵ ଶ ௫ ଶ In Figure 10, the local cavitation number σ*<sup>L</sup>* is defined by Equation (2) and is proportional to the ratio of *<sup>p</sup>* <sup>−</sup> *<sup>p</sup><sup>v</sup>* and <sup>1</sup> 2 ρ *vmax* 2 [72].

$$
\sigma\_L \propto \frac{p - p\_v}{\frac{1}{2} \rho |v\_{\text{max}}|^2} \propto \frac{p - p\_v}{p\_{\text{max}} - p\_2} f(r) \tag{2}
$$

*σ* ∂*σ* ∂ ∂*σ* ∂ ≈ Here, *p* and *vmax* are the pressure and the maximum flow velocity on the impinging surface. The *vmax* has a maximum at a certain distance from the jet center, *r*, as it is zero at the jet center and at a further point from the jet center as shown in the lower figure of Figure 10. As *vmax* is mainly determined by the pressure difference, *pmax* – *p*2, σ*<sup>L</sup>* is described by a function of *f*(*r*), as shown in the right-hand term of Equation (2). When the ring vortex cavitation sheds on the surface, the cavitation develops the ∂σ*L*/∂*r* < 0 region, and it collapses at ∂σ*L*/∂*r* ≈ 0. This is why the ring treatment area is obtained by the impinging cavitating jet. The detail of the pressure distribution on the impinging surface was shown in references [72,73].

**Figure 10.** Schematic diagram of local cavitation number [72].

#### *3.2. Periodical Shedding of Cavitation Cloud*

*σ σ σ σ σ σ* As mentioned above, the cavitation cloud sheds periodically [16,18–27], and it forms the ring vortex cavitation; then, the ring vortex cavitation collapses, producing impacts. Thus, the periodical shedding of the cavitation cloud is an important phenomenon. Figure 11 reveals aspects of periodical shedding of the cavitation cloud of a cavitating jet [21]. In Figure 11, the jet flows from the left-hand side to the right-hand side. As shown in Figure 11a, the cavitation cloud breaks near the nozzle at 0.5, 0.9, 1.4 and 1.9 ms; thus, the shedding frequency *fshedd* is about 2 kHz for *p*<sup>1</sup> = 20 MPa and σ = 0.014. When *p*<sup>1</sup> is increased to 30 MPa at a constant cavitation number, i.e., σ = 0.014, *fshedd* is increased to 2.4 KHz, as shown in Figure 11b. When σ is increased from 0.014 to 0.02 at constant *p*1= 20 MPa, *fshedd* also increases from 2 to 3.2 KHz, as shown in Figure 11a,c. *fshedd* is affected by not only *p*<sup>1</sup> but also σ. As is well known, the cavitating length *Lcav* and the width *wcav* are affected by *p*1, σ and *d*. In a previous report, *fshedd* was measured experimentally, changing with *p*1, σ and *d*, and the following relations were found [21].

$$f\_{shell} \propto w^{-1} \tag{3}$$

$$\oint\_{\text{shell}} \propto \mathcal{p}\_1^{0.45 \pm 0.03} \tag{4}$$

$$f \; f\_{\text{shell}} \; \propto \; d^{-0.98 \pm 0.14} \tag{5}$$

$$\int\_{\cdot} f\_{\text{shedd}} \propto \sigma^{0.83 \pm 0.10} \tag{6}$$

௦ௗௗ ∝ .଼ଷ±.ଵ U When the Strouhal number, *St*, is defined by *fsedd*, *w* and jet velocity at the nozzle exit, *U*, which is calculated from *p*1, the following experimental equation is obtained [21].

$$St \neq \frac{f\_{\text{shell}} \cdot w}{U} \approx 0.18 \pm 0.02\tag{7}$$

This result suggests that the cavitation cloud shedding is a phenomenon governed by a constant Strouhal number. On the other hand, the nozzle outlet geometry affects the aggressive intensity of the cavitating jet [39,74]. Details are described in Section 4.6. The *fsedd* is also affected by the nozzle outlet geometry. For example, when the guide pipe with the cavitator was installed, the aggressive intensity of the jet was four times larger than that without the guide pipe and the cavitator [39], and *St* became nearly a quarter of that of the jet without the guide pipe and the cavitator [75], as *fsedd* was decreased. Namely, the constant value of *St* should be unique to the nozzle outlet geometry. The investigation of *St* for various nozzles would be a future work, to enhance and/or control the aggressive intensity of the cavitating jet.

*σ σ σ* **Figure 11.** Aspects of periodical shedding of cloud cavitation: (**a**) *p*<sup>1</sup> = 20 MPa, σ = 0.014; (**b**) *p*<sup>1</sup> = 30 MPa, σ = 0.014; (**c**) *p*<sup>1</sup> = 20 MPa, σ = 0.02 [21].

#### **4. Key Parameters of Cavitating Jet**

#### *4.1. Type of Cavitating Jet*

As mentioned above, a cavitating jet normally means a submerged high-speed water jet with cavitation, i.e., a cavitating jet in water. Soyama developed the cavitating jet in air by injecting a high-speed water jet into a low-speed water jet without a water-filled chamber [33,35], for cavitation peening treatment outside of a tank and/or pipes [34,76]. When the residual stress of stainless steel was measured, it was reported that the cavitating jet in water introduced compressive residual stress in a deeper region, and the cavitating jet in air introduced large but shallow compressive residual stress on the surface [4]. The cavitation peening using the cavitating jet in water corresponds to shot peening using large shots, and that of the cavitating jet in air corresponds to shot peening using small shots at high velocity. Thus, the characteristics of the treated surface strongly depend on the type of the cavitating jet.

#### *4.2. Stando*ff *Distance*

As shown in Figure 9b, as cavitation is generated inside and/or outside of the nozzle, it becomes cloud cavitation and forms ring vortex cavitation on the impinging surface, then collapses, producing the impacts. In view of the practical applications of the cavitating jet, the working mechanism strongly depends on the standoff distance, which is the distance from the nozzle to the surface. Figure 12 illustrates the weight loss as a function of standoff distance at constant injection pressure, *p*<sup>1</sup> = 120 bar (12 MPa) [77]. In Figure 12, the weight loss means a kind of aggressive intensity of the jet. Two peaks are observed at each condition in Figure 12. For convenience, the peak near the nozzle side and the peak at the further nozzle side are named the 1st peak and 2nd peak, respectively. The 1st peak results from the impacts produced by water column collisions in the jet center. Even for the submerged water jet, the similar mechanism of water jet cutting is still active, whereas the affective region is limited near the nozzle. The 2nd peak is generated by the cavitation impact, as shown in Figures 6 and 9b. As cavitation is developed and then collapsed, a certain distance from the nozzle is required. As shown in Figure 12, the optimum standoff distance *sopt* of the 1st peak is scarcely affected by cavitation number, and that of the 2nd peak strongly depends on cavitation number. At *p<sup>2</sup>* = 2.4 bar (0.24 MPa), the weight losses at 1st peak and 2nd peak are 250 and 450 mg, respectively. They are 250 and 220 mg for *p<sup>2</sup>* = 3.0 bar (0.3 MPa) and 300 and 110 mg for *p<sup>2</sup>* = 3.6 bar (0.36 MPa). When the maximum values of the 1st peak and 2nd peak are compared, the value of the 2nd peak at *p<sup>2</sup>* = 0.24 MPa is 1.5 times larger than that of the 1st peak at *p<sup>2</sup>* = 0.36 MPa. Namely, at optimum cavitating conditions, the aggressive intensity due to cavitation impact, i.e., the 2nd peak, is larger than that of water column impact, i.e., 1st peak. Note that water jet peening and cavitation peening use the 1st peak and 2nd peak, respectively.

To avoid confusing cavitation peening and water jet peening, Soyama proposed a classification map for cavitation peening and water jet peening using standoff distance and cavitation number, as shown in Figure 13. Over 150 points were collected from references [39,74,77–89], and it was found that the line shown in Equation (8) distinguished between cavitation peening and water jet peening [90].

$$\frac{d\,^{S\_{opt}}}{d} = 1.8 \,\sigma^{-0.6} \tag{8}$$

One easy way to confirm the 2nd peak region, i.e., cavitating peening region, is as follows. Considering the aggressive intensity of the cavitating jet, the target materials are chosen. For example, in the case of the weak cavitating jet, a soft metal would be better. Then, the target is exposed to the jet. When a ring pattern is obtained, such as in Figure 8, it is a cavitation peening condition. The removal of paint can show the treatment area of cavitation peening [91]. A pressure-sensitive film also detects the treatment area of a cavitating jet [12].

**Figure 12.** Weight loss, i.e., the aggressive intensity of the jet, as a function of standoff distance changing with cavitation number at constant injection pressure. The peak at the near side of the nozzle, i.e., the 1st peak, was caused by water column impacts. The peak at the far side from the nozzle, i.e., the 2nd peak, was produced by cavitation impacts. The standoff distance of 2nd peak was changed by the cavitation number [77].

**Figure 13.** Classification map for cavitation peening and water jet peening considering standoff distance and cavitation number [90].

#### *4.3. Injection Pressure*

*β β β β* To show the effect of injection pressure on the processing capability of the cavitating jet, Figure 14 reveals the processing capability β at (a) the 1st peak, i.e., water jet peening, and (b) 2nd peak, i.e., cavitation peening, at the constant downstream pressure condition [92]. The processing capability β is defined by the arc height *h* of band steel made of the same material as Almen strip, considering the width of the steel *w<sup>s</sup>* and the peening width *wp*, as follows [92].

$$\alpha = \left| 1 - \frac{w\_p}{w\_s} \right| \tag{9}$$

$$\frac{1}{\rho} = \frac{8h}{L^2} \tag{10}$$

$$
\beta = \frac{1}{\rho} \left( 1 + a \right) \tag{11}
$$

Here, *L* is the length to measure *h*. In Figure 14, the processing capability, i.e., a kind of aggressive intensity, was obtained by the arc height of the peened plate, as the arc height using the Almen strip is commonly used to measure the peening intensity [93]. In the case of water jet peening, β increases with the injection pressure *p*1, as the peening effect is produced by water column impacts, which increases with *p*1. On the other hand, in the case of cavitation peening, β has a maximum at *p*<sup>1</sup> = 40 MPa at a constant downstream condition. When the maximum values of the 1st peak and 2nd peak are compared: β at the 2nd peak is 1.7 times larger than that at the 1st peak. As the jet power of 60 MPa is 1.8 times larger than that of 40 MPa, the peening efficiency of cavitation peening is about three times higher than that of water jet peening. Note that too high an injection pressure reduces the peening intensity of cavitation peening, as shown in Figure 14b. The reason the peening intensity of cavitation peening decreases at *p*<sup>1</sup> > 40 MPa is discussed in "5. Estimation of Aggressive Intensity of Cavitating Jet". *β β β*

α = ฬ1 െ

 1 = 8ℎ ଶ

௦ ฬ

**Figure 14.** Effect of injection pressure on the aggressive intensity of cavitating jet at constant downstream pressure condition; (**a**) 1st peak; (**b**) 2nd peak [92].

#### *4.4. Cavitation Number*

As the cavitating jet is the cavitating flow, the cavitation number is one of the key parameters of the cavitating jet. To reveal how important the cavitation number is in the cavitating jet characteristics, the relationship between the cavitation number and the optimum standoff distance is shown in Figure 15 [72]. In Figure 15a, the cavitating length is also added, and the length and the distance are normalized by effective nozzle diameter *D<sup>e</sup>* , which is defined by the nozzle throat diameter and discharge coefficient [72]. The data of the references are put in together in Figure 15 [77,81,86,87,89,94], as well as the relationship on each line, which is a straight line on a log–log scale. That is, the relationship can be described by Equation (12) [72].

$$\frac{s\_{opt}}{d} = c\_1 \sigma^{-c\_2} \tag{12}$$

 = ଵ Here, *c*<sup>1</sup> and *c*<sup>2</sup> are constants, and they depend on the nozzle geometry.

In Figure 15b, the expected standoff distance considering the pressure distribution and the local cavitation number on the impinging surface is also shown. The expected standoff distance is on the line of the log–log scale, and it is very close to the experimental result. The local cavitation number on the impinging surface is also an important parameter when considering the flow pattern of the impinging cavitating jet.

**Figure 15.** Relation between cavitation number and optimum standoff distance; (**a**) experimental result of optimum standoff distance and cavitating length; (**b**) estimated standoff distance from local cavitation number [72].

#### *4.5. Sound Velocity in Cavitating Flow Field*

In order to consider the key factors in cavitation impacts produced by the cavitating jet experimentally, sound velocity in the cavitating flow field has been investigated [95]. As erosion rate increases with acoustic impedance [96,97], the cavitation impact might increase with sound velocity in the cavitating flow field. In view of luminescence [62,98–100] and erosion [101], the effect of sound velocity of dissolved gas has been considered; however, there is no example of evaluating the sound velocity of the cavitating flow field itself. As is well known, the sound velocity changes drastically with the void ratio [102]. In the present review, the sound velocity in the cavitating flow field is considered.

To reveal the measuring method of the sound velocity of the cavitating flow field, Figure 16 shows the aspect of the cavitating flow through a Venturi tube. In Figure 16, the water flows from the left-hand side to the right-hand side. The cavitation occurs at the throat, and the vortex cavitation is observed at the end of the cavitating region. Further downstream of the vortex cavitation, many tiny bubbles are observed: they are residual bubbles [58] after cavitation collapses. As shown in Figure 16, the densities of residual bubbles differ greatly between the left and right sides of the white arrow. When this phenomenon is observed by a high-speed video camera as shown in Figure 17, the boundary of the density, which is marked by the yellow arrow in Figure 17, moves downstream. In the experiment, the upstream pressure in the absolute pressure of the throat was 0.6 MPa, and the ratio of the cross-sectional area between the throat and the tube was 9; thus, the flow speed at the downstream of the throat was about 3.5 m/s. On the other hand, the moving speed of the boundary marked by the yellow arrow was over 600 m/s. As shown in Figure 17, at the starting point of the boundary, the vortex cavitation shrank. This aspect suggests that the pressure wave is produced by vortex cavitation collapse, and the boundary movement reveals the pressure wave, as the residual bubbles are collapsed by the pressure wave. The sound velocity in the cavitating flow field can be estimated by measuring the moving speed of the pressure wave, i.e., the boundary of the density of the residual bubbles [95].

**Figure 16.** Aspects of vortex cavitation in the Venturi tube observed with a flash lamp. A pressure wave, which is indicated by the white arrow, was observed by collapsing residual bubbles.

**Figure 17.** Aspects of vortex cavitation in Venturi tube observed by a high-speed video camera. The recording speed of the camera was 51,999 fps. A pressure wave, which is indicated by a yellow, arrow was observed after the vortex cavitation collapsed.

To explore the effect of cavitation number on the sound velocity, Figure 18 shows the aspect of the pressure wave changing with the cavitation number. The pressure wave is indicated by the yellow arrow in Figure 18. In Figure 18, the velocity of the pressure wave, i.e., the sound velocity, is shown in the right-hand side of the aspect. Figure 19 illustrates the relation between cavitation number and the sound velocity. The sound velocity *v<sup>s</sup>* increases with cavitation number. At a relatively low void ratio, the sound velocity increases with a decrease in void ratio; thus, the tendency of the relation in Figure 18 is reasonable. As mentioned above, the erosion rate increased with the acoustic impedance [96,97], and the acoustic impedance is expressed as the product of the sound velocity, the density and the speed of the microjet; thus, the sound velocity in the cavitating flow field is a parameter of the cavitating jet. This result suggests that the increase in the sound velocity is one of reasons why the aggressive intensity of the cavitating jet increases with cavitation number, as the sound velocity increases with the cavitation number, as shown in Figure 19. On the other hand, the aggressive intensity of the

cavitating jet decreases with an increase in cavitation number, as the cavitating region and/or bubble size decreases with an increase in the cavitation number. These two conflicting tendencies are the aggressive intensity of the cavitating jet has a peak for the cavitation number. More details are provided in "5. Estimation of Aggressive Intensity of Cavitating Jet". Note that the aggressive intensity of the cavitating jet had a peak at σ = 0.01 − 0.014 [4,103,104], and that of the hydrodynamic cavitation through Venturi tube had a peak at σ = 0.4 − 0.7 [99]. These peaks might depend on the flow pattern of the cavitating flow and/or the density of the residual bubbles. *σ* − *σ* −

**Figure 18.** Aspects of pressure wave in Venturi tube changing with cavitation number. The recording speed of the camera was 109,999 fps.

**Figure 19.** Sound velocity in cavitating flow field through Venturi tube changing with cavitation number.

#### *4.6. Nozzle Geometry and Diameter*

As mentioned above, the vortex cavitation in the cavitating jet is an important phenomenon, and the nozzle geometry, especially nozzle outlet geometry, is one of the key factors of the aggressive intensity of the cavitating jet, as the vortex is strongly affected by the nozzle outlet geometry. Figure 20 shows a schematic diagram of nozzle outlet geometry and the relative aggressive intensity of the cavitating jet from data published in previous reports [28,39,75,79]. Whereas a similar figure was introduced by Soyama [3], the data of nozzles K and L [28,75] were added in Figure 20. Nozzles A–E in Figure 20 are conventional nozzles for a water jet, and nozzle F is the standard nozzle for a standard test method for the erosion of solid materials by a cavitating jet [15]. Nozzle G is the nozzle obtained by optimizing the outlet bore and length experimentally [74]. Nozzle I and J had a guide pipe and a cavitator, as these enhance the aggressive intensity by about two times [39]. When both the guide pipe and the cavitator were installed, the aggressive intensity became four times larger that of nozzle J [39]. For nozzle K, when water flow holes were made near nozzle outlet, the aggressive intensity improved by 34%. When a long guide pipe with holes and water flow holes near the nozzle outlet were installed for nozzle L [28], the aggressive intensity of its cavitating jet L was 2.5 times larger than that of nozzle J and nearly 60 times larger than that of conventional water jet nozzles, as shown in Figure 20. The effect of nozzle geometry was also discussed in [14]. Resonating nozzles were also proposed and investigated [105,106].

**Figure 20.** Effect of nozzle outlet geometry on aggressive intensity of cavitating jet.

The other important factor regarding nozzles is nozzle size. As mentioned above, vortex cavitation is an important phenomenon in the cavitating jet. In view of the Reynolds number, which is a key parameter of vortical flow, the larger velocity and the larger size are effective for the cavitating jet. However, too high a speed, i.e., too high an injection pressure, decreases the aggressive intensity of the cavitating jet, as shown in Figure 14b, as the sound velocity is decreased at too low a cavitation number, i.e., too-high speed condition. Thus, in the case of practical applications of the cavitating jet, the cavitating jet using a large nozzle at a relatively low injection pressure is better than that of a small nozzle at a high injection pressure [80]. The scaling law of nozzle size on the aggressive intensity of the cavitating jet is discussed in reference [107].

*σ*

#### *4.7. Water Qualities*

As cavitation is a phase-change phenomenon from liquid phase to gas phase, as mentioned above, the water temperature affects the aggressive intensity of the cavitating jet [108–110]. It was reported that peening intensity using the cavitating jet at 288–308 K was nearly constant [4]. Whereas cavitation nuclei are required for the cavitating jet, too many air bubbles reduce the aggressive intensity of the cavitating jet due to the cushion effect [2] and the decrease in the sound velocity. When the water-filled chamber for the cavitating jet is too small, the suction vortex caused by the jet reduces the aggressive intensity of the cavitating jet. Thus, in the report, the effects of water depth and the chamber size on the aggressive intensity were investigated experimentally [111]. In the report [111], the effect of gas content on the peening intensity using the cavitating jet was also investigated using degassed water. ඥ<sup>ଵ</sup> െ <sup>ଶ</sup> ௩ ∝ ሺ௩ሻ ଷ ∙ ሺ<sup>ଶ</sup> െ <sup>௩</sup> ሻ ∙ ሺ<sup>௦</sup> െ ௦ ௧ሻ ∙ ሺඥ<sup>ଵ</sup> െ ଶሻ <sup>ଵ</sup> െ <sup>ଶ</sup> ඥ<sup>ଵ</sup> െ <sup>ଶ</sup> ඥ<sup>ଵ</sup> െ <sup>ଶ</sup>

ଷ

∙ ሺ<sup>ଶ</sup> െ ௩ሻ

௩ ∝ ሺ௩ሻ

#### **5. Estimation of Aggressive Intensity of Cavitating Jet** ඥ<sup>ଵ</sup> െ <sup>ଶ</sup>

As mentioned above, in the constant downstream pressure condition, too high an injection pressure, i.e., too low a cavitation number, reduces the aggressive intensity of the cavitating jet. In this section, the mechanism is discussed, and a method to estimate the aggressive intensity of the cavitation jet as a function of cavitation number is proposed. ௩ = <sup>ଷ</sup> ିଵ.଼ ∙ ሺ<sup>ଶ</sup> െ <sup>௩</sup> ሻ ∙ ሺ<sup>௦</sup> െ ௦ ௧ሻ ∙ ሺඥ<sup>ଵ</sup> െ ଶሻ ଷ 

In Figure 21, the experimental results of the aggressive intensity of the cavitating jet as a function of cavitation is revealed by blue closed circles [112]. The aggressive intensity is normalized by the maximum value, and it has a peak at σ = 0.016. As the energy of cavitation is proportional to the volume of the cavitation and the pressure difference of the bubble [111], the aggressive intensity of the cavitating jet *Icav* can be assumed as follows. <sup>௦</sup> = <sup>ସ</sup> <sup>ହ</sup>

$$I\_{\rm cav} \propto \left(L\_{\rm cav}\right)^3 \cdot \left(p\_2 - p\_v\right) \tag{13}$$

**Figure 21.** Aggressive intensity of cavitating jet changing with cavitation number.

As *Icav* is affected by the acoustic impedance [95,111], the term of the sound velocity *v*<sup>s</sup> is included in Equation (13). As *Icav* is also affected by flow velocity, which is defined by the pressure difference, i.e., √ *p*<sup>1</sup> − *p*2, the velocity term is also included in Equation (13).

$$I\_{\rm cav} \propto (L\_{\rm env})^3 \cdot (p\_2 - p\_v) \cdot (v\_s - v\_{s\
th}) \cdot \left(\sqrt{p\_1 - p\_2}\right)^n \tag{14}$$

Here, *vs th* is the threshold level of the sound velocity considering the threshold level of *Icav* [43]. The *n* is put as an exponent in Equation (14) to consider the power law of the velocity on *Icav* [76,105,111–113]. If it was assumed that *Icav* was proportional to the flow energy, which was defined by the product of the pressure and the flow rate, the flow energy is proportional to *<sup>p</sup>*<sup>1</sup> <sup>−</sup> *<sup>p</sup>*<sup>2</sup> and <sup>√</sup> *<sup>p</sup>*<sup>1</sup> <sup>−</sup> *<sup>p</sup>*2, as <sup>√</sup> *p*<sup>1</sup> − *p*<sup>2</sup> was proportional to the flow velocity. Namely, the flow energy is proportional to cube of √ *p*<sup>1</sup> − *p*2. Thus, in the present review, *n* = 3 is chosen. From Equations (8) and (14), Equation (15) was obtained.

$$I\_{\rm cav} = c\_3 \, \sigma^{-1.8} \cdot (p\_2 - p\_v) \cdot (v\_s - v\_{s\
th}) \cdot \left(\sqrt{p\_1 - p\_2}\right)^3 \tag{15}$$

Here, *c*<sup>3</sup> is constant. When Equation (15) is assumed, *Icav* can be estimated by obtaining *c*3, *c*4, *c*<sup>5</sup> and *vs th* by a least-square method.

$$
v\_s = c\_4 \,\,\sigma + c\_5 \,\,\,\tag{16}$$

The estimated *Icav* is shown in Figure 21 by empty squares. The correlation coefficient between experimental data and estimated values is 0.924. As the number of the datasets was 11, the probability of non-correlation is less than 0.005%. If the probability of a non-correlation is less than 1%, it can be concluded that the relationship is highly significant. Thus, it can be concluded that the relationship between the experimental data and the estimated values is highly significant. *Icav* can be estimated by Equation (15).

#### **6. Applications of Cavitating Jet**

The cavitating jet has been applied for drilling and cutting rocks [16,105,106,113]. The cavitating jet in air can also dig concrete structures for the maintenance of infrastructure [114]. In the case of mechanical surface treatment, the use of a submerged water jet was attempted to mitigate stress corrosion cracking (SCC) of nuclear power plants by using impinging impacts of a water column in the jet center at the beginning [84]. Soyama et al. found that the cavitating jet could introduce compressive residual stress into stainless steel by using cavitation impacts [13]; then, the cavitating jet was successfully applied to mitigate SCC in nuclear power plants [115]. Based on experimental results of the introduction of compressive residual stress into metallic materials, the improvement in the fatigue strength of metallic components by cavitation peening was proposed using a pressurized chamber to enhance the aggressive intensity of the cavitating jet [116–118], and it was demonstrated for forging die [119], gears [120,121], continuous valuable transmission CVT elements [122] and rollers [123]. After enhancing the aggressive intensity of the cavitating jet by optimizing the nozzle geometry, the improvement in fatigue strength by cavitation peening with an open chamber was demonstrated [17,124–127]. Cavitation peening using oil was also proposed [128]. Cavitation peening also improves tribological properties such as fretting fatigue properties [129,130]. The improvement in fatigue strength of metallic materials using a cavitating jet in air was also demonstrated [36,131,132]. In view of environmental-assisted cracking, the suppression of hydrogen-assisted fatigue crack growth in austenitic stainless steel and delayed fracture resistance on chrome molybdenum steel by cavitation peening were reported [133,134].

The cavitating jet can be applied in the semiconductor industry. It can be used not only for cleaning [135], but also gettering [136]. When erosion rates using ultrasonic vibratory apparatus [137] and cavitating jet apparatus [15] were compared, the erosion rate of the cavitating jet apparatus was larger than that of ultrasonic vibratory apparatus. As conventional ultrasonic cleaning devices use ultrasonic cavitation, cleaning using the cavitating jet is more powerful. The gettering technique is a method to remove unwanted impurities from active device regions in a silicon wafer into the back side of the wafer by introducing oxidation-induced stacking faults (OSF) [138–142]. In order to produce OSF, the introduction of strain into the wafer is required, and shot blasting is used in a conventional way [143,144]. However, the cleaning of shots is required. In the case of gettering using the cavitating jet, the cavitating jet can introduce strain for OSF [145,146] and cleaning at the same time. This is a great advantage for semiconductor processes.

One of the applications that was proposed in the bioengineering area is oral cleaning using the cavitating jet [147–149]. Dental implants have been used as the solution for the loss of teeth; however, peri-implant mucosis and peri-implantitis are new dental diseases affecting implants [150,151]. The most effective treatment for these diseases is cleaning dental plaque, which is a kind of biofilm that adheres to the surface of teeth and implants [152]. Conventional cleaning methods of the dental plaque are oral brushing, ultrasonic scaling and rubber cup cleaning. Unfortunately, micro-textured roughness is made on the implant surface to improve biocompatibility [153], so oral toothbrushes and the tips of scaling instruments cannot reach the bottom of the rough surface. The cleaning of dental plaque on the rough surface of dental implants using the cavitating jet has been successfully demonstrated [147–149].

As is well known, sonochemistry is a research area of chemistry using ultrasonic cavitation [154]. Hydrodynamic cavitation such as the cavitating jet and cavitating flow through orifices can be applied for wastewater treatment [155–162] and to oxidize organic compounds [163]. The dispersion of spilled oil by a cavitating jet at sea has also been proposed [164]. It was reported that the efficacy of hydrodynamic cavitation was 20 times better than that of ultrasonic cavitation when the efficiency of the hydrodynamic cavitation on the pretreatment of biomass was compared with that of ultrasonic cavitation [165].

#### **7. Conclusions**

To use the cavitating jet for practical applications, the unsteady behavior of the cavitating jet, i.e., a submerged water jet with cavitation, was reviewed. The key factors on the aggressive intensity of the cavitating jet were also summarized. In the present review, the aggressive intensity of the cavitating jet was investigated by erosion rate and/or peening intensity. The main topics reviewed in the paper are summarized as follows.


**Funding:** This research was partly supported by JSPS KAKENHI grant numbers 18KK0103 and 20H02021. **Conflicts of Interest:** The author declares no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Insight into the Kelvin-Helmholtz Instability Appearance in Cavitating Flow**

**Peter Pipp , Marko Hoˇcevar and Matevž Dular \***

Faculty of Mechanical Engineering, University of Ljubljana, Askerceva 6, 1000 Ljubljana, Slovenia; peter.pipp@fs.uni-lj.si (P.P.); marko.hocevar@fs.uni-lj.si (M.H.)

**\*** Correspondence: matevz.dular@fs.uni-lj.si; Tel.: +386-1-477-1453

**Abstract:** Recently the development of Kelvin-Helmholtz instability in cavitating flow in Venturi microchannels was discovered. Its importance is not negligible, as it destabilizes the shear layer and promotes instabilities and turbulent eddies formation in the vapor region, having low density and momentum. In the present paper, we give a very brief summary of the experimental findings and in the following, we use a computational fluid dynamics (CFD) study to peek deeper into the onset of the Kelvin-Helmholtz instability and its effect on the dynamics of the cavitation cloud shedding. Finally, it is shown that Kelvin-Helmholtz instability is beside the re-entrant jet and the condensation shock wave the third mechanism of cavitation cloud shedding in Venturi microchannels. The shedding process is quasi-periodic.

**Keywords:** cavitation; Kelvin-Helmholtz instability; microchannel; numerical simulation

**Citation:** Pipp, P.; Hoˇcevar, M.; Dular, M. Numerical Insight into the Kelvin-Helmholtz Instability Appearance in Cavitating Flow. *Appl. Sci.* **2021**, *11*, 2644. https:// doi.org/10.3390/app11062644

Academic Editor: Florent Ravelet

Received: 12 February 2021 Accepted: 12 March 2021 Published: 16 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

When the pressure drops below the vapor pressure, cavitation, that is, the formation of vapor cavities within a homogeneous liquid medium occurs. The vapor structures are precarious, and they frequently collapse vigorously as they encounter an area of elevated pressure [1]. The mechanism of cavitation cloud shedding is among the most distinctive features of evolved cavitation. Two unique mechanisms are well established to influence the process of cloud shedding [2–7].


We recently published on the visualization experiment of cavitation within a microchannel [8]. Although the initial objective of the analysis was to create supercavitating conditions within it, we discovered that due to the development of a Kelvin-Helmholtz instability, which induces a semi periodic collapse of the attached cavity, this regime is suppressed. Detailed experimental studies using high-speed photography with visible light and X-rays, showed that this is indeed a third process contributing to the shedding of cloud cavitation in addition to the re-entrant jet and the shock wave mechanisms.

The purpose of this paper is to provide greater insight into the occurrence of Kelvin-Helmholtz instability of cavitation structures in the Venturi microchannel. Above all, we will show how and why the phenomenon occurs, and how it affects the dynamics of cavitation. In doing so, we will focus primarily on showing the results of vapor volume fractions, velocity profiles, and velocity field, thus showing an insight into what is happening not only around and at the boundary, but also within cavitation structures.

We give a very brief summary of the experimental findings in the following sections and then use a computational fluid dynamics (CFD) study to peek deeper into the onset of the Kelvin-Helmholtz instability and its effect on the dynamics of the cavitation cloud shedding.

#### **2. Experimental Observations**

The experimental method is seen here only briefly, a more detailed study can be found in [8]. The set-up is illustrated in Figure 1. By sandwiching them between two acrylic glass plates, microchannels are constructed of 450 µm thick stainless-steel sheets that have a laser-cut convergent-divergent narrowing. The Venturi channel's convergent angle measures 18◦ and the divergent angle is 10◦ , with a throat height of 675 µm. Both the channel entrance and exit are perpendicular to the cross-section plane, with the channel inlet on the left and the outlet on the right side of the channel.

**Figure 1.** Microchannel geometry with the 675 µm throat gap and the 450 µm channel width.

High-speed cameras (Photron SA-Z and Photron Mini AX200, Photron, Tokio, Japan) have captured cavitation images at a framerate of typically 200,000 fps in both the visible and X-ray light spectrum. Typical flow features that imply the formation of the Kelvin-Helmholtz instability in the channel are shown in Figure 2. The flow is from the left side to the right.

**Figure 2.** Development of the Kelvin-Helmholtz instability in the microchannel ( · m = 9.15 g/s, ∆p = 4.00 bar, σ = 1.24).

Figure 2 illustrates common microchannel representations of the Kelvin-Helmholtz instability. In particular, these images are from an experiment with a mass flow ( · m) 9.15 g/s,

a pressure difference between inlet and outlet (∆p) 4 bar, and a cavitation number (σ) of 1.24. A closer look at the flow conditions shows that the attached cavity first forms in the shape of a large single cavitation bubble extending downstream from the starting point until the pressure rises well above the saturation pressure—the presence of cavitation is like the supercavitation state. Kelvin-Helmholtz instability forms when the bubble approaches its maximum size. The ripple of the interface increases and the instability encircles the entire cavitation cloud while rolling up and induces its very rapid dissolution in the flow. The pattern then continues periodically.

The numerical simulation that provides further information in the creation of the Kelvin-Helmholtz instability and the cavitation dynamics in the microchannel is presented in the next section.

#### **3. Numerical Procedure**

A relatively straightforward simulation was carried out using the commercial Fluent 2020 R2 software package (Ansys Inc, Canonsburg, PA, USA). The solution was obtained with the Unsteady Reynolds-averaged Navier-Stokes (URANS) equations, with Reboud correction modified SST k-ω turbulence model and the Schnerr-Sauer cavitation model. The next section provides some more information.

#### *3.1. Mesh*

We regarded the matter two-dimensionally to promote the computation of the case and to expedite the simulation calculation. Thus, as seen in Figure 3, the channel inlet and outlet are modified appropriately. The structured computational mesh was used. The mesh independence was tested on 3 meshes while observing the average pressure difference ∆p and the average cavity length l. The discretization error of less than 0.6% was calculated by Richardson extrapolation [9] for the medium mesh with around 160,000 cells at a flow rate of · m = 9.03 g/s (Table 1).

**Figure 3.** The geometry of the microchannel computational domain and the mesh detail in the wedge apex area.



Along with the channel height, the final mesh had a constant number of 67 cells. Therefore, 10 µm is the height of an individual cell at the wedge apex. In relation to the geometry, the cells grow in the length to height ratio towards the inlet and the outlet of

the domain. We have considered the boundary layer, where the overall thickness of the boundary layer of 50 µm was defined, with 10 layers and a growth rate of 1.2, which is a sufficiently small mesh to maintain the y+ values below 5.

#### *3.2. Reynolds-Averaged Navier-Stokes (RANS) Equations*

The Navier-Stokes equations are a set of partial differential equations, describing the conservation of momentum and mass for a viscous fluid flow motion. A direct numerical simulation, in which the Navier-Stokes equations are solved without a turbulence model, requires a huge amount of computing power because the whole range of spatial and temporal scales of the turbulence must be resolved. As this is not practical, averaging of variables is used. The idea behind averaging is the Reynolds decomposition, an instantaneous quantity decomposition into its time-averaged Φ and fluctuating quantities ϕ'(t) with zero mean value. Reynolds-averaged Navier-Stokes equations will be used in the form, that is generally accepted in the commercially available computational fluid dynamics (CFD) solvers. Here, the overline indicates a time-averaged variable, while the tilde indicates a density-weighted averaged or Favre-averaged variable [10]. The Favre-averaging method was selected because in the present case turbulent fluctuations lead to sizeable density fluctuations. The continuity equation is as follows, where *ρ* is the fluid density, *t* time, and *U* velocity vector, defined as *U* = [*U*, *V*, *W*] T :

$$\frac{\partial \overline{\rho}}{\partial t} + \nabla \cdot \left(\overline{\rho}\tilde{\mathbf{U}}\right) = 0 \tag{1}$$

The momentum equations in the *x* and *y* directions are as follows, where *P* is pressure, *µ* dynamic viscosity, *u* ′ , *v* ′ and *w* ′ velocity fluctuations of velocity vector *U*, while *SMx* and *SMy* represent momentum source terms in *x* and *y* directions, respectively

$$\begin{aligned} \frac{\partial(\overline{\rho}\overline{\mathcal{U}})}{\partial t} + \nabla \cdot (\overline{\rho}\overline{\mathcal{U}}\overline{\mathcal{U}}) &= -\frac{\partial \overline{\mathcal{V}}}{\partial x} + \nabla \cdot (\mu \nabla \overline{\mathcal{U}}) + \left[ -\frac{\partial \overline{(\overline{\rho}u'^2)}}{\partial x} - \frac{\partial (\overline{\overline{\rho}u'v'})}{\partial y} - \frac{\partial (\overline{\overline{\rho}u'w'})}{\partial z} \right] + S\_{Mx} \\ \frac{\partial(\overline{\rho}\overline{\mathcal{V}})}{\partial t} + \nabla \cdot (\overline{\rho}\overline{\mathcal{V}}\overline{\mathcal{U}}) &= -\frac{\partial \overline{\mathcal{V}}}{\partial y} + \nabla \cdot (\mu \nabla \overline{\mathcal{V}}) + \left[ -\frac{\partial (\overline{\overline{\rho}u'v'})}{\partial x} - \frac{\partial (\overline{\overline{\rho}v'^2})}{\partial y} - \frac{\partial (\overline{\overline{\rho}v'w'})}{\partial z} \right] + S\_{Mx} \end{aligned} \tag{2}$$

For the accurate representation of cavitating flow details, some authors perform compressible flow simulations [11–13]. In the compressible case, the energy conservation equation is added to the above set of continuity and momentum equations. Also, equations of state of vapor and liquid are introduced. On the other hand, some authors have shown, that sufficient flow representation accuracy was obtained in the incompressible flow case, much decreasing the computational load requirements (for example [2,14]).

#### *3.3. Turbulence Modeling*

Increased available computational power has made possible advances in computational fluid dynamics, among them Large-Eddy Simulation (LES) and Detached-Eddy Simulation (DES). These relatively novel approaches are in research used to accurately resolve all scales of turbulent fluctuations, while for engineering purposes there is usually no need for such complexity. In above mentioned Reynolds-averaging models, turbulent fluctuations are not resolved, and the turbulence model is used instead. The k-ω and k-ε turbulence models are the most common models used to simulate mean flow characteristics for turbulent flow conditions. Both are two-equation models that provide a general description of turbulence utilizing two transport equations. The Menter's Shear Stress Transport (SST) k-ω turbulent model [15], combines both k-ω and k-ε turbulence models such that it merges the robust, accurate, and less y+ sensitive performance of the k-ω model in the near-wall region with the solid freestream operation of the k-ε model in the far-field. Nowadays the SST k-ω turbulent model is the preferred two-equation turbulent model and the preferred model for a wide range of applications.

For the modeling of nonstationary cavitation flows, among them, periodic cavitation cloud shedding, two-equation turbulent models in the original form are not the preferred option. Two equation models prevent the reentrant jet from appearing at the rear of the cavitation cloud, all this being caused by too high prediction of turbulent viscosity. For a better non-stationary turbulence modeling performance, Reboud et al. [16] proposed a modification of the turbulent model. In the modification, the turbulent viscosity of the mixture is reduced in regions with high void fractions and the equation of turbulent viscosity becomes as follows, with *k* representing turbulence kinetic energy and *ω* specific turbulence dissipation rate:

$$
\mu\_t = f(\rho) \frac{k}{\omega} \tag{3}
$$

The density function is set by Equation (4), here indices *l*, *v, m* represent liquid, vapor, and mixture.

$$f(\rho) = \rho\_{\upsilon} + \frac{\left(\rho\_{\text{m}} - \rho\_{\upsilon}\right)^{n}}{\left(\rho\_{l} - \rho\_{\upsilon}\right)^{n-1}} \text{  $n \gg 1$ } \tag{4}$$

Various values were proposed for the exponent n. Coutier-Delgosha et al. [17] proposed value *n* = 10, while the described turbulent model extension was used and validated for various applications, among them in Venturi channel and a hydrofoil [2,3,18–20].

#### *3.4. Two-Phase Flow Modeling*

We more commonly use the concept of the homogeneous flow of the mixture for cavitation simulation, where the two-phase flow is assumed to be a single-phase flow of the liquid-vapor mixture. This helps one to solve only one motion equation since we consider the problem as a single-phase, but with the mixture's properties. Accordingly, the properties of the liquid-vapor mixture are described by the proportion of the vapor phase, using the model suggested by Bankoff [21]. The mixture's density is written:

$$
\rho\_{\mathcal{W}} = \mathfrak{a}\rho\_{\mathcal{U}} + (1 - \mathfrak{a})\rho\_{l} \tag{5}
$$

where *α* denotes vapor volume fraction, and mixture's dynamic viscosity is written as:

$$
\mu\_m = \alpha \mu\_\upsilon + (1 - \alpha)\mu\_l \tag{6}
$$

The equations of mass and momentum conservation are resolved in the model of the homogeneous flow of the mixture by the properties of the mixture, and the equation of phase fraction conservation must be solved:

$$\frac{\partial}{\partial t}(\rho\_{\upsilon}\alpha) + \nabla \cdot (\rho\_{\upsilon}\alpha \mathbf{u}\_{\upsilon}) = R\_{\varepsilon} - R\_{\varepsilon} \tag{7}$$

where *R<sup>e</sup>* and *R<sup>c</sup>* represent mass transfer source terms for evaporation and condensation, which account for the mass transfer between the phases of liquid and vapor in cavitation and are thus related to the growth and collapse of the vapor bubbles. According to the cavitation model used, their formulation varies.

#### *3.5. Cavitation Model*

We know multiple models of cavitation in which we can model liquid and vapor evaporation and condensation. Bubble dynamics models, introduced by Kubota et al. [22], where he used a linear portion of the Rayleigh-Plesset equation to explain the evolution of bubble radius as a function of surrounding pressure have been the most used cavitation models in CFD in recent years. The proportion of the vapor phase and hence the density of the mixture is calculated by bubble radius and bubble number density. Other cavitation models were derived relying on the pressure and bubble radius dependencies [23], based on the Rayleigh-Plesset equation. However, with all cavitation models, it is conditioned to specify parameters, such as the bubble number density or the initial size of the bubble, which are very hard to define. The Schnerr-Sauer [24] model is the cavitation model based on bubble dynamics that is most widely used in commercial CFD software packages, where the mass transfer source terms *R<sup>e</sup>* and *R<sup>c</sup>* (see Equation (7) are defined as:

when local pressure is equal to or smaller than vapor pressure—*p* ≤ *pv*:

$$R\_e = F\_{evap} \frac{\rho\_v \rho\_l}{\rho} a (1 - a) \frac{3}{R\_b} \sqrt{\frac{2}{3} \frac{(p\_v - p)}{\rho\_l}} \tag{8}$$

when local pressure is larger than vapor pressure—*p* > *pv*:

$$R\_c = F\_{cond} \frac{\rho\_v \rho\_l}{\rho} \mathfrak{a} (1 - \mathfrak{a}) \frac{3}{R\_b} \sqrt{\frac{2}{3} \frac{(p - p\_v)}{\rho\_l}} \tag{9}$$

where the default values of empirical calibration coefficients of evaporation *Fevap* and condensation *Fcond* are 1 and 0.2, respectively. To link the fraction of the vapor volume to the number of bubbles per liquid volume *n<sup>b</sup>* , the Schnerr-Sauer cavitation model uses:

$$\alpha = \frac{n\_b \frac{4}{3} \pi R\_b^3}{1 + n\_b \frac{4}{3} \pi R\_b^3} \tag{10}$$

The parameters that must be determined in this model are bubble radius *R<sup>b</sup>* and bubble number density *n*. The compressibility of water and water vapor was not considered in the modeling.

#### *3.6. Boundary Conditions*

Modeling cavitating flow in the microchannel was performed at several mass flow rates. The boundary condition was set at the inlet of the Venturi channel computational domain as average velocity. The inlet turbulence intensity was set to zero. The vapor volume fraction at the inlet was set to zero also.

At the outlet of the Venturi channel computational domain the second boundary condition was set prescribing the absolute pressure 1 bar and backflow turbulence intensity 5% for all mass flow rates. The vapor volume fraction at the outlet was set to zero.

On all walls, the no-slip stationary wall boundary condition was used with a standard wall roughness model. The temperature was in all cases set to 20 ◦C.

#### *3.7. Physics and Solver Settings*

It is important to set proper simulation settings, as are the initial conditions from which to start the simulation, so we first performed simulations as single-phase simulations for each mass flow under stationary conditions. The outcomes were then considered as the initial conditions for further transient simulations. Numerical calculations were carried out using time-dependent Reynolds-averaged Navier-Stokes equations. Homogeneous water and water vapor mixture was considered and a Schnerr-Sauer cavitation model [24] with an evaporation pressure of 2340 Pa was used. A modified SST k-ω model, using the turbulent viscosity correction suggested by Reboud et al. [16] and Coutier-Delgosha et al. [17], mentioned above, was used for the turbulent model. For pressure and velocity coupling, the PISO algorithm [25] was used. We used the second-order upwind scheme [26] for spatial numerical discretization of solving hyperbolic partial differential equations for all but the pressure and volume fraction, providing more precise results with marginally higher computer resource usage compared to the first-order upwind scheme. The PRESTO interpolation scheme was used to discretize the pressure [27], and the first-order upwind scheme [26] was used for the discretization of the volume fraction. Time integration was modeled for transient solutions using the implicit transient formulation of bounded second-order.

The convergence criteria were established by evaluating the progression of various flow parameters (absolute pressure at the inlet and outlet, and velocities at the microchan-

nel inlet, outlet, and throat). Following the sum of the imbalance of transport equations between iterations for all cells in the computational domain, the observed flow parameters were converged when falling below 10−<sup>5</sup> of the iterative numerical solution of the individual equations at each time step of the simulation. The iteration error was estimated to be less than 0.02%. By assessing its effect against the average pressure difference and the cavity length, the size of the time step was obtained. There was little difference in these parameters if the time step was shorter than 5 µs, but a shorter one—1 µs was ultimately chosen for observation of the Kelvin-Helmholtz instability. We conducted 50 ms of computational modeling for each case, where the last 30 ms was applicable for further analysis.

#### **4. Results**

First, we assess the capability of the numerical simulation to capture the general cavitation dynamics which was observed inside the microchannel—more specifically the onset of the Kelvin-Helmholtz instability (Figure 4). The flow is from the left to the right, the time difference between the images is 100 µs.

**Figure 4.** Separation of cavitation clouds and Kelvin-Helmholtz instability formation (experiment: · m = 9.15 g/s, ∆p = 4.00 bar, σ = 1.24; simulation: · m = 9.03 g/s, ∆p = 3.71 bar, σ = 1.26). The time difference between the images is ∆t = 100 µs.

The experiment and simulation both reveal the same story, but in the case of simulation, the Kelvin-Helmholtz instability is more pronounced and can be understood better and described more in detail.

In Figure 5, the development of the Kelvin-Helmholtz instability is shown with velocity profiles drawn over at cross-sections at horizontal intervals of 3 mm in the downstream direction, starting at the wedge apex.

Fluid flow velocity at the throat of the Venturi microchannel is 30 m/s, where liquid flows over the cavity and descents afterward towards the bottom wall of the channel, where the flow divides into the upstream and downstream flow (1). Vapor fraction α profiles for both phases are approximately uniform in both the fluid and vapor layers. They are however discontinuous at the straight interface between the two (step 1). Shear stress acts on the fluids due to discontinuity of the velocity gradient (velocity gradients are shown by the black line in Figure 5), while vapor gradient locally promotes the generation of instabilities and vorticity. The backflow due to the presence of the adverse pressure gradient at 6 mm downstream from the throat in the vapor phase increases from a peak value of 5 m/s to 13 m/s (step 2), which further increases shear stress. Increased shear stress between the liquid and cavity interface results in the formation of ripples on the interface, while the reentrant jet progresses further upstream. The attached cavity expands and reaches its full size (step 2) and shortly thereafter (100 µs) in step (3) one can see the first beginnings of ripples on the upper edge area of the vapor phase (approx. 6 mm from the Venturi throat). This promoting the formation of instabilities (step 3) in addition to the already present vertical shear stress between the liquid flow and the cavity interface enables the formation of Kelvin Helmholtz instability.

**Figure 5.** Development of Kelvin-Helmholtz instability with velocity profiles at individual crosssections. The time difference between the images is ∆t = 100 µs.

In step (4), 300 µs after the start of the instability formation process, we already observe the formation of unstable local turbulent structures approx. 6 mm from the Venturi throat, which are now confined to the shear layer. Due to their limited size, velocity gradients in step (4) are still unchanged in comparison to step (1) and the rolling motion of turbulent structures is not yet recognized. With time, the surface becomes more and more deformed from the original straight boundary interface, and turbulent structures grow in size and number (step 5). It seems that the vapor region is affected by the formation of turbulent structures, here velocity, density, and momentum are low. Velocity profiles in the vapor region change with instabilities continuously protruding deeper in the vapor region, finally leading to the situation in step (5).

Instabilities in the vapor phase continue to grow from step (5), and the region of high pressure between the cavities moves in the opposite direction to the main flow towards the Venturi throat, from approx. 4 mm in step (5) to approx. 2 mm in step (6). Advance in the direction opposite is related to the presence of the adverse pressure gradient (step 6), with instabilities and turbulent eddies in later steps reaching the Venturi throat section and causing full cavity separation and partial destruction of the vapor region in the vicinity of the throat. The vapor phase in the vicinity of the throat does not reappear until step (10) 900 µs after the start of the process. In steps (7–9) we notice the change of the velocity profiles, caused by the separation of the cavity from the Venturi throat and hence smoothing of the velocity profiles in the horizontal direction. The region of constant velocity near the upper Venturi section wall is reduced and the velocity profiles in horizontal direction become skewed. Regions of backflow are reduced in size and intensity. In steps (7–9) generation of Kelvin Helmholtz instabilities increase further downstream of the detached cavitation cavity. Finally, the entire detached cavitation cloud becomes

engulfed in the border instabilities. In addition to the change of velocity gradient, the vapor fraction α in the vertical direction becomes non-homogeneous and areas of slightly higher vapor fraction α can be traced to the region where instabilities were initially formed. The detached cavitation cloud together with all instabilities continues to flow downstream and finally dissipates (not shown in this sequence).

In step (10), 900 µs after the start of the instabilities formation, cavitation cavity is again formed at the Venturi throat section and the entire process is about to continue in a quasi-periodic way. More insight into the instability formation process may be gained from Figure 6, which shows an enlarged segment of the Venturi section relative to Figure 5.

**Figure 6.** A close-up view of the initial development of Kelvin-Helmholtz instability. The time difference between the images is ∆t = 20 µs. The black line shows a 10% vapor fraction region limit.

Velocity vectors reveal the presence of the turbulent eddies within the vapor phase and their interaction with the main flow above the discontinuity of the horizontal velocity. Figure 6 also shows contours of 10% vapor region location.

Step (1) shows the initial situation with a turbulent eddy on the right side of the observation region. Just next to the lower limit of 10% vapor region is the turbulent eddy impingement point on the shear layer. On the left of the 10% vapor fraction region, we observe the backflow, while on the right the turbulent eddy rotation causes the flow velocity to be in the direction of the main flow. Impingement point marks the position of the large variation of the shear stress gradient, on the left of the impingement point, the stress is higher than on the right. The horizontal 10% vapor fraction region limit of the velocity discontinuity is just above the impingement point. The transition from step (1) to step (2) is marked with an additional depression of the horizontal 10% vapor fraction region limit. The location of the impingement point and the horizontal 10% vapor fraction region has moved to the right in step (2) relative to step (1), while the lower 10% vapor fraction region has widened from the original location in step (1) to the newly established impingement point of step (2).

As the instability process continues in step (3) we observe gradual folding of the lower 10% vapor fraction limit, even before we notice significant variations in the velocity within the shear layer. We, therefore, speculate that the vapor fraction properties influence the instability process through the variation of the flow density and the transfer of the momentum. In step (3) we notice velocity variations within the main flow, corresponding well to the vapor fraction variations in intensity and location. The process, first observed in step (3) continues with step (4) showing increased spatial variations of the vapor fraction, followed by velocity variations in the main flow. As a consequence, turbulent eddy, first responsible for impinging the shear layer, is decayed into two turbulent eddies with an additional one formed near the location of the original impingement point. The location of the lower 10% vapor fraction region limit has moved in the direction of the Venturi section throat.

Steps (5) to (8) show markedly increasing folding of the shear layer, coupled with the additional breaking of the turbulent eddy into several smaller ones, including interactions among neighboring turbulent eddies. Gradually location of the lower 10% vapor fraction region moves further left towards the throat finally causing the detachment of the cavitation cloud as described above for Figure 5. Folding of the shear layer is coupled with its gradual thickening. Location, where the gradient of shear layer folding is high, is accompanied by the regions of main flow deceleration and hence variations of the velocity and pressure in the main flow.

#### **5. Conclusions**

Numerical simulation has provided additional insights into the formation of Kelvin-Helmholtz flow instabilities in the Venturi section. The process is being able to cause the shedding of cavitation clouds similar to the re-entrant jet or shockwave cases. Results conform well to the available experimental data.

The formation of Kelvin-Helmholtz instabilities is determined by the large vertical variations in the shear stress, vapor fraction, and the different physical properties of liquid and vapor regions. The vapor with its much lower density and momentum is more impacted by the instabilities creation as the liquid region. This promotes perturbations and the development of the turbulent eddies in the vapor region. Turbulent eddy impingement point on the shear layer point marks the position of the large variation of the shear stress gradient. We observed the folding of the lines of constant vapor fraction within the shear layer even before any velocity changes were noticed. Thus, the vapor fraction properties may influence the Kelvin-Helmholtz instability formation process through the variation of the density and the transfer of the momentum. As instabilities grow further, adverse pressure gradient and velocity lead to the detachment of the cavitation cloud together with all instabilities.

In both layers, Kelvin-Helmholtz instabilities for the case of cavitation in the Venturi section express less pronounced rolling motion and are less periodic and symmetric in comparison with their more well-known occurrences. This perhaps was the reason that Kelvin-Helmholtz instability was not recognized before as a possible source of the cavitating flow instabilities and cavitation cloud shedding.

More experimental and numerical studies will be required in the future to better understand the formation of Kelvin-Helmholtz flow instabilities in cavitation flow and their role in the cavitation shedding process.

**Author Contributions:** Conceptualization, M.D. and P.P.; methodology, P.P., M.H. and M.D.; validation, P.P., M.H. and M.D.; Investigation, P.P., M.H. and M.D.; resources, M.H. and M.D.; writing original draft preparation, P.P., M.H. and M.D.; writing—review and editing, P.P., M.H. and M.D.; visualization, P.P. and M.D.; supervision, M.H. and M.D.; funding acquisition, M.H. and M.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** European Research Council (ERC) under the European Union's Framework Program for research and innovation, Horizon 2020 (grant agreement n 771567—CABUM), the Slovenian Research Agency (Core funding No. P2-0401).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge the financial support from the European Research Council (ERC) under the European Union's Framework Program for research and innovation, Horizon 2020 (grant agreement n 771567—CABUM), the Slovenian Research Agency (Core funding No. P2-0401) and the PhD grant for young researchers (P.P.).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


*Article*

## **Development of Attached Cavitation at Very Low Reynolds Numbers from Partial to Super-Cavitation**

**Florent Ravelet \* , Amélie Danlos, Farid Bakir, Kilian Croci , Sofiane Khelladi and Christophe Sarraf**

Arts et Metiers Institute of Technology, CNAM, LIFSE, HESAM University, 75013 Paris, France; amelie.danlos@ensam.eu (A.D.); farid.bakir@ensam.eu (F.B.); kilian.croci@ensam.eu (K.C.); sofiane.khelladi@ensam.eu (S.K.); christophe.sarraf@ensam.eu (C.S.)

**\*** Correspondence: florent.ravelet@ensam.eu

Received: 1 October 2020; Accepted: 16 October 2020; Published: 20 October 2020

**Abstract:** The present study focuses on the inception, the growth, and the potential unsteady dynamics of attached vapor cavities into laminar separation bubbles. A viscous silicon oil has been used in a Venturi geometry to explore the flow for Reynolds numbers ranging from *Re* = 800 to *Re* = 2000. Special care has been taken to extract the maximum amount of dissolved air. At the lowest Reynolds numbers the cavities are steady and grow regularly with decreasing ambient pressure. A transition takes place between *Re* = 1200 and *Re* = 1400 for which different dynamical regimes are identified: a steady regime for tiny cavities, a periodical regime of attached cavity shrinking characterized by a very small Strouhal number for cavities of intermediate sizes, the bursting of aperiodical cavitational vortices which further lower the pressure, and finally steady super-cavitating sheets observed at the lowest of pressures. The growth of the cavity with the decrease of the cavitation number also becomes steeper. This scenario is then well established and similar for Reynolds numbers between *Re* = 1400 and *Re* = 2000.

**Keywords:** partial cavitation; super-cavitation; laminar cavitation; cavitation instabilities

#### **1. Introduction**

The presence of small separation bubbles close to the leading edge of a profile or to the summit of a wedge provides favorable conditions for the attachment of vaporous cavities when lowering the absolute pressure. The fundamentals of sheet cavitation inception can be found for instance in the classical books of C. E. Brennen [1] or J.-P. Franc [2]. The inception of such "sheet" cavitation has been for instance studied in Refs. [3–7] on smooth axi-symmetric bodies, on propellers blades or on foils. Moreover, a recent review of the different sheet cavitation inception mechanisms can be found in Ref. [8].

Once the inception of cavitation has occurred, for lower pressures the attached sheet cavities usually grow and become unstable, leading to periodical cloud shedding [9–11]. Partial and cloud cavitations are associated with a vapor region that extends over a part of the cavitating body [2]. At very low pressures, one can observe super-cavitation, consisting of a large and stable vapor cavity that extends beyond the body and closes in the liquid. In that case the pressure fluctuations are usually small [12,13]. For turbomachinery applications, the unsteady dynamics of attached sheets may be responsible for an increase of noise and vibrations and is of crucial importance for system instabilities [14,15]. The instabilities of attached sheet cavitation are attributed to two main phenomena: the formation of a re-entrant jet which is governed by inertia [16], and a bubbly shock propagation mechanism [17]. These two mechanisms have been recently subject to extensive studies, at least

for very large Reynolds numbers, in water [18–22]. The cavitation in that case arises on a turbulent flow background.

The current studies that are conducted in the LIFSE (Laboratory of Fluid Engineering and Energetic Systems) deal with cavitating flows at very low Reynolds numbers, i.e., in laminar or transitional regimes. To the best of our knowledge, only very few studies have been carried out on hydrodynamic cavitation [23,24] or flow-induced degassing [25,26] in laminar flows. The main goal of the present study is to identify and compare potential unsteady regimes in viscosity-dominated cavitating flows to the aforementioned regimes observed in turbulent flows. In the present study, viscous silicone oil is used in a well-documented Venturi-type geometry [10]. In a first study conducted in our laboratory (Croci et al. (2018, 2019) [27,28]), the main features of the multiphase structures that arise have been observed from a qualitative point of view, using silicone oil. Silicone oil can dissolve much more air than water under the same thermodynamical conditions [29,30]. In the work of Croci et al. (2018, 2019) [27,28] the oil was saturated with air at atmospheric pressure. Single-phase numerical simulations have been performed additionally to illustrate the flow topology and to estimate the threshold for cavitation inception. As a result, different types of cavities are observed: tadpoles that attach on the lateral wall [27], sheets that attach on the Venturi slope and top-tadpoles that attach on the opposite side of the vein. The numerical simulations, pressure measurements and high-speed visualizations suggest to attribute the tadpoles to degassing, the main cavity sheet to cavitation, and the top-tadpoles to a secondary flow separation that captures air [28].

In the present article, experiments were performed with degassed silicone oil, to remove air bubbles and to get rid of degassing phenomena [25,26,31]. As a result, only central attached sheet cavities are present in the flow. Quantitative measurements based on high-speed visualizations are performed in order to characterize the evolution of the shapes, lengths and temporal features of the attached sheets with a decrease of the ambient pressure at constant Reynolds numbers *Re*, both in steady laminar flow (at *Re* 6 1200) and in transitional flows (up to *Re* ≃ 2000).

Section 2 is dedicated to the presentation of the experimental setup. The test-bench and the protocol are described in Section 2.1; the main parameters and the post-processing of the high-speed movies are discussed in Section 2.2; and a few numerical results giving insight into the single-phase flow topology and about its transition are recalled in Section 2.3.

The main results are presented in Section 3. The differences between experiments with air-saturated and degassed oil are illustrated in Section 3.1. The different regimes are presented in Section 3.2 and are quantitatively characterized in Section 3.3. Concluding remarks are then given in Section 4.

#### **2. Experimental Study Overview**

#### *2.1. Experimental Test-Bench and Protocol*

The experiments are conducted in a test-bench of the *LIFSE* facilities especially designed to study cavitation in silicone oils (see Figure 1). A Pollard MPLN 142 volumetric pump (1) impels the silicone oil in the loop. A tank of 100 liters (2) is equipped with a temperature sensor ThermoEst PT100 (3). The tank is connected to compressed air supply or to a vacuum pump (4). The oil then flows in a 40 mm inner diameter pipe and through an ultrasonic flowmeter KHRONE Optisonic 3400 (5) placed 1.5 m upstream of the test section (7). The inlet and outlet pressures are monitored with two JUMO dTrans p30 pressure sensors (6).

The test section consists of a Venturi geometry with convergent/divergent angles of respectively 18◦ and 8◦ [10]. It presents an inlet section *<sup>S</sup>in* <sup>=</sup> <sup>20</sup> <sup>×</sup> 10 mm<sup>2</sup> and a square section at the throat *<sup>S</sup>throat* <sup>=</sup> <sup>10</sup> <sup>×</sup> 10 mm<sup>2</sup> . Both side and top visualizations are possible. Two round-to-rectangle contraction nozzles with an area ratio of 6.3 are placed on each sides of the test section to connect it to the pipe system.

All experiments are realized following a protocol to eliminate a maximum of dissolved air. The pressure in the setup is first decreased to an absolute pressure of ≃70 mbar, then the oil is kept

circulating for this pressure at slow velocity during two hours, while removing periodically the air that accumulates in the vertical tube situated above the outlet of the test section and visible on the top part of the sketch on the left in Figure 1. The measurements are then operated at roughly constant Reynolds number, increasing progressively the pressure in the test-bench.

**Figure 1.** Sketch of the Venturi geometry. The width of the test section is *w* = 10 mm and its height at the inlet section is *hin* = 20 mm. In the following of the article the axes origin is positioned at the Venturi throat edge in the middle of test section along the width. Reprinted with permission from K. Croci, Ph.D. dissertation [32].

#### *2.2. Control Parameters, Flow Visualization and Measured Quantities*

The inlet and outlet test section absolute pressure, named *P*<sup>1</sup> and *P*2, are measured with two pressure transducers located respectively at 103 mm upstream and 201 mm downstream of the Venturi throat. The discharge throat velocity *Vin* at *Sin* is computed from flow rate measurements. The viscosity and density have been measured as a function of the temperature [32].

The following dimensionless numbers are defined based on these values: a Reynolds number *Re*, an inlet and an outlet cavitation numbers *σ*<sup>1</sup> and *σ*<sup>2</sup> and a capillary number *Ca*. The definition of these parameters, together with the ranges that have been explored in the present article and the associated uncertainties, can be found in Table 1. Please note that the outlet cavitation number *σ*<sup>2</sup> has not been considered in the previous study [28]. It is used in the present article as it seems to be more significant for analyzing the results. We also introduce a dimensionless pressure coefficient *Kp*. Please note that *K<sup>p</sup>* is simply *K<sup>p</sup>* = *σ*<sup>1</sup> − *σ*2.

Two Optronis CR1000x3 high-speed cameras are used to visualize the cavities. The two cameras allow synchronized pictures from the top and from the side of the test section, with a field of view of 1280 × 512 pixels. Image sequences of 2 to 3 seconds are recorded at 1500 Hz or 1000 Hz. The flow is illuminated from the bottom and the backside with two continuous white LED plate from Phlox and consequently, the gaseous cavities appear as black features in the images. A typical result obtained at *Re* ≃ 1300 is displayed in Figure 2.


**Table 1.** Flow parameters, ranges and estimated uncertainties for silicone oil 47*V*50. The surface tension is *<sup>S</sup>* <sup>≃</sup> 20 mN·m−<sup>1</sup> and the vapor pressure is *<sup>P</sup><sup>v</sup>* <sup>≃</sup> 1.3 Pa according to the oil data file.

**Figure 2.** (**a**) Typical side (left column) and top (right column) instantaneous views. The flow is from left to right. *Re* ≃ 1300, *σ*<sup>1</sup> = 5.10 ± 0.04 and *σ*<sup>2</sup> = 0.70 ± 0.03. (**b**) Illustration of the image processing. The treatment is based on binarization, morphological closing and regions labeling (see text for details). The large region close to the throat is displayed in black in the figure, whereas the other small gaseous cavities are in white and the background in pale gray. A binary image (not shown) is generated that only contain the large region close to the throat. (**c**) Time average of 3000 binary images showing the average attached cavity on the side and top views. The colorbar corresponds to the proportion of time when the cavity is present.

The height of the Venturi throat *hthroat* = 10 mm is chosen as the reference length scale for the dimensionless coordinates *X* ∗ , *Y* ∗ and *Z* ∗ (see also the bottom right part in Figure 1 for the definition of the coordinate system). The origin of the positions is on top of the wedge, in the symmetry plane of the flow vein. The diverging part of the Venturi section thus extends from *X* ∗ = 0 ; *Z* ∗ = 0 to *X* ∗ = 1/tan (8 ◦ ) ≃ 7.1 ; *Z* ∗ = −1; the two side walls are at *Y* ∗ = ±0.5 and the upper flat wall is at *Z* ∗ = 1.

Time series of 3000 images are recorded for each experiment. A background reference image *Ire f* is also recorded with the same cameras setting, while the flow is at rest. The instantaneous images *I<sup>i</sup>* are then processed with Python, using the scipy.ndimage and scikit-image libraries [33] in the following way, based on Morphological Image Processing [34]:


**Figure 3.** Illustration of image processing for the characterization of the temporal features of the attached cavities. (**a**): Instantaneous normalized side view *I* ′ *i* at *Re* = 1600, *σ*<sup>1</sup> = 4.97 ± 0.04 and *σ*<sup>2</sup> = 1.17 ± 0.03. (**b**): Gray level along the dashed blue line displayed in (a); *s* ∗ is the dimensionless curvilinear abscissa along the line. (**c**): Spatio-temporal diagram along the line; space is from top to bottom and time is from left to right. (**d**): Corresponding temporal magnitude spectrum.

The length of the attached cavities is then measured on the mean images of the side view, using a threshold *I <sup>m</sup>* > 0.6. This value corresponds to the region of space where the attached sheet is present on 60% of the time series. For the case illustrated in Figure 2, it gives a length *L* <sup>∗</sup> ≃ 3.5. Varying the threshold between 0.5 and 0.7—see the colorbar in Figure 2c—leads to a variation of ±0.2 on the measure of *L* <sup>∗</sup> which is a relative variation of ±5%.

Spatio-temporal diagrams are also extracted along various lines in the normalized images *I* ′ *i* in order to characterize the dynamics of the cavities and to measure the frequencies of periodical shrinking or shedding of the cavities. The post-processing is illustrated in Figure 3, where the spatio-temporal diagram is extracted along a line parallel to the Venturi diverging wall, about 1 mm above. The spatio-temporal matrix of gray levels is averaged along the spatial coordinate to give a time signal, from which the magnitude spectrum is computed. The main result which is the frequency of the main peak in the spectrum—if existing—is insensitive to the precise choice of the line: similar diagrams and the same frequency are obtained for horizontal or inclined lines at various heights, providing these lines cross the attached sheet. The best signal-to-noise ratio is obtained for the line displayed in Figure 3.

#### *2.3. Numerical Highlights and Comparison with Experimental Data*

The topology of the single-phase flow and its nature as a function of the Reynolds number have been studied numerically [28]. The computation was performed with the Computational Fluid Dynamics code StarCCM+ (version 13.02). It uses a finite volume approach. The computation is three-dimensional; the numerical domain is similar to the experiment with an extrusion of 10*hthroat* upstream of the inlet of the test section and an extrusion of 20*hthroat* downstream of the outlet of the test section. The boundary conditions are a uniform inlet velocity and a uniform static pressure outlet. The domain is meshed with a structured grid. A grid convergence test was made at *Re* = 1500 with up to 7.2 × 10 6 cells. The final mesh contains 1.1 × 10 6 cells, with 51 points in the cross-stream direction and a near-wall cell size of 5 × 10 <sup>−</sup>3*hthroat* . The equations that are solved are the incompressible Navier–Stokes equations for a Newtonian fluid with constant properties. The spatial discretization for the convective flux is the second-order upwind scheme. Unsteady simulations with the Euler first-order implicit scheme have been performed. The main results are plotted in Figure 4. The experimental (obtained in single-phase flows) and numerical pressure loss coefficients fairly collapse, which is a clue in favor of the quality of the simulations (Figure 4a).

**Figure 4.** Main results of the single-phase numerical simulations. (**a**) Pressure loss coefficient *K<sup>p</sup>* as a function of the Reynolds number *Re*. Blue squares represent experimental values while black squares represent numerical simulations. The dashed line is a line of equation *K<sup>p</sup>* ∝ *Re* −1 . (**b**) Dimensionless pressure profile *P* ∗ = *P*/ 1 2 *ρV* 2 *in* with respect to axial coordinate *X* ∗ = *X*/*hthroat* . (**c**) Dimensionless velocity profile *V* ∗ = *V*/*Vin* with respect to *X* ∗ . For both images (**b**,**c**): *Re* = 1200, *σ*<sup>1</sup> = 6.57 and *σ*<sup>2</sup> = 2.58; Blue line is along the horizontal line *Z* ∗ = 0.5 and *Y* ∗ = 0; Red line is a broken line following the bottom of the vein, 0.5 mm above it at *Y* ∗ = 0.

The main results of this numerical study are the following:


#### **3. Results**

#### *3.1. Comparison between Degassed and Air-Saturated Oil*

In the previous study using silicone oil saturated with air at 1 bar (Croci et al. (2019) [28]), three types of gaseous/vaporous structures have been observed. A typical example obtained at *Re* = 1810, *σ*<sup>1</sup> = 4.81 and *σ*<sup>2</sup> = 1.52 is shown in Figure 5a. The first one is called "bottom tadpoles" (see label (I) in Figure 5a). They attach to the side walls close to the throat, in the region of the primary flow separation that has been numerically demonstrated, and can be distinguished only in the top view. These structures appear for *Re* 6 650 and at inlet pressures well above the estimated critical value for cavitation. They thus might be composed of degassed air. The second type of structure consists of an attached cavity in the middle of the vein, just downstream of the throat (see label (II) in Figure 5a). It has been observed for *Re* 6 800 and their inception fairly corresponds to the critical value of the cavitation number estimated with numerical simulations. This attached cavity thus probably contains vapor. The third type of structures is reminiscent of the tadpoles and has been called "top tadpoles" (see label (III) in Figure 5a): they appear for *Re* > 1000 and attach to the upper wall in the secondary boundary layer separation. The pressure in these regions is above the vapor pressure and they might also be composed of degassed air.

**Figure 5.** Comparison between two experiments performed at similar Reynolds numbers (*Re* ≃ 1800) and cavitation numbers (*σ*<sup>2</sup> ≃ 1.5) with a different degassing protocol. (**a**): Silicone oil is saturated with air at 1 bar before running the experiment. The three types of gaseous structures are spotted with (I) for bottom tadpoles, (II) for attached sheet cavity and (III) for top tadpoles. (**b**): Silicone oil is first degassed at 70 mbars; an attached sheet cavity is only present.

In the present article, with a different protocol of intense degassing at 70 mbars before running the experiments, the bottom and top tadpoles are not observed (see Figure 5b, obtained at *Re* = 1810, *σ*<sup>1</sup> = 4.91 and *σ*<sup>2</sup> = 1.57). This confirms that the top and bottom tadpoles are linked to degassing phenomena. Only the central attached cavity remains, and these cavities are moreover observed for cavitation numbers lower than the numerical critical value.

The following paragraphs are devoted to the study of the development of these cavities when varying the absolute pressure at constant Reynolds number.

#### *3.2. Developed Cavitation Regimes Identification in Degassed Oil*

Approximately one hundred of movies have been recorded at various Reynolds numbers and cavitation numbers. Only the cases with no spurious air bubbles coming from upstream and with cavities larger than 1 mm have been reported in Figure 6.

**Figure 6.** Explored map in the Reynolds—cavitation number (*Re* − *σ*2) space. Black circles (•) correspond to type-1 steady small cavities, red squares () to type-2 cavities with periodic shrinking, crosses (×) to type-3 unsteady cavities with no clear periodic behavior and blue diamonds () to type-4 supercavities with *L*\* > 7.1. The solid pale blue line corresponds to the critical cavitation numbers obtained with the single-phase numerical simulations. Images of the cases (a) to (d) are displayed in Figure 7.

Re

The data all fall below the numerical inception curves. Four different regimes have been identified. Typical instantaneous side and top views illustrating these four regimes are displayed in Figure 7. The first regime identified with a full circle (•) corresponds to steady cavities. For *Re* 6 1000, even at very small cavitation numbers, only steady cavities of this type are observed. The cavity viewed from side is quite flat and its bottom part is detached from the divergent bottom wall.

For larger Reynolds numbers (*Re* > 1200), the type-1 steady attached cavities are only observed at high cavitation numbers, corresponding to small cavities. When lowering the cavitation number, a second regime arises: the cavity periodically grows and shrinks, which corresponds to spatio-temporal diagrams looking like the one displayed in Figure 3c. A typical cycle of shrinking and growing is illustrated in Figure 8. This second regime, characterized by a "sawteeth" spatio-temporal diagram is called "type-2 cavities with periodic shrinking" and is displayed with red squares (). In this regime, the shape of the cavity is different and varies during the cycle. It consists of an elongated bubble with a narrow tail at maximum length and during the shrinking, one can observe that the tail of the cavity turns into a vertical front that collapses towards the throat. For the case illustrated in Figure 8, the length of the cavity varies between *L* ∗ *max* ≃ 4 and *L* ∗ *min* ≃ 0.2.

Lowering further the cavitation number, the attached cavity loses this periodic behavior: the spatio-temporal diagrams still show unsteadiness in the flow but with no clear periodic cycle. This third regime is displayed with crosses (×) in Figure 6. The variation of its length corresponds this time to the detachment of delta and hairpin vortices as can be seen in Figure 7c. The cavity in this case keeps a length greater than *L* ∗ *min* ≃ 1 and has a maximum length of the order of *L* ∗ *max* ≃ 6. The characteristic vertical front at the rear part of the cavity that has been noticed in the previous regime is no longer observed.

Finally, at the lowest cavitation numbers, one can obtain an almost steady cavity that appears to be completely transparent from both top and side view, and whose length is longer than the end of the divergent part of the Venturi profile (i.e., *L* ∗ > 7.1). (see Figure 7d). This fourth type is called "supercavities" and is displayed with blue diamonds () in Figure 6.

**Figure 7.** Side and top views in different regimes. (**a**): type-1 steady cavity (•) at *Re* = 992, *σ*<sup>1</sup> = 5.96 and *σ*<sup>2</sup> = 0.51. (**b**): type-2 cavity with periodic shrinking () at *Re* = 1614, *σ*<sup>1</sup> = 4.87 and *σ*<sup>2</sup> = 0.92. (**c**): type-3 unsteady cavity with no clear periodic behavior (×) at *Re* = 1600, *σ*<sup>1</sup> = 4.89 and *σ*<sup>2</sup> = 0.73. (**d**): to type-4 supercavity () at *Re* = 1810, *σ*<sup>1</sup> = 4.66 and *σ*<sup>2</sup> = 0.09.

**Figure 8.** One cycle viewed from side in the type-2 periodic shrinking regime () at *Re* = 1614, *σ*<sup>1</sup> = 4.87 and *σ*<sup>2</sup> = 0.92. The measured frequency of the cycle is 34.6 *Hz*.

#### *3.3. Quantitative Characterization of the Different Regimes and Effects of the Reynolds Number*

The type-2 cavities with periodic shrinking are first observed in a very narrow region of the parameter space for 1200 6 *Re* 6 1350 and are then present in a wider range of cavitation numbers for *Re* ∈ [1400; 2100]: they are found to lie in *σ*<sup>2</sup> ∈ [0.9; 1.7].

The type-4 supercavities have been observed only for one point at very low pressure at *Re* ≃ 1300 (*σ*<sup>2</sup> ≃ 0.27). If supercavities exist for *Re* = 1200 it would be at *σ*<sup>2</sup> 6 0.3. However, for higher Reynolds numbers, i.e., for *Re* > 1400 the zone where supercavities are observed progressively extends up to a constant value of the cavitation number *σ*<sup>2</sup> 6 0.57.

There thus seems to be a progressive transition in the dynamics of attached sheet cavitation at very low Reynolds number between *Re* ≃ 1200 and *Re* ≃ 1400. Quantitative measurements of the temporal and spatial features of the attached cavities are analyzed and discussed in the next paragraph.

First, the frequency of the periodic shrinking for all the type-2 cavities that have been identified ( in Figure 6) are plotted as a function of the throat velocity in Figure 9. The two quantities are highly correlated: the frequency seems to be proportional to the throat velocity, following a law *f* = 3.83*Vthroat* . A Strouhal number based on the throat velocity and on the throat height as a length scale is *S<sup>t</sup>* = *f hthroat Vthroat* . The type-2 regime is thus characterized by a constant Strouhal number of the order of *S<sup>t</sup>* ≃ 0.0026. Such a low value is reminiscent of the "sheet cavity oscillation" regime described for a turbulent flow in a similar water experiment by Danlos et al. (2014) [10] as opposed to the "cloud cavitation" regime that is characterized by a Strouhal number greater by almost one order of magnitude. Moreover, in that study [10] the sheet cavity oscillation was observed to dominate the dynamics for short cavities and to disappear as the cavity length increases.

The evolution of the average length of the cavities as a function of the outlet cavitation number is plotted in Figure 10. For all Reynolds numbers, the overall trend is an increase of the length with a decrease of the cavitation number, which is not surprising. However, the shape of the curve exhibits a dependence on the Reynolds number. One can indeed notice that the curve becomes steeper as the Reynolds number increases, from *Re* ≃ 800 to *Re* ≃ 1300. At very low Reynolds number the cavity grows relatively slowly with a decrease of the cavitation number and the transition to super-cavitation is progressive at *Re* ≃ 1300. The data obtained for *Re* > 1400 then follow a different trend, with first a slow growing until *L* <sup>∗</sup> ≃ 1 and then a faster growing of the cavity. Moreover, one can notice that for *Re* > 1300, the transition between the type-2 regime () and the type-3 regime (×) corresponds to an average length *L* <sup>∗</sup> ≃ 2, around *σ*<sup>2</sup> ≃ 1.

**Figure 9.** Measured frequency as a function of the throat velocity for type-2 regime. The dashed line is a linear fit of equation *f* = 3.83*Vthroat* .

**Figure 10.** Dimensionless average length of the attached cavity *L* ∗ as a function of the outlet cavitation number *σ*2. For *Re* > 1300, the symbols refer to the regimes identified in Figure 6: type-1 steady cavity (•), type-2 cavity with periodic shrinking (), type-3 unsteady cavity with no clear periodic behavior (×) and type-4 supercavity (). The lines are a guide to the eye based on spline interpolation of the points at *Re* = 800, *Re* = 1000, *Re* = 1200 and of all the points for *Re* > 1400.

2

The dimensionless pressure drop *K<sup>p</sup>* in the experiment without cavitation, i.e., at *σ*1,2 >> 1, is plotted in Figure 4a as a function of the Reynolds number. Let us recall that *K<sup>p</sup>* = *σ*<sup>1</sup> − *σ*2. As the cavitation develops while decreasing the ambient pressure at constant discharge velocity, the presence of gaseous structures of increasing volume may change the pressure drop. It might thus be interesting to explore the variation of both the inlet and outlet pressure, i.e., in a dimensionless form of both *σ*<sup>1</sup> and *σ*2. The experiments are plotted in this parameter space in Figure 11. The points at constant Reynolds number are plotted with the same symbol.

The development of cavitation associated with a decrease of the ambient pressure corresponds to a path from top right to bottom left in this parameter space. For the seven Reynolds numbers that are displayed in the Figure 11, one can first notice that the points follow a similar trend: the two pressures are decreasing together. Moreover, one can distinguish two main different shapes in the curves. On the one hand, for *Re* ≃ 1000 and *Re* ≃ 1200, the points align quite on a straight path. On the other hand, for *Re* > 1300, one can notice some kind of break in the slope of the curve around *σ*<sup>2</sup> ≃ 1. For *σ*<sup>2</sup> > 1, the slope is of the order of magnitude of the former (roughly −2), whereas for *σ*<sup>2</sup> 6 1, the slope is steeper, of the order of −10: a very small decrease in the inlet pressure causes a

dramatic decrease of the outlet pressure corresponding to the fast increase of the cavity length observed in Figure 10.

**Figure 11.** Explored map in the inlet cavitation number—outlet cavitation number (*σ*<sup>1</sup> − *σ*2). The different symbols stand for different Reynolds numbers. The arrows are guides to the eye.

These observations support the idea of a transition in the cavitating flow that takes place around *Re* ≃ 1300.

Looking at the attached cavities from top (see the right column in Figures 2, 5 and 7), one can notice that the morphology of the cavity also exhibits a dependence in the transverse direction *Y* ∗ . At very low Reynolds number (see Figure 7a) the interface at the leading edge presents a single smooth and convex cap. At larger Reynolds numbers (see Figure 7d) this interface presents several indentations. These structures resemble the "divots" that have been discussed for instance in Ref. [35] and that are attributed to a local breakdown of the two-dimensionality of the laminar boundary layer. We report in Figure 12 the number of indentations that are observed on the cavities of an average length of *L* <sup>∗</sup> > 4. The first indentation is observed at *Re* ≃ 1200 and then their number increases almost linearly with the Reynolds number.

**Figure 12.** Number of transverse indentations at the leading edge of the developed attached cavities as a function of the Reynolds number.

Re

#### **4. Conclusions**

In the present article, the development of attached cavities under cavitating conditions in a highly degassed silicone oil has been explored, at Reynolds numbers in the range *Re* ∈ [800; 2000]. The regime map in terms of Reynolds and cavitation numbers have been established. Four different regimes have been observed: type-1 small cavities that are almost steady, type-2 periodically shrinking cavities of intermediate sizes, type-3 vortex bursting with no clear period and type-4 supercavities. For *Re* 6 1000, the background flow is steady and laminar and the cavitating flow gives rise to steady attached cavities whose length slowly increases with a decrease of the cavitation number. The first footprint of unsteadiness in the cavitating flow consists of the type-2 periodic regime. It appears for a sufficiently great Reynolds number (*Re* > 1200), first in a very narrow range of cavitation numbers and then in a wider range as *Re* > 1400. The frequency is almost proportional to the velocity of the upstream flow and gives a constant Strouhal number of very small value. More modal analysis with POD (Proper Orthogonal Decomposition) or DMD (Dynamic Mode Decomposition) is in progress to characterize better this regime. These experiments will require higher temporal resolution. The analysis of the length of the cavity, of the transition to super-cavitation and of the evolution of the pressure drop are consistent with a smooth transition in the cavitating flow for 1200 > *Re* > 1400, which corresponds to the range where the single-phase (non-cavitating) flow is transitional. It would be interesting to extend further the range of Reynolds numbers, using an oil of lower viscosity, in order to give a lower bound in Reynolds number for the periodical cloud shedding regime that is known to take place in this geometry in fully developed turbulent flows.

**Author Contributions:** Conceptualization, F.R., F.B., A.D. and S.K.; methodology, F.R.; software, F.R.; validation, F.R., F.B., A.D. and S.K.; formal analysis, F.R.; investigation, F.R. and A.D; resources, K.C. and C.S.; data curation, F.R. and A.D.; writing–original draft preparation, F.R.; writing–review and editing, F.R.; visualization, F.R.; supervision, F.B. and S.K.; project administration, F.R.; funding acquisition, F.B. and S.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors acknowledge Marc Joulin and Jocelyn Mistigres for technical support, and Michaël Pereira for fruitful discussion.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Compressible Two-Phase Viscous Flow Investigations of Cavitation Dynamics for the ITTC Standard Cavitator**

**Ville M. Viitanen 1,\* , Tuomas Sipilä <sup>2</sup> , Antonio Sánchez-Caja <sup>1</sup> and Timo Siikonen <sup>3</sup>**


Received: 4 September 2020; Accepted: 28 September 2020; Published: 7 October 2020

**Featured Application: The present method can be applied to study characteristics of hydrodynamic cavitation and its consequences, such as performance variations and noise, in propellers and rotating turbomachinery. Compressible formulation of the numerical method ensures a more precise treatment of the multiphase flow dynamics, and can reveal some interesting physical phenomena due to the cavitation. Accurate prediction of the multiphase flow is a requisite for optimal, efficient and environmentally friendly design and operation of such devices.**

**Abstract:** In this paper, the ITTC Standard Cavitator is numerically investigated in a cavitation tunnel. Simulations at different cavitation numbers are compared against experiments conducted in the cavitation tunnel of SVA Potsdam. The focus is placed on the numerical prediction of sheet-cavitation dynamics and the analysis of transient phenomena. A compressible two-phase flow model is used for the flow solution, and two turbulence closures are employed: a two-equation unsteady RANS model, and a hybrid RANS/LES model. A homogeneous mixture model is used for the two phases. Detailed analysis of the cavitation shedding mechanism confirms that the dynamics of the sheet cavitation are dictated by the re-entrant jet. The break-off cycle is relatively periodic in both investigated cases with approximately constant shedding frequency. The CFD predicted sheet-cavitation shedding frequencies can be observed also in the acoustic measurements. The Strouhal numbers lie within the usual ranges reported in the literature for sheet-cavitation shedding. We furthermore demonstrate that the vortical flow structures can in certain cases develop striking cavitating toroidal vortices, as well as pressure wave fronts associated with a cavity cloud collapse event. To our knowledge, our numerical analyses are the first reported for the ITTC standard cavitator.

**Keywords:** hydrodynamic cavitation; compressible two-phase flow; turbulence modelling; system instabilities

#### **1. Introduction**

Cavitation, when it occurs, is a cause of many detrimental effects on marine propellers and turbomachinery. Apart from performance degradation, it may trigger erosion, cause vibration problems, or result in high-intensity noise emission on a broad frequency range. Different types of cavitation involve complex features such as bubbly flow, bubble coalescence, nucleation, free-surface instability, and turbulence. These features do not necessarily obey a single but multiple length and time scales, which renders their reliable prediction, either experimentally or numerically, a major undertaking. For instance, when investigating the performance of marine propellers, depending on their operating conditions, many different forms of cavitation can take place simultaneously, such as steady attached sheet cavitation near the blade's leading edge, cloud cavitation in the propeller's wake, and vortex cavitation originating from the blade tips and the propeller hub.

Even in the case of steady uniform inflow conditions, we can observe unsteady behavior of the cavitation. In particular, the shedding behavior of sheet cavitation (Figure 1) and its subsequent transformation to cloudy and bubbly cavities are important in the generation of the aforementioned erosion and noise issues, since upon collapsing these cavitation types can yield strong pressure fluctuations and even shock waves in the bubbly mixture or in the liquid phase [1,2]. A sketch of the shedding cycle of two-dimensional sheet cavitation is shown in Figure 1. The first frames depict the growth of the attached sheet cavity following the break-off of the previous sheet. The growth will cease after a certain length of the sheet cavity, and a strong re-entrant jet is formed at the closure of the sheet cavity, as depicted in Frame 3. The following frame shows the impingement of the cavity interface by the jet, as the jet reaches the cavity near the leading edge. The rear part of the sheet cavity breaks off and reduces to a bubble cloud, which is convected downstream as shown in Frame 5. As the cloud travels downstream, parts of it may undergo a violent collapsing stage, possibly causing noise and material erosion. Meanwhile, as seen in Frame 6, the newly formed sheet cavity on the foil behind the LE will begin to grow, initiating a subsequent break-off cycle.

**Figure 1.** Sketch of the shedding mechanism of sheet cavitation. The figure illustrates one break-off cycle, adapted from [3].

The dynamics of the sheet cavitation are mainly determined by the re-entrant jet, although a bubbly shock wave can in certain conditions affect this development [1]. If the sheet cavity is thick, the break-off cycle is relatively periodic with approximately constant shedding frequency. In these cases, a Strouhal number *St* = *f lre f* /*Ure f* can be determined. The length scale of the Strouhal number is chosen as the maximum length of the cavity. This can be argued to be a typical length scale of the cavity since the shedding frequency depends on the distance the re-entrant jet has to travel. The choice of the velocity scale is somewhat problematic, but can be taken as the free-stream velocity *Ure f* = *U*∞, or based on an approximation of the re-entrant jet velocity *Ure f* = *U*<sup>∞</sup> √ 1 + *σ*, where *σ* is the cavitation number, which acknowledges the presence of a cavity structure in a better way than *U*∞ alone [4]. The frequency, *f* , is based on global shedding of the cavity clouds, determined with a visual inspection of the cavity shedding periodicity, with the help of time histories of the measured lift and drag forces. Typically, the Strouhal numbers for sheet-cavitation shedding vary as 0.20 . . . 0.50, depending on the investigated conditions and the definition for the reference quantities [4–7].

In frames 5–6 of Figure 1 we observe that the detached cavity clouds convect downstream toward a higher-pressure region. Upon reaching the regions of sufficiently higher ambient pressure, the detached cavity clouds likely undergo a violent collapsing stage. Collapse events of cavitation structures behind a ship's propeller, for instance, can lead to very high-pressure peaks, and such peaks can exhibit very rapid temporal variation. As noted by Matusiak [8], a primary source of high-frequency noise is the collapse of the free cavitation bubbles. The shedding behavior is mainly governed by the development of the re-entrant jet, vortical structures and turbulent flow near the foil, i.e., the unsteady dynamics of especially the liquid phase. As Reynolds-averaged Navier–Stokes (RANS) equations are mainly derived for steady single-phase flows, complex unsteady phenomena may be severely under-predicted when using an unsuitable turbulence-modelling approach. High-fidelity turbulence modelling is typically needed to capture a detailed flow field for the prediction of the dynamics of cavitation and its consequences such as erosion and noise [9–12]. Therefore in this study, we employ two different turbulence closures: an unsteady RANS-based SST *k* − *ω* turbulence model and a hybrid RANS/LES (large eddy simulation) closure called the delayed-detached eddy simulation (DDES). We analyze the two-phase flow using a compressible multiphase model. Cavitating problems have been typically treated with incompressible flow solvers [13–17] (especially in the case of marine propellers), and recently either isothermal compressible or fully compressible flow formulations have been used as well [18–21]. If we aim for a precise treatment of hydrodynamic cavitation in its various forms, a compressible multiphase flow model is a requisite for reliable numerical predictions [18]. Peculiarly, the sound speed in a mixture can drop to a fraction of that in the liquid and even the vapor phases. Prediction of the correct acoustic signal speeds calls for a complete flow model.

In this paper, the dynamics of the unsteady sheet cavitation is studied in detail. We employ a compressible two-phase flow model to account for the cavitation phenomena, and the computational case is the ITTC standard cavitator [22]. The results are compared to the model tests carried out in a cavitation tunnel. The flow simulations have been conducted with the general purpose CFD solver FINFLO. The cavitation tunnel tests, which include hydrophone measurements of the underwater noise, have been conducted at SVA Potsdam. We have carried out simulations of two different cases. The inflow velocity for both cases was 8 m/s, and the cavitation numbers were *σ* = [0.53, 1.27]. To the knowledge of the authors, our results include the first reported numerical studies of the ITTC standard cavitator.

This paper is organized as follows. Section 2 describes the flow solver used, as well as the turbulence and multiphase flow modelling. In Section 3 we describe the test case. Section 4 presents the results in terms of cavitation observation and its dynamics, and their comparison to the experimental findings. Lastly, we present conclusions and suggestions for future work in Section 5.

#### **2. Flow Solution**

#### *2.1. Governing Equations*

Our flow model is based on a compressible form of the RANS equations where an assumption is made on the homogeneity of the fluid mixture among the different phases involved in the computation [23,24]. Including the energy equations for the solution allows prediction of the correct acoustic signal speeds. The governing continuity and momentum equations are

$$\begin{aligned} \frac{\partial a\_k \rho\_k}{\partial t} + \nabla \cdot a\_k \rho\_k \mathbf{V} &= \Gamma\_{k'}\\ \frac{\partial \rho \mathbf{V}}{\partial t} + \nabla \cdot \rho \mathbf{V} \mathbf{V} + \nabla p &= \nabla \cdot \boldsymbol{\tau}\_{ij} + \rho \mathbf{g}\_{\mathbf{v}} \end{aligned} \tag{1}$$

where *p* is the pressure, **V** the absolute velocity in a global non-rotating coordinate system, *τij* the stress tensor, *α<sup>k</sup>* a void (volume) fraction of phase *k*, *ρ<sup>k</sup>* the density, *t* the time, Γ*<sup>k</sup>* the mass-transfer term, and **g** the gravity vector. The void fraction is defined as *α<sup>k</sup>* = V*k*/V , where V*<sup>k</sup>* denotes the volume occupied by phase *k* of the total volume, V . For the mass transfer, ∑*<sup>k</sup>* Γ*<sup>k</sup>* = 0 holds, and consequently

only a single mass-transfer term is needed. The energy equations for phase *k* = *g* or *l*, the indices referring to gas (*g*) and liquid (*l*) phases, are written as

$$\begin{split} \frac{\partial a\_k \rho\_k (e\_k + \frac{V^2}{2})}{\partial t} + \nabla \cdot a\_k \rho\_k (e\_k + \frac{V^2}{2}) \mathbf{V} &= \\ - \nabla \cdot a\_k \mathbf{q}\_k + \nabla \cdot a\_k \mathbf{r}\_{lj} \cdot \mathbf{V} + q\_{lk} + \Gamma\_k (h\_{ksat} + \frac{V^2}{2}) + a\_k \rho\_k \mathbf{g} \cdot \mathbf{V} . \end{split} \tag{2}$$

Here, *e<sup>k</sup>* is the specific internal energy, **q***<sup>k</sup>* the heat flux, *qik* interfacial heat transfer from the interface to phase *k*, and *hk*sat saturation enthalpy. The interfacial heat transfer coefficients are based on the mass-transfer model by assuming a saturated temperature for the gas phase, and the heat transfer is then determined by the mass-transfer model. In the case of a homogeneous flow, the resulting sound speed *c* for a two-phase mixture is obtained from the eigenvalues of the Jacobian of the flux vector as

$$\frac{1}{\rho c^2} = \frac{\kappa}{\rho\_\mathcal{S} c\_\mathcal{S}^2} + \frac{1-\kappa}{\rho\_l c\_l^2} \quad \text{and} \quad \frac{1}{c\_k^2} = \frac{\partial \rho\_k}{\partial p} + \frac{1}{\rho\_k} \frac{\partial \rho\_k}{\partial h\_k}. \tag{3}$$

In the expressions above, *h<sup>k</sup>* denotes the enthalpy of phase *k*.

The momentum and total continuity equations in the homogeneous model do not change, except for the material properties like density and viscosity. They are calculated for the mixture as a weighted sum by the phasic volume fractions. The turbulence effects are those of single-phase models for the mixture with material properties derived from the pressure and temperatures of the individual phases. The determination of the material properties is described in [25].

#### *2.2. Turbulence Modelling*

Nominally, a Reynolds-averaged form of the Navier–Stokes equations is used, and the DDES approach [26] that combines RANS and LES is also applied in the same form. In this study, the calculations are performed up to the wall both in the model- and full-scale simulations, which avoids the use of wall functions. The height of the first cell was adjusted such that the non-dimensional wall distance *y* <sup>+</sup> = *ρuτy*/*µ* . 1 for the first cell, with *u<sup>τ</sup>* = p *τw*/*ρ* being the friction velocity, *τ<sup>w</sup>* the wall shear stress, and *y* is the normal distance from the solid surface to the center point of the cell next to the surface.

The base turbulence closure applied in the present calculations is the SST *k* − *ω* model [27]. That is a zonal model, referring to the formulation where the *k* − *ω* equations are solved only inside the boundary layer, and the standard *k* − *ε* equations, transformed to the *ω*-formulation, are solved away from the walls. The DDES is also based on the SST model [28], and described in [12]. DDES is a slightly modified version of the detached-eddy simulation (DES). Both are hybrid RANS/LES models, and function as an LES subgrid-scale model in regions where the local turbulent phenomena are of greater size than the local grid spacing, and reduce to a RANS model in regions where the largest turbulent fluctuations are of a smaller size than the local grid spacing [29]. A time-accurate solution is made to resolve turbulent fluctuations.

#### *2.3. Mass and Energy Transfer*

Several mass-transfer models have been suggested for cavitating problems. Usually, the mass-transfer rate is proportional to a pressure difference from a saturated state or to a square root of that. In this study, we employ a mass-transfer model that is similar to that of Choi and Merkle [30], described in [24,31]. The empirical parameters of the cavitation model are calibrated by Sipilä [3]. The saturation pressure is based on the free-stream temperature, and the gas phase is assumed to be saturated.

#### *2.4. Solution Algorithm*

The governing equations are discretized in a finite volume manner with viscous fluxes and pressure terms centrally differenced. For the convective part, a third-order upwind-biased monotonic upstream-centered scheme for conservation laws (MUSCL) interpolation is used for the variables on the cell surfaces. A compressive flux limiter that has been shown to improve the predicted tip vortex cavitation pattern [24] is employed for the void fraction in the present study. For the time derivatives, a second-order three-level fully implicit method is used.

The solution method is a segregated pressure-based algorithm where the momentum equations are solved first, and then a pressure-velocity correction is made. The basic idea in the solution of all equations is that the mass balance is not forced at every iteration cycle, but rather the effect of the mass error is subtracted from the linearized conservation equations. A pressure correction equation was derived from the continuity equation, Equation (1), linked with the linearized momentum equation. The velocity-pressure coupling is based on the corresponding algorithm for a single-phase flow [23], and is described in more detail in [12,24]. A physical time step of <sup>∆</sup>*<sup>t</sup>* = 2.5 × 10 −5 s is used in the simulations, which resolves the cavitation dynamics cycles in the investigated cases with about 250 (*σ* = 1.27) and 1200 (*σ* = 0.53) steps depending on the conditions. Approximately 100 inner iterations were required within each physical time step.

#### **3. Test Case**

A schematic view of the ITTC standard cavitator is shown in Figure 2, and the cavitator inside the cavitation tunnel is seen in Figure 3. The rectangular cross-section of the tunnel has the dimensions of <sup>600</sup> <sup>×</sup> <sup>600</sup> mm<sup>2</sup> . Two cases with different cavitation numbers were investigated. The inflow velocity for both cases was *U*∞ = 8 m/s, and the cavitation numbers were

*σ* = *p* − *psat* 1 2 *ρ*∞*U*<sup>2</sup> ∞ = [0.53, 1.27] .

**Figure 2.** Schematic view of the cavitator.

Three grid densities were used in the calculations. The details are summarized in Table 1. The structured grid has an O-O topology around the foils. The finest grid level consists of roughly 7.3 million cells in 80 grid blocks. The grid resolution around the leading and trailing edges is reasonable, and there are about 15 cells around the leading-edge radius. Due to the O-O topology, the same resolution is applied around the blade tip and the trailing edge as well. In this work, we emphasized the resolution of the sheet cavitation on the wings, and for this reason no clustering of cells was targeted for the cavitating tip vortex, for instance. The grid is refined normal to the viscous surfaces such that *y* <sup>+</sup> <sup>≈</sup> 1 at the finest grid level. A detail of the grid is shown in Figure 4, and the grid resolution on the surface of the foil is shown in Figure 5.

**Figure 3.** Placement of the cavitator inside the cavitation tunnel.

**Table 1.** Summary of the details of the grid used in the calculations for the three grid levels. The number of surface cells are given for one side of the foil.


**Figure 4.** A detail of the grid near the foil.

**Figure 5.** Grid resolution on the surface of the foil. Fine grid.

The foils and the body are modelled as no-slip solid surfaces. A velocity boundary condition is applied at the inlet, and a pressure boundary condition is applied at the outlet. Slip boundary condition is applied at the tunnel walls. The inflow velocity is set based on the case considered, and the background pressure level is set to either a large value for the wetted flow simulations, O(10 6 ), or set based on the cavitation numbers in the cavitating cases. The calculations can be typically performed on three grid levels. On the coarser grid levels, every second point is taken into account compared to the finer level grid. The computations on a coarser grid are used as the initial guess for the computations performed on a finer grid. Initially, we made the unsteady computations for the coarse and medium size meshes to keep the computational time within reasonable limits. However, for the *σ* = 1.27 case, the cavitation extent is not as wide as in the *σ* = 0.53 case, and the dynamics appear at a higher frequency. To capture better cavitation dynamics, we considered it adequate to also use a finer mesh in the simulations.

#### **4. Results**

#### *4.1. Cavitation Observation*

Cavitation observations made in the experiments in the case with *U*∞ = 8 m/s and *σ* = 0.53 are shown in Figure 6, and a snapshot of the simulations with the SST model is given in Figure 7, where the pressure coefficient is defined as *C<sup>p</sup>* = 2*pdi f* /*ρ*∞*U*<sup>2</sup> <sup>∞</sup>, where *pdi f* denotes the pressure difference. Generally, attached sheet cavitation covers roughly half of the suction side in the experiments with cavitation incepting at the leading edge. Almost the entire suction side is covered by unstable detached sheet cavitation in the simulations. Strong and periodic shedding of the sheet cavitation are observed in the simulations, with cloud-like cavities shedding to the wake of the foils. In the experiments, the sheet cavitation is observed to break off as cloud cavitation. We observed streak cavitation at the hemispherical head of the cavitator body in the experiments.

**Figure 6.** Cavitation extent observed in the experiments with *U*∞ = 8 m/s and *σ* = 0.53. Left picture shows a photograph of the head of the cavitator in the experiments, and right picture shows a sketch of the mean cavitation extent on the cavitator. The pictures are reproduced with permission from Rhena Klose (SVA).

In the simulations, however, we see attached sheet cavitation at the hemispherical head of the cavitator body. Strong hub vortex cavitation extending far behind the cavitator body is also observed both in the experiments as well as in the simulations. We notice that in the experiments, a fine tip vortex cavitation extends roughly a chord length in the slipstream, while in the simulations tip vortex cavitation is not visible. In the simulations, strong disturbances created by the unsteady suction side sheet cavities are observed both in the pressure coefficient as well as in the non-dimensional *v*-velocity. In particular, a collapse event of a cavity cloud gave rise to a considerable pressure peak and a subsequent shock-wave-like front in the wake of the foil. The disturbances in *C<sup>p</sup>* and *v*-velocity, however, dissipate quickly near the end of the cavitator body as the grid resolution deteriorates.

**Figure 7.** A snapshot of the simulations with *U*∞ = 8 m/s and *σ* = 0.53. Fine grid with the SST turbulence model. Top frame presents an overview of the cavity extent with the transparent iso-surface of *α* = 0.1 colored by blue. Lower left frame shows the pressure coefficient on the center *z* plane of the cavitator. Lower right frame shows the non-dimensional *v*-velocity on the center *z* plane of the cavitator.

Cavitation observations made in the experiments in the case with *U*∞ = 8 m/s and *σ* = 1.27 are shown in Figure 8, and a snapshot of the DDES is shown in Figure 9. Generally, sheet cavitation covers roughly half of the suction side in the experiments, and roughly one third of the suction side in the simulations with cavitation incepting at the leading edge. Periodic shedding of the sheet cavitation into smaller cavity structures is observed in the simulations. These smaller cavities convected a distance along the foil before disappearing. In the experiments, parts of the sheet cavitation are observed to break off as cloud cavitation. A cavitating hub vortex is visible both in the experiments as well as in the simulations. We noticed that in the experiments, a very fine tip vortex cavitation extends roughly a chord length in the slipstream, while in the simulations tip vortex cavitation is not visible. In the simulations, disturbances created by the unsteady suction side sheet cavities are observed both in the pressure coefficient as well as in the non-dimensional *v*-velocity. A collapse event of the smaller cavity structures gave rise to a pressure peak and a subsequent shock-wave-like front in the wake of the foil. The disturbances, however, dissipate quickly near the end of the cavitator body as the grid resolution deteriorates.

**Figure 8.** Cavitation extent observed in the experiments with *U*∞ = 8 m/s and *σ* = 1.27. Left picture shows a photograph of the right foil of the cavitator in the experiments, and right picture shows a sketch of the mean cavitation extent on the cavitator. The pictures are reproduced with permission from Rhena Klose (SVA).

**Figure 9.** A snapshot of the simulations with *U*∞ = 8 m/s and *σ* = 1.27. Fine grid DDES. Top frame presents an overview of the cavity extent with the transparent iso-surface of *α* = 0.1 colored by blue. Lower left frame shows the pressure coefficient on the center *z* plane of the cavitator. Lower right frame shows the non-dimensional *v*-velocity on the center *z* plane of the cavitator.

#### *4.2. Global Forces*

The mean values of the lift and drag coefficients were extracted from the simulations. In addition, we carried out a Fourier analysis for the lift and drag forces on the wings to determine the unsteady characteristics of especially the sheet-cavitation shedding. The results are shown in Table 2 for all investigated cases.


**Table 2.** Comparison of the mean lift and drag coefficients (*C*L and *C*D , respectively), and their dominant frequencies (*f*0).

We observed a dominant frequency of around 30 Hz for the case with and *σ* = 0.53. We show an example of the time histories and corresponding frequency contents in Figure 10 obtained from DDES, and in Figure 11 obtained from the SST simulations. A frequency of 29.0 Hz for both the lift and drag coefficients was seen with coarse and medium grids when using the SST model. A dominant frequency of 29.0 Hz was observed with coarse grid DDES for both the lift and drag coefficients, while a dominant frequency of 33.4 Hz was observed for lift and 35.0 Hz was observed for the drag with medium grid DDES. Animations from fine and medium grid DDES and SST simulations

revealed a frequency of approximately 30 Hz for the sheet-cavitation shedding. Comparing the DDES and SST simulations, we see that the shedding frequency increases slightly with DDES on the medium grid. The time histories of the SST simulations appear slightly more regular, which is seen also as the sharper frequency peaks whereas these are more spread when using DDES. If we approximate the reference length of the Strouhal number as the chord length, the Strouhal numbers' range is *St* = *f c*/*U*<sup>∞</sup> = [0.22 . . . 0.26]. If we instead take the reference velocity as *Ure f* = *U*<sup>∞</sup> √ 1 + *σ*, the Strouhal numbers would range as *St* = [0.18 . . . 0.21].

**Figure 10.** Time history and frequency content of the lift and drag coefficients. Medium grid, *σ* = 0.53, DDES.

**Figure 11.** Time history and frequency content of the lift and drag coefficients. Medium grid, *σ* = 0.53, SST simulations.

For the case with *σ* = 1.27 a dominant frequency in between 159.6–160.2 Hz for both the lift and drag coefficients was observed with medium and fine grid DDES. As we discuss later in Section 4.6, the SST model yielded a steady-state solution for this case. The time histories and corresponding frequency contents are shown in Figure 12. Also, the animations created from fine and medium grid DDES revealed a frequency of approximately 160 Hz for the sheet-cavitation shedding. If we approximate the reference length of the Strouhal number as a quarter of the chord length, the Strouhal number would be *St* = *f c*/*U*<sup>∞</sup> = 0.30. If instead we take the reference velocity as *Ure f* = *U*<sup>∞</sup> √ 1 + *σ*, the Strouhal number would be *St* = 0.20. Compared to the case with the lower cavitation number, we now see more distinct peaks even up to the 4th harmonics of the lift and drag.

**Figure 12.** Time history and frequency content of the lift and drag coefficients. Fine grid, *σ* = 1.27, DDES.

We note that the conducted simulation time was not enough with the fine grid for the case with *σ* = 0.53 to provide a sufficient time history to conduct a reliable Fourier analysis. The resolution on the coarse grid was not enough to obtain a reasonable frequency for the cavitation shedding with *σ* = 1.27.

There was a small deviation between the mean lift coefficients of the left and right foils with the same cavitation number, while the mean drag coefficients were the same for both foils on each grid. We noted small differences also in the frequencies exhibited in the time histories. A decrease of approximately 40% of the mean lift and an increase of approximately 40% of the mean drag are seen in the results with *σ* = 0.53 compared to the results of *σ* = 1.27. This is because there is considerably more suction side sheet cavitation for the case with *σ* = 0.53, since the cavitation number is lower. The increase in the amount of the vapor phase on the suction side in turn decreased the lift and increased the drag of the foil. This difference was visible in the results using each of the three different grid densities.

#### *4.3. Analyses of Cavitation Dynamics*

Next, we compare the simulations to the experiments in Section 4.3.1 in terms of the hydrophone measurements. The dynamics of the investigated cases with different cavitation numbers are then analyzed in Sections 4.3.3 and 4.3.4. Cavitation collapse events and resulting pressure variations are investigated in Section 4.4. Later in Section 4.5 we demonstrate a formation of a cavitating ring vortex, and in Section 4.6 a comparison of the DDES and SST simulations is shown.

#### 4.3.1. Comparisons with Hydrophone Measurements

Different hydrophone configurations were tested in the cavitation tunnel of SVA Potsdam, depicted in Figure 13 and summarized in Table 3. The hydrophones were installed on the starboard side of the cavitation tunnel window.


**Table 3.** Hydrophone arrangements in the cavitation tunnel.

**Figure 13.** Hydrophone arrangement in the starboard side of the cavitation tunnel. Figure adapted from [32], and reproduced with permission from Rhena Klose (SVA).

Figure 14 shows the results of the conducted hydrophone measurements. In the figures, we show six different hydrophone arrangements. HF1 and HF2 denote the hydrophones in an acoustic chamber, with HF1 mounted in the wall of the cavitation tunnel at a distance of 0.500 m from the origin of the cavitator (see Figure 13), and HF2 mounted a distance of 0.414 m from the cavitator. HF3 and HF4 are located longitudinally to the flow, and the distances to the cavitator are 0.543 m for HF3 m and 0.352 m for HF4. HF5 and HF6 are located transversally to the flow, and the distances to the cavitator are 0.522 m for HF5 m and 0.420 m for HF6. The sensors HF1, HF3 and HF5 are located upstream of the cavitator, and the sensors HF2, HF4 and HF6 are located downstream of the cavitator.

Figure 14a shows the case with *σ* = 0.53, and Figure 14b shows the case with *σ* = 1.27. In both instances, the predicted sheet-cavitation shedding frequencies and their three harmonics are shown as well, cf. Section 4.2. We observe that in the case of *σ* = 0.53, the shedding frequency is picked up by hydrophones HF1, HF2, HF3 and HF4. Moreover, the measurements conducted with HF2, HF3 and HF4 show a smaller peak at the second harmonic of the sheet-cavitation shedding frequency. In the case of *σ* = 1.27, we see that the CFD predicted cavitation shedding frequency peak is visible only in the HF1 measurements. Other frequencies are also seen in the hydrophone measurements that are of greater magnitude than the ones induced by cavitation dynamics in the CFD simulations. For instance, the shedding cavities formed into cloud-like cavitation structures both in the experiments as well as in the simulations, in both the investigated operation points. This should manifest itself as a high-frequency source of noise, the type of which is seen in the hydrophone signals in Figure 14: at frequencies between approximately 100 Hz and 1 kHz for the case with *σ* = 0.53 and at frequencies around 1–2 kHz for *σ* = 1.27, visible especially in the hydrophones HF1 and HF2. The high-frequency broadband noise is likely due in part to very rapid collapse events of the cavity clouds once they reach a region of sufficiently high pressure in the slipstream of the cavitator. A cavitation collapse event is assessed in Section 4.4. The best arrangement of hydrophones in the tests for the investigations with the standard cavitator turned out to be the one in the acoustic chamber, i.e., hydrophones HF1 and HF2.

**Figure 14.** Sound pressure level measurements inside the cavitation tunnel with different hydrophone arrangements. The magenta vertical lines denote the predicted cavitation shedding frequency and its harmonics. The measurement results were reproduced with permission from Rhena Klose (SVA).

#### 4.3.2. Cavitation Dynamics and Shedding Mechanisms

Below in Sections 4.3.3 and 4.3.4, we discuss observations of the sheet-cavitation dynamics and shedding mechanisms. Figure 15 shows a sheet-cavitation shedding cycle of the case with *σ* = 0.53. Figure 16 shows a sheet-cavitation shedding cycle of the case with *σ* = 1.27. Former The results are obtained for the former on the medium grid, and the latter on the fine grid. We study here the cavitation dynamics by visualizing on the left side of a frame the rate of evaporation together with velocity vectors on a cut plane at mid-span, as well as the contour of void fraction value *α* = 0.1. A negative rate of evaporation denotes the rate of condensation. Additionally, on the right side of the frame a perspective view of the cavitation extent is shown by an iso-surface of *α* = 0.1. In each figure, six frames are presented to represent one cavitation shedding cycle. The shedding cycle is taken as the corresponding frequency shown earlier in Table 2. It is notable that the flow features and cavitation dynamics, especially in the case with the lower cavitation number, were observed as highly three-dimensional over the span of the foil. Nevertheless, we focus on a qualitative sense on the flow features of the aforementioned cut plane at mid-span.

**Figure 15.** A sequence of snapshots visualizing one suction side sheet-cavitation shedding cycle. DDES on medium grid, *σ* = 0.53. Left frame depicts a cut plane of the left foil at mid-span, with velocity vectors denoted by arrows (only 25% of the velocity vectors is drawn on the plane for clarity), contour of the void fraction *α* = 0.1 denoted by the grey curves, and the plane is colored by the rate of evaporation. The right frame shows the cavity extent near the left foil by an iso-surface of *α* = 0.1. The frame number denotes the respective stage of the cycle.

**Figure 16.** *Cont.*

**Figure 16.** A sequence of snapshots visualizing one suction side sheet-cavitation shedding cycle. Fine grid, *σ* = 1.27, DDES. Left frame depicts a cut plane of the left foil at mid-span, with velocity vectors denoted by arrows (only 25% of the velocity vectors is drawn on the plane for clarity), contour of the void fraction *α* = 0.1 denoted by the grey curves, and the plane is colored by the rate of evaporation. The right frame shows the cavity extent near the left foil by an iso-surface of *α* = 0.1. The frame number denotes the respective stage of the cycle.

#### 4.3.3. Cavitation Dynamics: *σ* = 0.53

In the first frame of Figure 15 (case *σ* = 0.53) a cavity cloud is shed in the wake of the foil, and an attached sheet cavity has grown from the leading edge of the foil. There is strong evaporation in the layers close to the wall of the sheet cavity and on the upper surface of the cavity. Some condensation takes place at the closure of the sheet cavity. Parts of the shed cavity clouds convect still farther downstream, visible in the right side of the first frame. The frames 1–2 depict the growth of the attached sheet cavitation to super-cavitation, where the sheet cavity momentarily covers the entire suction side of the foil. A re-entrant jet is observed at the closure line of the sheet cavitation in frame 2. Once we reach Frame 3, the re-entrant jet has impinged the sheet cavitation, and continues to a draw liquid phase underneath it. We also notice strengthening condensation at the closure of the sheet cavitation in frames 2–3, as well as strong vorticity as the shedding initiates. By Frame 4, almost all sheet cavity has detached at mid-span, and the re-entrant jet has further intensified. We observe that the impingement does not occur simultaneously over the whole span, but rather exhibits a highly three-dimensional shape and structure with the root and tip parts of the foil still completely covered by the sheet cavitation. As we reach Frame 5, a large cavity cloud is broken off the sheet cavity, which is convected to the wake of the foil. Strong evaporation begins again on the foil as we notice at the cut plane. Finally, once we reach Frame 6 the sheet cavitation at the LE continues to grow, and the broken-off cavity clouds are convected farther in the wake by the external flow over the foil. A new shedding cycle is then initiated at Frame 6. Again we note that while the visually observed shedding cycle was confirmed by the dominant frequencies of the lift and drag coefficients, the two-phase flow over the span is three-dimensional and, for instance, separate and phase-lagging breaking cycles were observed near the root parts of the foil.

#### 4.3.4. Cavitation Dynamics: *σ* = 1.27

Similarly as above, the first frame of Figure 16 (case *σ* = 1.27) depicts the initial stages the attached sheet cavity growth. Cavity clouds are shed in the wake of the foil, and an attached sheet cavity begins to grow from the leading edge of the foil. There is strong evaporation in the layers close to the wall of the sheet cavity, and condensation takes place on the closure of the forming cavity. A detached cavity cloud is convected to the slipstream and vortical flow is associated with it. In frames 2–3, the attached sheet cavitation grows in length along the foil. In Frame 2, we observe evaporation at the center of the attached cavity cloud, whereas condensation takes place around the closure of the sheet cavity. As we reach Frame 3 a re-entrant jet is visible. The jet draws the liquid phase toward the LE, underneath the sheet cavity and splits the evaporation zone of the previous frame in two parts. Frame 4 depicts a partial break-off stage as the re-entrant jet begins to impinge the sheet cavity. Condensation now takes place both in the LE and the TE of the detaching cavity, although evaporation is present at its core. A tiny cavity at the LE of the foil begins to form. Once we reach Frame 5, the impingement of the re-entrant jet is more complete, and we observe that the sheet cavity has broken off on a wide span. The cavity cloud starts to rise away from the foil to be convected by the external flow. Stronger vortical flow forms around the core of the detached cavity. In Frame 6 we observe this cavity cloud being convected farther in the wake, as a new shedding cycle is about to begin. We furthermore observe condensation beginning to dominate over evaporation in the detached cavity. A strong recirculating flow is visible around the detached cavity cloud. Again, the flow is three-dimensional over the span of the foil, and the break-off does not occur simultaneously over the span.

#### *4.4. Cavitation Collapse Events and Pressure Waves*

In Figure 17 we show a set of snapshots for a stage of the shedding cycle at which the cavity undergoes a rapid collapsing stage, and the associated instantaneous pressure distributions. The figure shows the case with the lower cavitation number, *σ* = 0.53. The cavity cloud has shed off at *t*0, and we observe an increase in pressure on the plane at mid-span, just below the cavity. As the collapse stage advances, the pressure values grow larger, and the pressure front starts to propagate in all directions quite rapidly. We note very high values for the *C<sup>p</sup>* due to the major collapse events *t*<sup>0</sup> + 0.8 ms. The maximum value recorded for the investigated plane was ≈3. After this, the pressure wave front grows and intensifies, and then also quickly dissipates. Next, we are left with small regions of negative pressure in the place of the shed cavity clouds at the two latest investigated time instants. Looking at the cavitator body, the pressure wave due to the major collapse at *t*<sup>0</sup> + 0.8 ms reaches it a few tenths of ms later. The cavitator body on the upstream side of the wing shows areas of increased positive pressure, and the downstream side of it in turn has a momentarily lower negative pressure. At the moment of the collapse event, also high pressures are recorded on the suction side of the wing near the trailing edge, as the cavity collapses just above it, at the same time the pressure side of the wing is masked from this effect. The dynamics of the bulk cavitation on the wing, then again, appears to be unaffected by the pressure waves that result from the cavity collapse. The attached sheet cavitation continues to grow on to enter the shedding phase described in Section 4.3.3. Also, at subsequent time instants, the shed-off cavity continues to form the cavitating ring vortex. We assess this further in Section 4.5. The cavitating vortex due to the hub of the cavitator body also appears to be not affected by the collapse event.

**Figure 17.** A sequence of snapshots visualizing a collapse event of the cavity. DDES with medium grid, *σ* = 0.53. Top frames show the left wing from the side and the top, and the bottom frame shows a cut plane at mid-span. The cavitator and the cut plane are colored by *Cp*.

#### *4.5. A Formation of Cavitating Ring Vortex*

In Figure 18 we demonstrate the formation of a cavitating ring vortex. The figure depicts six snapshots at 2 ms intervals of the DDES on a medium grid with *U*∞ = 8 m/s and *σ* = 0.53. On the left side of each frame, the isosurfaces of *α* = 0.1 and Q = 40,000 are shown. Moreover, the *Q*-criterion is colored by helicity H = *<sup>u</sup>i*Ω*i*/(|*u<sup>i</sup>* ||Ω*<sup>i</sup>* |). The helicity denotes the cosine of the angle between the velocity and the vorticity vectors, and tends to ±1 in the cores of free vortices, the sign indicating the direction of swirl. On the right side of each frame, three cut planes are shown at spans of *y*/*s* = [0.25, 0.50, 0.75], which are colored by the rate of evaporation R. Also, on the right side of each frame, the iso-surface of *α* = 0.1 is shown as the transparent orange surface. At the beginning of this series of snapshots (*t* = *t*0), we are approaching the end of a shedding cycle (cf. Figure 15) with cavitation and vortical flow structures being shed in the wake of the foil. A complex interaction with the shedding vapor clouds and liquid-phase flow introduce vast vortical flow structures when the sheet cavity detachment takes place between *t* = *t*<sup>0</sup> and 2 ms. As separate and partially cavitating vortex filaments convect downstream, some of these join together to form a ring-shaped structure, or a toroidal-type vortex between *t* = *t*<sup>0</sup> + 2 and . . . 6 ms. Observing the helicity of these toroids, we note that once it is a true vortex ring or even a piece of such, the vorticity is strong enough for the flow inside these rings to cavitate in most parts clearer and at a higher intensity. Once these structures convect farther downstream (*t* = *t*<sup>0</sup> + 8 . . . 10 ms), dissipation owing to the coarsening grid in the wake diminishes the strength of the toroid-like shape, leading to dissipation of the cavitation as well, despite the used compressive flux limiter. Moreover, as was mentioned above, we observe the strongly three-dimensional cavitation structure in Figure 18. For instance, the rates of evaporation on the three cut planes are vastly different from each other, albeit the sheet cavitation covering the entire foil superficially appears as rather steady at several time instants. Also chord- and spanwise fluctuations can be observed on the iso-surface of the *Q*-criterion.

**Figure 18.** *Cont.*

**Figure 18.** A sequence of snapshots depicting the formation of a partial cavitating ring vortex in the wake of the foil. DDES with medium grid, *σ* = 0.53. On the left in each frame: left foil with iso-surface of void fraction *α* = 0.1 colored as orange, together with an iso-surface of *Q* = 40, 000, which is colored by the helicity H. On the right in each frame: three cut planes at *y*/*s* = [0.25, 0.50, 0.75] which are colored by the rate of evaporation R, with iso-surface of void fraction *α* = 0.1 colored as orange.

#### *4.6. A Comparison of Different Turbulence Closures*

A peculiar modelling issue regarding the cavitation dynamics was observed during the simulations. The case with *σ* = 1.27 was initially begun with the SST model, but the predicted suction side sheet cavitation remained stable throughout the simulations. We observed no natural shedding of the sheet cavitation, and the time histories of the global forces were flat. We carried out the simulations with the SST model until *t* ≈ 0.155 s, after which the DDES approach was activated. The outcome is illustrated in Figure 19 in terms of the time histories of the lift and drag coefficients. It was observed that almost immediately after the DDES approach was activated, the suction side sheet cavitation began its shedding behavior, leading to high-frequency oscillations also in the global force coefficients.

A comparison of the DDES and the SST simulation is shown in Figure 20, made from visualizations of the eddy-viscosity levels near the LE of the foil on a cut plane at mid-span, cavitation extents and vortical flow structures. We observe a significant difference in the eddy-viscosity distributions. The SST simulation predicts several orders of magnitude larger eddy-viscosity levels on the suction side of the foil. This not only suppresses any vortical structures around the foil, but also prevents the natural shedding mechanisms of the sheet cavitation to take place. Instead, the predicted suction side cavitation remained stable as consequently did the lift and drag forces as well, as we noted in Figure 19. By contrast, when using DDES, similar levels of eddy-viscosity are not predicted by the model at the closure of the sheet cavity. We also note considerably more unsteady vortical flow structures with DDES, shed in the wake of the foil. Due to these reasons, the natural shedding of the sheet cavity did occur as was demonstrated in Sections 4.2 and 4.3.

**Figure 19.** Illustration of an effect of DDES. Simulations were carried out using SST *k* − *ω* until *t* ≈ 0.155 s, after which the DDES approach was activated. Medium grid, *U*<sup>∞</sup> = 8 m/s, *σ* = 1.27.

**Figure 20.** Snapshots of simulations with DDES and SST. Medium grid, *σ* = 1.27. Left frame depicts a cut plane of the left foil at mid-span with the contour of void fraction *α* = 0.1 denoted by the grey curves, and the plane is colored by the eddy-viscosity ratio. The right frame shows the cavity extent near the left foil by an iso-surface of the void fraction *α* = 0.1, colored by light blue, and vortical flow structures by an iso-surface of |Ω*<sup>i</sup>* | = 300 1/s as the transparent grey color.

We note here the natural shedding behavior of the case with lower cavitation number was reasonably captured also with the used URANS approach, as discussed in Section 4.2 and seen in Figures 10 and 11. This observation might be because the foil was, in an average sense, at a super-cavitating state where the attached sheet cavitation momentarily covered the entire suction side of the foil. It appears that in this case the shedding mechanism is strong enough to overcome the, perhaps excess, diffusion originating from the eddy-viscosity modelling even on the medium grid.

#### **5. Conclusions**

We studied the ITTC Standard Cavitator numerically in a cavitation tunnel at different cavitation numbers. To our knowledge, our numerical analyses are the first reported ones for the ITTC standard cavitator. Emphasis was put upon the assessment of sheet-cavitation dynamics. A compressible two-phase flow model was used for the flow solution, and two turbulence closures are employed: a two-equation RANS model, SST *k* − *ω*, and a hybrid RANS/LES model, called the delayed-detached eddy simulation (DDES). A homogeneous mixture model was used for the two phases.

We compared the simulations with observations made in the cavitation tunnel experiments. Good agreement was achieved with the mean cavitation extent, in terms of the sheet cavitation on the foils, and the hub vortex cavitation which formed behind the cavitator body. The numerically predicted sheet-cavitation shedding frequencies can also be observed with some hydrophones in the acoustic measurements, which were carried out in the cavitation tunnel. Streak cavitation at the hemispherical head of the cavitator body was predicted as an attached sheet cavitation in the simulations. The numerical grid resolution did not permit the capturing of the very fine tip vortex cavitation that was observed in the experiments. Also, the prediction of flows with bubble cavitation is difficult with a homogeneous two-phase mixture model.

Detailed analysis of the cavitation shedding mechanism confirmed that the dynamics of the sheet cavitation are dictated by the re-entrant jet. Once the jet impinges the sheet cavitation, a detachment takes place together with complex interaction between the re-entrant jet, vortical flow structures and mass transfer near the phase boundary. We observed that the break-off cycle was relatively periodic in both investigated cases with an approximately constant shedding frequency. The Strouhal numbers lay within the usual ranges reported in the literature for sheet-cavitation shedding. Furthermore, we noted that a strong vortical flow formed around the detached cavity. In the case of the lower cavitation number, we demonstrated that the vortical flow structures convected to the slipstream of the foils, where they developed striking cavitating toroidal vortices. An intricate interplay with the shed vapor clouds and the liquid-phase flow introduced vast vortical flow structures when the sheet cavity detachment took place. We showed that once the toroidal vortex ring was formed, the vorticity was in

most parts strong enough for the flow inside the vortex rings to cavitate. Furthermore, a collapse event of a cavity cloud in the slipstream of the foils gave rise to a considerable pressure peak and subsequent wave-like fronts.

DDES and simulations with the URANS model revealed a difference of approximately 5 Hz in the sheet-cavitation shedding for the case with the lower cavitation number. The shedding behavior was reasonably captured also with a two-equation turbulence model. This observation might be because the foil was, in an average sense, at a super-cavitating state where the attached sheet cavitation momentarily covered the entire suction side of the foil. It appears that in this case the shedding mechanism is strong enough to overcome the, perhaps excess, diffusion originating from the eddy-viscosity modelling even on the rather coarse grids. For the case with the higher cavitation number, the two-equation turbulence model predicted a steady attached sheet cavitation, whereas the hybrid RANS/LES model predicted a highly transient sheet-cavitation shedding being on a par with the experiments. We attributed this difference to the eddy-viscosity modelling. In this case, detailed vortical structures around the foil and the shedding mechanisms of the sheet cavitation were inhibited due to the unsteady RANS modelling approach.

Recently we have investigated the Reynolds number sensitiveness of different cavitation types observed on a ship's propeller [31]. For the case of static hydrofoils exhibiting similar phenomena as shown here, this remains to be assessed. Moreover, the ability of an asymptotic mass-transfer model coupled to the compressible flow solver to predict the shock wave dynamics, which can have an effect on the upstream cavitation dynamics [1], remains an open issue. One alternative is to include the Rayleigh–Plesset equation, which takes the bubble dynamics better into account, to the flow model for a more complete physical description of the two-phase flow. Other options to achieve better results comprise more holistic modelling approaches that can better account for the phase change physics, such as inhomogeneous or multi-scale flow models. An extended model of the former type is under development by the authors, and we have obtained promising initial results of propeller cavitation using such an approach.

**Author Contributions:** Conceptualization, V.M.V., T.S. (Timo Siikonen), A.S.-C. and T.S. (Tuomas Sipilä); methodology, V.M.V. and T.S. (Tuomas Sipilä); software, V.M.V. and T.S. (Timo Siikonen); validation, V.M.V.; formal analysis, V.M.V.; investigation, V.M.V.; resources, V.M.V. and T.S. (Tuomas Sipilä); data curation, V.M.V. and T.S. (Timo Siikonen); writing—original draft preparation, V.M.V.; writing—review and editing, T.S. (Timo Siikonen), A.S.-C. and T.S. (Tuomas Sipilä); visualization, V.M.V.; supervision, T.S. (Timo Siikonen), A.S.-C. and T.S. (Tuomas Sipilä); project administration, T.S. (Tuomas Sipilä); funding acquisition, T.S. (Tuomas Sipilä). All authors have read and agreed to the published version of the manuscript.

**Funding:** V.M.V. would like to express his gratitude for the support granted by the Finnish Funding Agency for Innovation (TEKES) and the German Federal Ministry for Economic Affairs and Energy within the project "PropNoise". Without the funding the research project could not have been realized. The APC was funded by VTT.

**Acknowledgments:** The authors would like to express their gratitude to Rhena Klose and Lars Lübke from SVA Potsdam GmbH for providing the experimental results and photographs.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


RANS Reynolds-averaged Navier–Stokes

SST Shear stress transport

TE Trailing edge

URANS Unsteady RANS

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Cavity Formation during Asymmetric Water Entry of Rigid Bodies**

**Riccardo Panciroli 1,\* and Giangiacomo Minak <sup>2</sup>**


**Abstract:** This work numerically evaluates the role of advancing velocity on the water entry of rigid wedges, highlighting its influence on the development of underpressure at the fluid–structure interface, which can eventually lead to fluid detachment or cavity formation, depending on the geometry. A coupled FEM–SPH numerical model is implemented within LS-DYNA, and three types of asymmetric impacts are treated: (I) symmetric wedges with horizontal velocity component, (II) asymmetric wedges with a pure vertical velocity component, and (III) asymmetric wedges with a horizontal velocity component. Particular attention is given to the evolution of the pressure at the fluid–structure interface and the onset of fluid detachment at the wedge tip and their effect on the rigid body dynamics. Results concerning the tilting moment generated during the water entry are presented, varying entry depth, asymmetry, and entry velocity. The presented results are important for the evaluation of the stability of the body during asymmetric slamming events.

**Keywords:** slamming; fluid-structure interaction; fluid detachment; cavitation

**Citation:** Panciroli, R.; Minak, G. Cavity Formation during Asymmetric Water Entry of Rigid Bodies. *Appl. Sci.* **2021**, *11*, 2029. https://doi.org/ 10.3390/app11052029

Academic Editor: Florent Ravelet

Received: 30 January2021 Accepted: 22 February 2021 Published: 25 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The phenomenon of a structure impacting on the water generally implies large forces due to the massive fluid volume displaced in a small time. The duration of the slamming event is in the order of milliseconds. These loads might damage the structure or, because of their short duration, excite a dynamic response of the local structure of the hull and cause vibrations. Non-symmetric impacts might further introduce a tilting moment and influence the impact dynamics and structural stability.

The first analytical solution to solve the impact dynamics of rigid bodies entering the water was presented by Von Karman [1], who developed a formula capable of predicting the maximum force acting on a rigid body entering the water, to make a stress analysis on the members connecting the fuselage with the floats of a seaplane. Many analytical methods have been proposed to extend Wagner's method to different shapes (e.g., [2–7]) and most of them are very effective in predicting the water entry of simple-shaped structures impacting the surface with pure vertical velocity. Some of these solutions are even capable of accounting for oblique impacts (e.g., [8–12]). It is reported that particular conditions of entry velocity, deadrise angle, and tilt angle might lead the fluid to detach at the wedge apex, introducing difficulties in evaluating the pressure at the interface by analytical formulations. A typical water entry event that largely deviates from a symmetric case is ditching [13–15], where fluid detachment and cavitation arise if the advancing velocity is sufficiently large. Another area of interest for asymmetric impacts is the analysis of the response of planing hulls during the maneuvering operation since, depending on the conditions, restoring or capsizing moments can take place [16]. Xu [17] defined two types of asymmetric impact. Type A flow is the one when there is small asymmetry and the flow moves outward along the contour on both sides of the vertex. Type B flow occurs when there is large asymmetry and the flow detaches from the body contour at the vertex on

one side. Chekin [18] concluded that there was only one unique combination of wedge angle and impact angle at which no separation of flow from the vertex would occur. For a given wedge shape and impact direction, any other impact angle would force separation. Defining *U* as the horizontal velocity and *V* as the vertical velocity, the advancing ratio *U*/*V* at which the flow separation appears is less for bodies of larger deadrise angles. For small asymmetry impacts, the cavity flow during the water-entry is limited to a very small region. Furthermore, the flow that separates from the apex quickly re-attaches to the wedge. A symmetric body impacting with horizontal velocity will produce a flow similar to asymmetric impact with only vertical velocity when a tilting motion is not permitted. Judge et al. [8] performed experiments on wedges where asymmetry and horizontal impact velocity were present and compared the results with an analytical solution, showing good agreement for low angles of asymmetry and small advancing ratios.

Xu et al. [17] observed that the pressure near the tip of asymmetric wedges might attain negative values with respect to the gauge pressure, and hence be lower than the ambient pressure. They also observed that the peak pressure increases with the horizontal velocity, and with it, an increase in the wave elevation on the upstream side and a decrease on the downstream one. However, the simulations are based on the assumption of attached flow. In this work, we also consider the eventual detachment at the wedge tip, similarly to that described in [19] but utilizing a different numerical scheme. Semenov [20] studied the effect of the horizontal component of the entry velocity for various wedge orientations. Results show that certain configurations of the impact might induce a negative pressure along the whole wetted length of the downstream side of the wedge, leading to flow separation at the wedge apex. Semenov and coauthors [21–23] derived an analytical solution for the asymmetric/oblique water entry of a wedge which does not assume flow separation from the wedge vertex. It was shown that at some wedge orientations, the total force acting on the wedge side with the larger deadrise angle may become zero, a condition corresponding to flow separation.

Experimental evidence of the onset of cavity formation at high asymmetry ratios has been provided by Shams et al. [24], who studied asymmetric water impacts experimentally through particle image velocimetry.

In this manuscript, we initially define the problem of asymmetric water impacts in Section 2 and we detail the numerical model utilized for the analysis in Section 3. Asymmetric water entry events are then divided into categories and studied separately in Sections 4–6, to finally draw some final remarks in Section 7.

#### **2. Problem Statement: Water Entry of Asymmetric Wedges**

In the most general case, wedges enter the water with combined vertical and horizontal velocity components, whose ratio defines the advancing ratio *ǫ*, and are not symmetric due to an initial tilt angle *γ*. As mentioned in the introduction, this kind of impact might introduce ventilation into the fluid flow [8] due to fluid detachment at the wedge tip. A sketch of the problem is shown in Figure 1: a wedge with nominal deadrise angle *β* enters the water at the velocity *V*<sup>0</sup> and is rotated by a tilt angle *γ* with respect to the water surface.

**Figure 1.** Sketch of the asymmetric impact of a wedge.

In the following, the analysis will always refer to wedges designed as symmetric structures with nominal deadrise angle *β* but rotated by a tilt angle *γ*, hence presenting different deadrise angles on the two sides (*β*<sup>1</sup> = *β* − *γ* and *β*<sup>2</sup> = *β* + *γ*). We also define the advancing ratio *ǫ* as the ratio between the horizontal velocity *U*<sup>0</sup> and the vertical velocity *V*0. Utilizing such definitions, we here classify the water entry event into three main categories:


In the case of vertical impact, fluid detachment is prone to happen only if the condition *β* + *γ* ≥ *π*/2 is met. In the presence of a horizontal velocity component, it is reported in the literature [20] that the flow separation from the wedge vertex may occur when *<sup>π</sup>* <sup>−</sup> *<sup>γ</sup>*inf <sup>&</sup>lt; *<sup>β</sup>L*, being *<sup>β</sup><sup>L</sup>* is the leeward-side deadrise angle, and *<sup>γ</sup>*inf the water-entry angle, as in these cases, the pressure near the tip might become negative [25]. Although this phenomenon is well known, it is usually neglected due to the difficulties in treating the fluid detachment. This work instead focuses on those cases where the fluid detaches at the wedge vertex.

#### **3. Coupled FEM/SPH Numerical Model**

In this study, we assume that cavitation at the wedge vertex does not occur and flow separation depends on the flow characteristics enabling ventilation. Gravity is always neglected, while the entry velocity and tilt motion either follow the actual impact dynamics or is set to a constant value, for in the latter case the solution is self-similar in time, hence the results might be normalized by the entry depth. In the case of full fluid–structure interaction simulations (hence modelling the whole impact dynamics) the solution is not self-similar in time not only because of the role of acceleration but also because the tilt motion modifies the deadrise angle. As a result, the jet of water piling up over the counterclock side thickens, while in the opposite side, the water is pushed by the rotating wedge and the free surface profile near the wedge deforms significantly.

The 2D numerical model , presented in Figure 2, is based on the SPH solutor available within the commercial software LS-DYNA. In all the numerical solutions, the fluid is modelled by SPH particles with a constant diameter of 0.25 mm covering a region 0.8 m wide and 0.3 m deep. Non-reflecting boundaries define the bounding box. The fluid is modelled through an equation of state following the Gruneisen model and the fluid particles are solved utilizing the SPH-enhanced fluid formulation (FORM = 16). Both wedge's sides are 0.3 m long and are modelled as a rigid shell composed by 100 elements on each side, leading to an average of four SPH particles that make contact with a single element in the wet region. Such a strategy allows to smooth out the marked pressure fluctuations otherwise attained at the fluid–structure interface. Deadrise angle, initial rotation, and entry velocities are parametrically varied, while the wedge mass is 100 kg per unit width.

**Figure 2.** Two-dimensional numerical model with highlighted boundary and initial conditions.

The numerical results will be compared against established analytical solutions. It has already been shown that the numerical scheme utilized within this work is in line with the analytical predictions provided by Wagner's model for the case of pure vertical velocity [26,27] which, in turn, was shown to provide results in line with experimental results conducted on rigid wedges [28–31]. In Wagner's expanding plate model, which has been formulated for symmetric impacts only, the pressure distribution along the wedge edge is given by:

$$\frac{p}{\rho} = \dot{V}\sqrt{r^2 - \mathbf{x}^2} + \frac{V \ r \ \dot{r}}{\sqrt{r^2(t) - \mathbf{x}^2}} - \frac{1}{2}\frac{V^2 \mathbf{x}^2}{r^2 - \mathbf{x}^2} \tag{1}$$

while the Modified Logvinovich Model (MLM), before flow separation, reads [32]

$$\frac{p}{\rho} = \frac{1}{2}V^2 \left[ \frac{\dot{r}}{V} \frac{r(t)}{\sqrt{r^2(t) - x^2}} - \cos^2 \beta \frac{r^2(t)}{cr^2(t) - x^2} - \sin^2 \beta \right] \tag{2}$$

where 0 < *x* < *r* defines the wet portion of the edge. Within this scheme, *x* is an abscissa oriented horizontally, i.e., following the direction of the free surface, and *r* is a function of the entry depth *ξ* and equals *<sup>π</sup>* 2 *ξ* tan *β* . In the case of constant entry velocity, the horizontal projection of the wet length reads *r* = *<sup>π</sup> V t* 2 tan *β* . MLM can be also utilized to study asymmetric impacts. However, a simplified method considering an equivalent deadrise angle equal to the average of the left- and the right-side deadrise angles is sufficiently accurate to predict the overall impact dynamics [32]. Within this numerical scheme, the distribution of the non-dimensional hydrodynamic pressure *p*¯, which is given by <sup>2</sup>*<sup>p</sup>* tan *<sup>β</sup> ρV*<sup>2</sup> , over the horizontal abscissa *x* is given by

$$p(\mathbf{x}) = \frac{\pi}{\sqrt{1 - \left(\frac{\mathbf{x}}{r}\right)^2}} - \tan\beta \frac{\left(\frac{\mathbf{x}}{r}\right)^2}{1 - \left(\frac{\mathbf{x}}{r}\right)^2} \tag{3}$$

The normalized pressure distribution for varying deadrise angle is shown in Figure 3.

**Figure 3.** Evolution of the normalized pressure along the wetted edge for variable deadrise angles.

Following Wagner's solution, the vertical force (per unit width) of the wedge can then be estimated as follows:

$$F(t) = 2 \int\_0^1 \tilde{p} d\mathbf{x} \, \frac{\rho V^2}{2 \tan \beta} r \cos \beta \quad \text{or} \qquad F(t) = F \frac{\pi}{2} \frac{\rho V^3}{\tan \beta} \frac{t \cos \beta}{\tan \beta} \tag{4}$$

This shows that the ratio between the vertical force and the entry depth, which in the case of constant velocity equals *Vt*, is

$$\frac{F(t)}{V \ t} = F \frac{\pi \rho V^2}{2 \tan \beta} \frac{\cos \beta}{\tan \beta} \tag{5}$$

which is time-independent, being a function of the deadrise angle only. Here, vertical force, horizontal force, and tilting moment will be normalized as

$$\bar{F}\_{\mathbf{x},\mathbf{y}} = F\_{\mathbf{x},\mathbf{y}} \frac{2\tan^2\beta}{\cos(\beta)\pi\rho V\_{\mathbf{y}}^3 t} \quad \text{and} \quad \bar{M} = M \frac{2\tan^2\beta}{\cos(\beta)\pi\rho V\_{\mathbf{y}}^2 V\_{\mathbf{x}} t^2}. \tag{6}$$

#### **4. Symmetric Wedges with Horizontal Velocity Component—***γ* = **0 and** *ǫ* 6= **0**

In this section, we investigate the effect of a horizontal velocity component on the pressure distribution over a rigid wedge. At first, we concentrate on rigid wedges with constant velocity and null tilt motion, since such approximation leads to a self-similar solution and the impact time can be taken out of the analysis.

As a general result, the horizontal velocity component alters the fluid motion beneath the wedge up to the point a cavity is formed. Indeed, a threshold value of it exists for the onset of cavity formation. Figure 4 shows three images of the cavity predicted by the numerical solutions on a wedge with *β* = 25◦ and varying advancing ratios *ǫ*. Notably, fluid detachment is attained as *ǫ* exceeds 2 and the void region increases with the advancing ratio.

**Figure 4.** Fluid detachment at the vertex of a *β* = 25◦ wedge for varying advancing ratio *ǫ* = <sup>6</sup> 3 , 7 3 and <sup>8</sup> 3 (left to right).

Indeed, pressure drops to zero in the dry region, the extension of which increases with the advancing ratio *ǫ*. Figure 5 details the predicted pressure distribution over a symmetric wedge with deadrise angle *β* = 30◦ entering the water with varying advancing ratios *ǫ*. These examples report wedges running at a constant speed and the tilt motion is inhibited. Each graph reports the normalized pressure over the normalized wet portion of the wedge at various instants of the impact. Notably, the solution remains self-similar for all the advancing ratios considered, as the normalized pressure and void region are not influenced by the entry depth. Interestingly, the leeward wedge side remains completely dry as *ǫ* exceeds 2.

**Figure 5.** Normalized pressure versus wet portion of a wedge (*β* = 30◦ ) at various instants of the impact for varying *ǫ*.

The self-similarity of the results allows for the effect of the advancing ratio on the pressure distribution along the body to be compared. Figure 6 reports results similar to the ones presented in Figure 5, but stacked on a single graph to ease their comparison. Results about *β* = 20◦ and *β* = 30◦ wedges with varying advancing ratios are presented. In the case of pure vertical impact, such deadrise angles experience a hydrodynamic load that is almost constant over the entire wet surface, as also predicted by the Wagner's solution reported in Figure 3. However, when a horizontal velocity component is introduced, the pressure on the right side (that is, the windward side) is found to increase almost linearly moving from the vertex to the free edge, while it remains constant on the leeward side, except for the region experiencing fluid detachment. The graph reports a single solution for each case due to the self-similarity of the results. In the case of *β* = 30◦ , fluid detachment initiates when *ǫ* becomes greater than 1 and the magnitude of the void region with respect to the overall wet length increases with the advancing ratio, to eventually remain fully dry when *ǫ* is greater than 2. Notably, for both the deadrise angles considered, the cavity incipient condition is met as soon as *β* + arctan *ǫ* ≈ 1.3, while the cavity is completely formed as soon as *β* + arctan *ǫ* ≈ *π*/2, a result which matches with the analytical predictions [25], validating the capability of the numerical model to correctly capture the cavity formation event. The pressure on the windward side is found to increase with the horizontal component while it decreases on the opposite side. At the higher advancing ratio considered in this study (*ǫ* = 2.5), the left side never makes contact with the fluid. Figure 6 also highlights, with particular reference to the *β* = 30 case, how the pile-up on the windward side of the wedge increases with the advancing ratio, while it decreases on the leeward side. These results should be ascribed to the increasing horizontal component of the velocity.

**Figure 6.** Normalized pressure versus normalized wetted surface on a symmetric wedge with *β* = 30◦ (**left**) and *β* = 20◦ (**right**), entering the water with different advancing ratios *ǫ* = *Vx*/*Vy*.

Having imposed a constant entry velocity further allows normalizing the impact forces and the tilting moment, whose magnitude remain constant during the whole impact event. The normalized results for all cases considered are collected in Figure 7. The normalized horizontal force *F<sup>x</sup>* is found to linearly increase with the advancing ratio, as the horizontal velocity component increases with it—while both the vertical force *F<sup>y</sup>* and the momentum *M<sup>z</sup>* remain constant (hence proportional to the vertical force only)—to eventually diverge as *ǫ* exceeds 2, which is the full cavity generation condition.

**Figure 7.** Normalized horizontal force, vertical force, and momentum versus advancing ratio *ǫ*.

#### **5. Asymmetric Wedges with Pure Vertical Velocity—***γ* 6= **0 and** *ǫ* = **0**

This section firstly proposes to adapt simple analytical formulations to investigate the impact dynamics of asymmetric wedges. Later, the hydrodynamic pressure evaluated numerically will be compared against Wagner's solution to assess its validity in asymmetric water impacts.

#### *5.1. Impact Dynamics of Asymmetric Wedges*

The impact dynamics of the asymmetric wedge impacting with pure vertical velocity will be here approximated through Wagner's formula, by adjusting the added mass term to account for the different deadrise angles, hence the wet length, of the two sides of the wedge. The mass of the expanding plate in Wagner's approach is here approximated by

$$m = \rho \frac{\pi}{2} \left[ \frac{1}{2} \left( \frac{\mathfrak{J}}{\tan(\beta\_0 + \gamma)} + \frac{\mathfrak{J}}{\tan(\beta\_0 - \gamma)} \right) \right]^2 \tag{7}$$

And the impact dynamics is estimated utilizing Equation (4) by substituting the term tan(*β*) with

$$\tan(\beta) = 2 \frac{\tan(\beta\_0 + \gamma) \cdot \tan(\beta\_0 - \gamma)}{\tan(\beta\_0 + \gamma) + \tan(\beta\_0 - \gamma)}\tag{8}$$

The comparison between analytical predictions and the numerical results shows a fairly good match during the initial stage of the water entry, to slightly diverge after the peak acceleration is reached. Such a difference should be ascribed to the tilt motion which modifies the deadrise angle over time, an effect that is not accounted for in the analytical model. As an example, Figures 8–11 show the comparison between numerical results and analytical predictions of the impact dynamics of a *β* = 25◦ wedge varying tilt angle and entry velocity.

**Figure 8.** Impact dynamics of an asymmetric wedge entering the water at 6 m/s. Deadrise angle *β* = 25◦ , tilt angle 5◦ .

**Figure 9.** Impact dynamics of an asymmetric wedge entering the water at 8 m/s. Deadrise angle *β* = 25◦ , tilt angle 5◦ .

**Figure 10.** Impact dynamics of an asymmetric wedge entering the water at 4 m/s. Deadrise angle *β* = 25◦ , tilt angle 15◦ .

**Figure 11.** Impact dynamics of an asymmetric wedge entering the water at 6 m/s. Deadrise angle *β* = 25◦ , tilt angle 15◦ .

#### *5.2. Hydrodynamic Pressure*

It has been found that the impact dynamics during the initial stages of the water entry process can be effectively predicted by a simple modification of the expanding plate model. To evaluate the effect of the tilt angle on the pressure distribution, we now concentrate on the hydrodynamic pressure at the interface and compare the results against Wagner's predictions.

Figure 12 shows the effect of the tilt angle on the pressure distribution for the case of a *β* = 30◦ wedge with constant entry velocity and null tilt motion. In such a condition, the solution is found to remain self-similar in time and the trend of the hydrodynamic pressure follows Wagner's prediction on each side of the wedge. The pileup coefficient [33] is not influenced by the initial tilt angle and it remains approximately constant to *π*/2, matching with the analytical value predicted for symmetric wedges.

**Figure 12.** Effect of the initial tilt angle on the pressure distribution on a *β* = 30◦ wedge entering the water with pure vertical velocity for varying initial tilt angle.

Indeed, if the velocity is not prescribed a priori, the simulation loses the self-similarity and the solution changes over time. Figure 13 reports, as an example, the time evolution of the pressure at the fluid–structure interface of a wedge with deadrise angle *β* = 25◦ entering the water with an initial tilt angle of 5◦ and initial velocity of 6 m/s. The hydrodynamic pressure compares well during the initial stage of the impact but soon diverges in both sides, to eventually switch the side experiencing the higher pressure. As already mentioned before, such a difference should be ascribed to the tilt motion, which is not taken into account in the analytical model. Please note that the sudden pressure jump at the keel, the magnitude of which increases with *ǫ*, is a graphical artefact introduced by the normalization of the results.

**Figure 13.** Asymmetric wedge entering the water at 6 m/s with pure vertical velocity. Deadrise angle *β* = 25◦ , tilt angle 5◦ . Pressure distribution at 4, 8, 12, and 16 ms from the impact.

The presented results show that, on one side, an eventual asymmetry during the water impact does not have a remarkable effect on the impact dynamics and on the maximum impact force, since the peak force is attained in the early stage of the impact, prior to a substantial tilt motion. On the other side, the tilt motion is found to largely alter the hydrodynamic pressure distribution at the larger entry depths.

The next section presents some numerical solutions for asymmetric wedges entering the water with combined vertical and horizontal velocity to investigate the combined influence of tilt angle and horizontal component.

#### **6. Asymmetric Wedges with Horizontal Velocity Component—***ǫ***,** *γ* 6= **0**

In this section, asymmetric wedges with combined vertical and horizontal impact velocities are studied, which represents a combination of the impact conditions studied in Sections 4 and 5. The analysis considers only the cases where the advancing ratio is in the same direction of the tilt angle (left direction in Figure 1, and counter-clock wise rotation). Advancing ratio, deadrise angle, and tilt angle have been parametrically varied in the numerical analysis, but only some reference cases are reported here for brevity.

As an example, two cases representative of incipient cavity formation are presented: Figure 14 compares the numerical results and analytical predictions of the pressure over a 25◦ wedge with an initial tilt angle of 15◦ entering the water with vertical velocity *v<sup>y</sup>* = 4 m/s and horizontal velocity *v<sup>x</sup>* = 2 m/s. Figure 15 shows a wedge with a deadrise angle of 25◦ and tilt angle of 5◦ entering the water with vertical velocity *v<sup>y</sup>* = 4 m/s and horizontal velocity *v<sup>x</sup>* = 4 m/s (*ǫ*=1).

**Figure 14.** Asymmetric wedge entering the water at 4 m/s (y) + 2 m/s (x). Deadrise angle *β* = 25◦ , tilt angle 15◦ .

**Figure 15.** Asymmetric wedge entering the water at 4 m/s (y) + 4 m/s (x). Deadrise angle *β* = 25◦ , tilt angle 5◦ .

In the case shown in Figure 15, the high horizontal velocity component induces the fluid to detach from the wedge tip, as visible in Figure 16, reporting the fluid velocity. Since air is not considered in the numerical model, the numerical solution might not represent the real behaviour of the fluid. Nevertheless, the results show that even if there is a region where the pressure drops to zero, the maximum pressure compared well with Wagner's model, even if the overall shape changed: the maximum pressure on the leeward face is now located where the fluid reattach to the wedge rather than in the fluid jet.

**Figure 16.** Velocity field in the fluid. Asymmetric wedge entering the water with combined horizontal (4 m/s) and vertical (4 m/s) velocity. Deadrise angle *β* = 25◦ , tilt angle 5◦ .

Overall, the results show that the combination of leeward deadrise angle and advancing ratio resulting in incipient cavity formation is *β* + arctan *ǫ* > 1.2, and full cavity formation is attained when *β* + arctan *ǫ* > *π*/2. These results are in line with the ones obtained for the symmetric wedges with horizontal velocity component presented in Section 4.

#### **7. Conclusions**

This work demonstrates that the proposed numerical solution scheme based on a coupled FEM–SPH approach allows for prediction of the hydrodynamic pressure during asymmetric water entry of rigid bodies and also in the case of cavity formation at the wedge vertex. The following three cases were considered: symmetric wedges entering the water with a velocity vector with components normal and tangential to the fluid surface, asymmetric wedges falling vertically, and, finally, the combination of asymmetry and a generic velocity vector. In the first scenario, the results show that the horizontal component of the velocity influences the fluid motion beneath the wedge and that a threshold value of the ratio between the horizontal and the vertical velocity components exists, leading to the onset of the cavity formation. This threshold depends on the deadrise angle. Cavity formation has an influence on the impact forces and on the tilting moment that increases significantly for advancing ratios above the threshold. On the other side, in pure vertical entry events, it is found that an asymmetry in terms of the initial tilt angle does not affect the maximum hydrodynamic load. This is very important because, as a consequence, the impact dynamics can be effectively predicted by simply introducing the tilt angle in the analytical formulations developed for symmetric impacts. However, the tilt motion alters the hydrodynamic pressure distribution at the larger entry depths. Finally, the combination of the two conditions shows that the maximum pressure is comparable to Wagner's solution, even if the pressure distribution is very different and presents zero pressure zones. As in the symmetric case, the cavity formation presents a threshold, which, in this case, is a function of the deadrise angle and the advancing ratio. As a further step, the proposed numerical scheme can be coupled with FEM to determine the effect of asymmetry on the hydro-elastic behaviour of flexible structures.

**Author Contributions:** R.P. contributed to carrying out the simulations and the data analysis. He is the main writer of the paper. G.M. contributed to the redaction of the manuscript and was responsible for the supervision. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not Applicable.

**Informed Consent Statement:** Not Applicable.

**Data Availability Statement:** The data presented in this study could be available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Analysis of Bulb Turbine Hydrofoil Cavitation**

**Andrej Podnar, Marko Hoˇcevar, Lovrenc Novak and Matevž Dular \***

Faculty of Mechanical Engineering, University of Ljubljana, Aškerˇceva 6, 1000 Ljubljana, Slovenia; andrej.podnar@gmail.com (A.P.); marko.hocevar@fs.uni-lj.si (M.H.); lovrenc.novak@fs.uni-lj.si (L.N.) **\*** Correspondence: matevz.dular@fs.uni-lj.si; Tel.: +386-1-477-1453

#### **Featured Application: Bulb water turbine blades with improved cavitation properties.**

**Abstract:** The influence of a bulb runner blade hydrofoil shape on flow characteristics around the blade was studied. Experimental work was performed on a bulb turbine measuring station and a single hydrofoil in a cavitating tunnel. In the cavitation tunnel, flow visualization was performed on the hydrofoil's suction side. Cavitation structures were observed for several cavitation numbers. Cavitation was less intense on the modified hydrofoil than on the original hydrofoil, delaying the cavitation onset by several tenths in cavitation number. The results of the visualization in the cavitation tunnel show that modifying the existing hydrofoil design parameters played a key role in reducing the cavitation inception and development, as well as the size of the cavitation structures. A regression model was produced for cavitation cloud length. The results of the regression model show that cavitation length is dependent on Reynolds's number and the cavitation number. The coefficients of determination for both the existing and modified hydrofoils were reasonably high, with R<sup>2</sup> values above 0.95. The results of the cavitation length regression model also confirm that the modified hydrofoil exhibits improved the cavitation properties.

**Keywords:** hydrofoil; bulb turbine; bulb turbine runner; cavitation; flow visualization; cavitation tunnel; regression model

#### **1. Introduction**

Historically, hydropower was used as a source of baseload generation for electric power, but now it is increasingly used for grid balancing. Among dispatchable renewable energy sources, hydropower plants are the most often used [1], achieving up to 15 s response times [2] for 90% of the total power response. The future potential of hydropower is shifting towards grid regulation tasks. Including nondispatchable renewable energy sources, e.g., inter alia wind and photovoltaics, and market deregulation, creates the requirement for rapid power generation fluctuations and power increase in the afternoon [3–5]. Any new hydraulic turbine design should therefore enable this type of operation and reach high efficiency and stable operation outside its optimal operating interval. Besides stable continuous turbine operation away from the best efficiency point, new designs will also need to function sufficiently in transient operation. Bulb turbines face challenges for such operation; designers will have to overcome the detrimental effects that flow fluctuations and cavitation have on turbine operation, while still ensuring adequate performance over a broad interval of operating points and in transient operation [3].

At the best efficiency point, an optimal combination of the runner frequency of the rotation, head, and discharge is available. At such an operating point, flow separation is usually low, and no cavitation is observed. For operation away from the optimum efficiency point, bulb turbines, as compared e.g., to Francis turbines [4], feature double regulation and may operate in a wide interval of operating conditions with reasonably high efficiency. Although all hydraulic turbine types are constantly improved, the gains in efficiency are in the per mille range and the cavitation in hydraulic turbines remains an

**Citation:** Podnar, A.; Hoˇcevar, M.; Novak, L.; Dular, M. Analysis of Bulb Turbine Hydrofoil Cavitation. *Appl. Sci.* **2021**, *11*, 2639. https://doi.org/ 10.3390/app11062639

Academic Editor: Florent Ravelet

Received: 12 February 2021 Accepted: 10 March 2021 Published: 16 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

issue. The cavitation in the hydraulic turbines [4,6,7] decreases efficiency and increases vibrations and erosion, thereby reducing the service life.

The procedures optimizing the performance of turbine hydraulics are most often based on the use of so-called inverse design methods [8] and computational fluid dynamics (CFD) [9].

Inverse design methods have been in use for decades, with the first 3D methods emerging around 30 years ago. These methods are based on the a priori specification of the pressure or blade loading distribution. Such distribution along the runner blade in the radial and meridional directions facilitates the calculation of the runner blade's angle. This reduces the computational CFD work required when insufficient information on the runner geometry is known in advance. Zangeneh [10] proposed a three-dimensional design method for radial and mixed turbines for the compressible case. Today his method is still often used for designing radial turbine machinery. Zhu et al. performed inverse design optimization of the reversible pump-turbine [11], while also considering cavitation characteristics. Leguizamón and Avellan [12] very recently provided an open-source implementation of an inverse Francis turbine design, together with a validation study confirming the approach's appropriateness. Cavitation is not often studied using inverse design methods; among the available studies, Daneshkah and Zangeneh [13] linked the blade loading and blade stacking to the efficiency and cavitation characteristics of a Francis runner, providing a valuable physical understanding and useful design guidelines.

The 3D design of the runners and the guide apparatus may be determined by performing consecutive CFD calculations, while constantly varying 3D geometry and design parameters [9,14]. Unfortunately, CFD alone is unable to establish an optimized 3D design. A 3D design study using CFD was performed by Xue et al. [15], considering the efficiency, stability of design, and cavitation of a pump turbine during pumping mode. A comprehensive hydraulic and parametric analysis of a Francis turbine runner was performed by Ma et al. [16]. They [16] performed a multipoint and multiobjective optimization design procedure, achieving efficiency improvements up to 0.91%. The optimization procedure included the 3D inverse design method, CFD analysis, and multipoint and multiobjective optimization technology.

Runner blades are central turbine elements, where energy is transferred from the fluid to the generator shaft. Blade sizes and angles, including hydrofoil selection, contribute to the blade design. The blade design is often a trade-off between efficiency, cavitation development, and cavitation erosion. A pressure difference between the suction and pressure sides is limited by cavitation, while at the same time the high pressure difference enables energy transfer from water to the generator. If the turbine is at full load and the water level is low in the reservoir, the pressure on the suction side may decrease and locally drop below the vapor pressure. Cavitation on the runner blades may appear in several locations. These are the turbine blade leading edge, the tip, the pivot, the gap between the blade tip and the discharge ring, the blade suction side close to the trailing edge, and the cavitation near the draft tube cone [7,17]. For the abovementioned cavitation types, the root causes are suboptimal runner blade design and hydrofoil selection for the selected operating conditions. Cavitation is extensively studied in hydraulic turbine machinery; among these, few studies have employed the visualization method [6,7].

On bulb turbines' blades, the initial cavitation most often appears at the location with the highest local speed i.e., at the blade tip [18]. We show an experimental analysis of a modification to the hydrofoil design of bulb turbine blades so as to resist cavitation during continuous operating conditions at high flows. The paper shows the cavitation characteristics of two bulb turbine hydrofoils. The improvement from the existing to modified hydrofoil is discussed based on visualization analysis in the cavitation channel. Furthermore, we show a comparison of the modified runner blades with the unmodified blades in terms of the turbine output. The paper shows marked improvement in cavitation occurrence between existing and modified hydrofoils, and a confirmed regression model of the void fraction.

#### **2. Hydrofoil**

Bulb turbine (Figure 1) blade hydrofoils are optimized for a mix of parameters, among which are maximum hydraulic efficiency, high output power, minimum turbine vibrations, high strength, reduced fatigue, and reduced cavitation phenomena. Reduced cavitation enables operators to reduce the service intervals required for cavitation damage repair and welding costs, as well as loss of money for nonoperation.

**Figure 1.** Bulb turbine runner.

In this study, we focus on a range of volumetric flows, which appear on a bulb turbine at a 25 ◦ blade angle (Figure 2). Here cavitation frequently occurs at approximately 15% higher volumetric flows than at the best efficiency point [6]. The 25 ◦ blade angle was chosen because it lies in the middle of the angle interval usually agreed between the producers and customers of bulb turbines. The blade angle 25 ◦ is often reinforced during performance or acceptance tests [18], because it corresponds to both common operating conditions and features increased probability for cavitation.

The bulb runner blade hydrofoil's shape determines the flow conditions in the runner. Well-shaped 3D blade hydrofoils reduce the inception and the growth of cavitation on the bulb turbine runner hydrofoils. Turbine hydrofoil profiles were used in designing the bulb runner blade's 3D shape, spanning from the runner hub toward the blade tip (Figure 2). For reliable 3D characterization of the turbine runner blade, and thus the entire bulb runner, several hydrofoils are usually selected. Later, a smooth surface is drawn over all hydrofoils to form the entire 3D geometry of the runner blade, further enabling CFD or experimental analysis.

**Figure 2.** Selection of a 25 ◦ angle on a modified bulb turbine hydrofoil.

' 's Although the bulb turbine blade can consist of many hydrofoils as explained above, in this paper we focus on a single hydrofoil. The most potential for energy conversion in bulb turbines is available near the tip, where tangential blade velocity and flow velocities are high. Due to high velocities, large pressure gradients, and the vicinity of the tip gap flow, this region is highly susceptible to cavitation occurrence and surface erosion. Here, several cavitation structures are usually found during acceptance tests. To focus on the cavitation properties of the hydrofoils, we have selected for further analysis the hydrofoil at 95% bulb turbine diameter D, measured from the runner rotation axis to the blade tip. The selected hydrofoil is thus located near the tip as shown in Figures 2 and 3. We made this selection on the basis of the energy transfer. Most of the energy is transferred where the turbine blade has a large diameter and the blade has a high surface area, and where the blade and flow velocities are high.

**Figure 3.** Bulb turbine blade design from individual 2D hydrofoils.

Hydrofoils have been studied extensively, including profile families such as NACA (among them 4–8 series), Göttingen, Munk, NPL, and others [19–21]. The available families of hydrofoils facilitate the turbine engineer's selection of the blade's hydrofoil to suit the required operating intervals and conditions, inter alia, for operation at high volumetric flows and runner blade angles. To enable operation in specific operating conditions, including cavitation properties, the turbine designer selects hydrofoil properties like leading

edge radius, camber, location of maximum thickness, maximum curvature, and trailing edge design.

The locations of cavitation occurrence on bulb turbine blades are well documented in [22]. In the hill diagram, we will focus on the high discharge region. Here, the cavitation is found in four regions [22]: the blade suction side, blade pressure side, fillet, and tip gap. The flow properties in and near the gap are dominated by turbulence and may appear in the entire hill diagram. The hydrofoil properties that most affect the cavitation on the blade suction and pressure side are the leading edge radius, the location of maximum thickness, and the curvature; we may also assume that the blade angle in bulb turbines is always optimal due to the turbine's dual regulation. We want to design a leading edge such that the flow around it is smooth so as to prevent any flow separation. As discussed above, we will focus on the hydrofoil at 95% of the blade radius, where we assume that in addition to the mentioned leading edge, maximum thickness and curvature influence, the cavitation properties are also influenced also by the presence of the gap flow. At 95% of the blade radius, the flow properties (velocity and pressure) are more pronounced because of higher blade velocity than near the hub. The sharp leading edges may introduce flow separation [23] and deteriorate the flow properties further downstream from the leading edge. Downstream from the leading edge the flow is determined by the camber shape, the position of maximum curvature, maximum thickness, etc. These parameters, among others, determine the velocity distribution in the blade channel from the pressure to the suction side. Velocity is related to absolute pressure; an increase in flow velocity decreases the local pressure [7], thus, the velocity distribution affects the cavitation properties. The pressure distribution may cause local cavitation inception, and both can influence boundary layer separation.

The hydrofoils in the bulb or any other water turbines are the base elements of the turbine runner blades and are thus rotating. The bulb turbine features pressure distribution with the pressure normally decreasing from the leading edge towards the trailing edge. According to the bulb turbine's triangles of velocity, downstream from the hydrofoil the relative velocity increases and reaches its maximum near the trailing edge [7], and the situation is remarkably different in comparison with the single hydrofoil in the water tunnel. Furthermore, the pressure distribution around the hydrofoil in the cavitation tunnel cannot fully represent the pressure properties around the bulb turbine blade. The situation in the bulb turbines is therefore very complex and the flow properties around individual hydrofoils cannot be directly correlated around the bulb turbine runner blades. Nevertheless, we will assume [6] that bulb runner cavitation can be related to the cavitation analysis of a single hydrofoil in the water tunnel.

#### *Hydrofoil Modifications*

The turbine blade hydrofoils are the primary component of the turbine blades. To improve performance, the turbine hydrofoil profiles should be smoothly and consistently remodelled in a spatial 3D redesign of the entire bulb turbine runner blade. The new blade differs from the old in the following parameters: the leading-edge radius has increased, the maximum thickness location has moved forward, and the location of the maximum hydrofoil curvature has moved forward. The characteristics of both turbine hydrofoils are shown in Table 1.


**Table 1.** Characteristics of hydrofoils as used on bulb turbine.

The modified hydrofoil as compared to the unmodified one features reduced absolute velocity in a section from the location of the maximum hydrofoil thickness and curvature to the trailing edge. The cavitation properties of the modified hydrofoil as compared to the unmodified one are in this section favorable. Because the points of maximum thickness and curvature for the modified hydrofoil profile were moved forward towards the leading edge, this section is relatively long.

The unmodified hydrofoil in a region from the leading edge to the location of the maximum thickness is very sensitive to the optimum angle of attack due to its sharp leading edge. Away from the optimum attack angle, the unmodified hydrofoil has disadvantageous cavitation characteristics in comparison with the modified one [6]. Contemporary water turbines often operate away from the best efficiency points due to the introduction of renewable energy sources and environmental changes. This trend will most probably increase in the future.

Two dimensional hydrofoils (Figure 4) were created from the 3D spatial bulb turbine blade hydrofoils (Figures 2 and 3). The 3D turbine blade hydrofoil was transformed from the cylindrical coordinate system of a turbine blade to a planar Cartesian x–y coordinate system as shown in Figure 3. –

**Figure 4.** Hydrofoils used in the cavitation tunnel analysis.

' As shown previously, both hydrofoils correspond to the runner blades' profiles at 95% R (Figure 3). The properties of both bulb runners are given in [4]. The turbine runners and both 2D hydrofoils were manufactured from brass using the same manufacturing techniques and both were given the same surface processing. The hydrofoils were later subjected to flow and image analysis in a cavitation tunnel and regression modelling.

The hydrofoil parameters for the experiment in the cavitation tunnel are shown in Table 2. The same manufacturing method as for the hydrofoils on the bulb turbine was used for the manufacture of the hydrofoils used in the cavitation tunnel.


**Table 2.** Hydrofoil parameters.

#### **3. Experiment**

Experiments were performed on the turbine measuring station for both bulb turbine runners and in the cavitation tunnel for both hydrofoils. A turbine measuring station was used to determine efficiency for a single runner blade angle. In the cavitation tunnel, cavitation properties were measured for analogous operating points.

#### *3.1. Turbine Measurements*

Measurements of efficiency were performed on a low-head axial and bulb turbines measuring station for low-head bulb turbines (Figure 5) in Kolektor Turboinštitut (Ljubljana, Slovenia). The measuring station was built according to standard [18] and was used for the model, witness, and acceptance tests. In agreement with standard [18], all measurements and data acquisitions were automatic. The nominal runner diameter of the test rig is 350 mm, and the recirculating flow is generated using two radial pumps with frequency regulation. The flow measurement was performed by using an electromagnetic flowmeter with a diameter of DN 400 and an accuracy of ±0.15%. The static pressure was measured using a differential pressure transducer (accuracy of ±0.025%), with the measuring taps located on the inlet and the outlet. Absolute suction pressure was measured using an absolute pressure transducer with an accuracy of ±0.025% to assure that the measuring station operated well above the incipient cavitation number. An electromagnetic incremental sensor was used to measure the rotational speed; the output was processed by a counter-timer module with an accuracy of ±0.01%. The shaft torque was measured on a horizontal shaft using a rotating torque transducer with an accuracy of ±0.01%. The friction torque was estimated by operation with air. 5) in Kolektor Turboinštitut

**Figure 5.** Turbine measuring station.

ψ Both bulb turbine runners were tested under similar conditions at six operating points. The blade angle was cases set to 25 ◦ in all cases (Figure 6). The blade angle was set using a control template. For comparison, both runners operated at the same specific energy E (J/kg) and the corresponding nondimensional energy numbers ψ were recorded. E was

determined from the total specific energy at the turbine inlet and the outlet cross-sections [7]. Energy numbers ψ were determined using the equation:

ψ

$$
\psi = \frac{\mathbf{2} \cdot \mathbf{E}}{\pi^2 \cdot \mathbf{N}^2 \cdot \mathbf{D}^2} \tag{1}
$$

Here N (min −1 ) is the rotational frequency, while D (m) is the discharge ring diameter of the bulb turbine. The turbine rotational frequency N (1000 min −1 ) was the same for all measuring points. Both runners had the same nominal diameter D (m). The guide vane opening coefficient A<sup>0</sup> was the same for each pair of measuring points used in the comparison: − −

$$\mathbf{A}\_0 = \frac{\mathbf{a}\_\mathbf{V} \cdot \mathbf{z}}{\mathbf{D}\_\mathbf{V}} \tag{2}$$

av (mm) is the guide vane's opening and z (/) is the number of guide vanes. Dv (mm) is the guide vane pivot diameter. Important manufacturing methods and characteristics were kept the same for reliable comparison between both runners. Control of runner tip gaps to the discharge ring and clearances of guide vanes to the guide vane pressure rings are very important for reliably estimating model turbine runners. Both runners were manufactured from the same material and using the same production methods and procedures. As such, for both runners, the tip gaps, guide vane clearances, and surface finishes were the same. D<sup>v</sup> e'

The measuring station was of a closed type. To maintain the ψ fixed, the main water circulation pump's frequency was controlled. The volumetric flow Q (m3/s) represented in its nondimensional form as flow number ϕ [18] is: The measuring station was of a closed type. To maintain the ψ fixed, the main water '

$$\boldsymbol{\upmu} = \frac{\mathbf{Q}}{\pi^2 / 4 \cdot \mathbf{N} \cdot \mathbf{D}^3} \tag{3}$$

The results from a previous study show that the increase in efficiency amounts to around 1% for the same runner geometry. Most of the increase according to the hill diagram is due to the small shift of the operating points towards low flow numbers and higher efficiency [6]. 

**Figure 6.** Turbine measuring points.

#### *3.2. Cavitation Tunnel Measurements*

The flow visualization measurements were taken on the cavitation measuring station. The approach to investigating the turbine blade flow properties in a separate tunnel experiment was similar to that used in [24]. A closed-type cavitation measuring station was used as shown in Figure 7. The flow was provided by a centrifugal pump driven

by a frequency variable drive. The volumetric flow measurement was performed using an electromagnetic transducer ABB COPA-XL. Absolute pressure pabs was measured on the suction side using an ABB 2600T 264NS pressure transducer. The temperature was measured by a Pt-100 4 wire temperature transducer, connected to the Agilent 34970A data acquisition unit. Valves in front of and behind the test section were fully open during experiments and they were only used during installation. Suction pressure was set using a vacuum pump, located on the top of a tank with a free water surface.

**Figure 7.** Cavitation measuring station with measuring and visualization equipment.

The cavitation tunnel test section cross-section was 50 cm × 100 cm. Hydrofoil width was 50 mm and chord length was l = 60 mm. The hydrofoil was installed at the angle of attack 8 ◦ .

A high-speed Photron SA-Z camera was used together with a Nikkor 105 mm 1:2.8 lens to visualize the cavitation in the cavitation section. The camera was mounted from the side. The camera was operated in a free-running mode. Images were recorded using PFV software in 8-bit BMP black and white mode at a resolution of 640 × 280. The framerate used was 25,000 frames/s. Shutter speed was equal to 10 µs. For each operating point, 2000 frames were recorded. For all experiments, all other camera settings were the same. A high intensity continuous LED illumination was provided by many CREE XLamp XM-L LED modules. LED modules were powered by an Instek DC power supply. Illumination was provided from the top of the test section above the hydrofoil on its suction side and arranged such that as few reflections as possible were present in the recorded sequences. Illumination intensity was not the same for all experiments; with the modified hydrofoil it was necessary to increase illumination intensity to visualize weaker cavitation structures. Keeping the illumination constant was impossible due to the low 8-bit depth of recorded images.

the Reynolds number Re and cavitation number σ. In Reynolds We have selected two dimensionless numbers to describe the cavitation in the cavitation tunnel, namely the Reynolds number Re and cavitation number σ. In Reynolds number Re

$$\text{Re} = \frac{\text{L} \cdot \text{v}}{\text{v}} \tag{4}$$

ν rofoil and ν is the L is hydrofoil length, v (m/s) is flow velocity in front of the hydrofoil and ν is the kinematic viscosity of water. To obtain the hydrofoil cavitation characteristics, the cavitation numbers σ were evaluated. Cavitation number σ was calculated according

ρ ∙ v

$$
\sigma = \frac{\mathbf{p}\_{\text{abs}} - \mathbf{p}\_{\text{v}}}{\rho \cdot \mathbf{v}^2} \tag{5}
$$

where p<sup>v</sup> (N/m<sup>2</sup> ) is the water evaporation pressure and pabs (N/m<sup>2</sup> ) is absolute pressure. Water evaporation pressure was estimate using the Equation

$$\mathbf{p}\_{\rm v} = 10^{2.7862 + 0.0312 \cdot \mathbf{T}\_{\rm W} - 0.000104 \cdot \mathbf{T}\_{\rm W}^2} \tag{6}$$

where T<sup>w</sup> ( ◦C) is water temperature. Equation (6) provides for a maximum difference of 0.2%, with a negligible influence on the measurements.

The cavitation σ-η curve is used to estimate the cavitation properties of bulb turbines. The cavitation point for the highest cavitation numbers is measured at ambient pressure, while for subsequent cavitation points suction pressure is decreased. For the tunnel measurements the cavitation number σ was used to observe the cavitation's development on both hydrofoils. For each experiment, the absolute pressure pabs was set as constant at the inlet of the cavitation tunnel. To decrease the cavitation number σ we increased the flow Q (m3/s) and flow velocity. Several such operating points were measured for increasingly lower absolute pressure pabs in the test section. This enabled comparison between both hydrofoils and among different absolute pressures.

The measuring points for both hydrofoils were selected according to cavitation structures sizes, shapes, and intensities. In general, cavitation intensity increased as cavitation number σ decreased. Instantaneous cavitation numbers σ were selected manually; we wanted to capture significant visible structural changes in cavitation. Such a procedure captured all important structures, starting with incipient cavitation on the leading edge, then sheet cavitation, cloud cavitation, and supercavitation across the entire hydrofoil [25,26]. Cavitation structure formation on both hydrofoils emerged at different cavitation numbers σ and, thus, the selection of cavitation number σ and corresponding measuring points for both hydrofoils was not the same.

#### *3.3. Cavitation Image Analysis*

Illuminated cavitation structures on the hydrofoils were recorded in a set of black and white images using the experimental setup from Section 3.2. The cavitation structures in the image have elevated grey level intensities, with vapor phase forming time and 2D space-dependent structures [22,27]. With computer-aided visualization of the process, the simultaneous time series of images were recorded. An assumption was made that the intensity of the illumination is proportional to the amount and intensity of the void fraction, being further proportional to the intensity of the cavitation structures in the observed regions [23]. This assumption is credible for cavitation intensity, low void fractions, and corresponding ratios of the intensity of the void fraction and grey level.

The time series of images from high-speed visualization in the cavitation tunnel were analyzed such that individual traveling bubbles were first filtered out. All bright objects in the individual image were deleted whenever they were very small. All other objects, like the attached cavitation or cavitation clouds, were retained. From these filtered images, the average grey level E<sup>g</sup> (i,j,t) and standard deviation were calculated as explained below.

The grey level variable E<sup>g</sup> (i,j,t) was acquired by processing intervals of 256 grey levels (8-bit camera image depth) from black (intensity = 0) to white (intensity = 255). The averaged grey level series were calculated in observation region. The standard deviation was estimated using the equation below:

$$\mathbf{S(i,j)} = \sqrt{\frac{1}{T} \sum\_{t=1}^{T} \left[ \mathbb{E}(\mathbf{i}, \mathbf{j}, \mathbf{t}) - \langle \mathbb{E}\_{\mathbf{g}}(\mathbf{i}, \mathbf{j}) \rangle \right]^2} \tag{7}$$

Here hE<sup>g</sup> (i,j)i is time and spatially averaged grey level, which is proportional to the average void fraction of the cavitation structure's.

$$
\langle \text{E}\_{\text{g}}(\mathbf{i}, \mathbf{j}) \rangle = \frac{1}{\text{T}} \sum\_{\mathbf{t} = 1}^{\text{T}} \text{E}(\mathbf{i}, \mathbf{j}, \mathbf{t}) \tag{8}
$$

The time interval t in Equations (7) and (8) is set from the duration and frequency of acquisition. 〈E<sup>g</sup> (i, j)〉 = T∑ E(i, j, t) t=1 

1

1

T

T

t=1

T∑ E(i, j, t)

(i, j)〉 =

〈E<sup>g</sup>

A sample of image analysis is shown in Figure 8 for the existing hydrofoil (p = 0.15 bar; Q = 61.7 m3/h; σ = 2.27) and in Figure 9 for the modified hydrofoil (p = 0.15 bar; Q = 61.7 m3/h; σ = 2.27). Cavitation length for the selected operating points was estimated to L = 22.74 mm for existing hydrofoil and L = 11.27 mm for modified hydrofoil. Such a method of estimating cavitation length was used for all operating points for both hydrofoils and later for regression cavitation length modelling. /h; σ /h; σ /h; σ /h; σ

/h, and σ **Figure 8.** Original hydrofoil, a sample image of cavitation (**top**), spatially averaged grey level (**middle**), and standard deviation of grey level (**bottom**), operating point p = 0.15 bar, Q = 61.7 m3/h, and σ = 2.27. /h, and σ

/h, and σ /h, and σ **Figure 9.** Modified hydrofoil, a sample image of cavitation (**top**), spatially averaged grey level (**middle**), and standard deviation of grey level (**bottom**), operating point p = 0.15 bar, Q = 61.7 m3/h, and σ = 2.27.

#### *3.4. Regression Cavitation Length Model*

Cavitation properties are a consequence of flow parameters, which are best described using Reynolds number Re and cavitation number σ. We will show that cavitation length on a hydrofoil in a cavitation tunnel may be well predicted based on these two dimensionless numbers. A good correlation is expected between the measured and calculated cavitation lengths. For this, we will use the following power-law regression model for cavitation length L σ

$$\mathbf{L} = \mathfrak{P}\_0 \mathbf{\cdot} \mathbf{R}^{\mathfrak{F}\_1} \mathbf{\cdot} \mathbf{s}^{\mathfrak{F}\_2} \ . \tag{9}$$

Model parameters β0, β1, and β<sup>2</sup> were calculated using the minimum mean square error method. Model parameters β , β and β

#### **4. Results and Discussion**

We have shown in Figure 6, that two bulb turbine runners with differently shaped blade hydrofoils (Figure 4) show very similar pressure characteristics. Because both runners are different, their hill diagrams and characteristics, in general, cannot be compared. The analysis of the cavitation properties of both hydrofoils will be shown in the following section. As an attempt to provide the background for the bulb runner evaluation with a focus on cavitation, two differently shaped hydrofoils were tested in the cavitation channel. We observed cavitation phenomena as a result of the hydrofoil shape design. The main tool for analyzing the cavitation structures appearing on the hydrofoil surface was computeraided visualization. From these cavitation images, the cavitation length was measured and used to make the regression model.

#### *4.1. Cavitation Dependence on the Cavitation Number σ*

L = β<sup>0</sup>

∙ Re

The visual comparison and analysis of the two hydrofoils as a direct result of decreasing cavitation number σ is shown in Figures 10–13. Cavitation numbers from high to low were selected such that a marked difference in cavitation structures occurrence was observed. The selected images represent to the greatest extent possible the most typical conditions for selected cavitation number σ. The analysis was limited to the hydrofoil suction side. A decrease in the cavitation number increases the intensity and amount of the cavitation structures. We will demonstrate that shape modification is essential for cavitation to occur on hydrofoils. A modified hydrofoil will demonstrate a decrease in cavitation phenomena intensity in comparison to an existing hydrofoil. The main area of visualization is over the hydrofoil's leading edge and the entire suction side of the hydrofoil. σ ing cavitation number σ is shown in – σ 's

**Figure 10.** Visualization results of cavitation structures at pabs = 0.65 bar.

**Figure 11.** Visualization results of cavitation structures at pabs = 0.45 bar.

**Figure 12.** Visualization results of cavitation structures at pabs = 0.30 bar.

– **Figure 13.** Visualization results of cavitation structures at pabs = 0.15 bar.

–

Figures 10–13 are limited to the intervals of the cavitation number that we were able to achieve in the cavitation tunnel with relevant cavitation intensities.

The measurements at absolute pressure pabs = 0.65 bar are shown in Figure 10. Cavitation structures (void fraction) start to appear on the modified hydrofoil at cavitation number around σ = 3.65. The cavitation starts behind the leading edge as an attached cavitation. With the further decrease of the cavitation number to around σ = 3.1 we observed an increase in the cavitation intensity and the quasiperiodic shedding of cavitation clouds. On the other hand, cavitation was observed on the original hydrofoil at cavitation numbers around σ = 4.45. At cavitation number σ = 3.65 we observed some cavitation cloud shedding, and at cavitation σ = 3.1 more than two-thirds of the hydrofoil was cavitating. At cavitation σ = 3.1, with the modified hydrofoil, the cavitation length was only 1/3 of the hydrofoil length. Since the hydrofoils in both cases were attached at the same angle of attack, we assume that the difference in cavitation occurrence is due to the shape and curvature of the leading edge. A sharp leading edge cuts through the flow such that it cannot follow the hydrofoil shape immediately downstream from the leading edge. The resulting flow separation induces cavitation inception due to the low absolute pressure and high velocity.

The measurements at absolute pressure pabs = 0.45 bar are shown in Figure 11. The cavitation on the modified hydrofoil profile starts behind the leading edge as an attached cavitation around σ = 3.45, and pabs = 0.65 bar. At cavitation σ = 3 we again observed cavitation cloud shedding as in Figure 10. The cavitation is much more intense for the case of the existing hydrofoil. We observed cavitation inception already at around σ = 4.5, while almost the entire hydrofoil was cavitating around σ = 3.

The measurements at the absolute pressure pabs= 0.30 bar are shown in Figure 12. Measurements commenced around σ = 3.55; at this operating point the incipient cavitation is already present. Incipient cavitation here is present as an attached cavitation. For the modified hydrofoil, the attached cavitation formed only at around σ = 3; it covered up to around half of the hydrofoil length. With the original hydrofoil, the cavitation inception began around the same cavitation number (σ = 3.55) as for the modified hydrofoil, although with slightly higher intensity. A remarkable difference in the intensity of cavitation was observed at cavitation number σ = 3.0. Here, the cavitation cloud is much larger for the existing hydrofoil. Below the cavitation number σ = 2.6 most of the hydrofoil is covered by a thick cavitation cloud in the case of the existing hydrofoil.

The measurements at absolute pressure pabs = 0.15 bar are shown in Figure 13. With modified hydrofoil, no cavitation was observed around cavitation number σ = 2.85, while with the modified hydrofoil the cavitation cloud at around σ = 2.25 covers up to one-third of the hydrofoil length. At around σ = 1.82 the cavitation cloud is again more intense than at cavitation number σ = 2.25. Below around cavitation number σ = 1.45 both hydrofoils are fully covered with supercavitation.

At cavitation numbers below around σ = 1.82 the large cavitation structures reduce the water passage area and increase the flow velocity in the cavitation tunnel at the selected flow. The results show favorable cavitation characteristics for the modified hydrofoil in comparison with the unmodified hydrofoil. We may also assume that the runner blade designed from such hydrofoils can perform well in the bulb turbine when operating under cavitation conditions.

#### *4.2. Cavitation Length Model*

Among several possible parameters of the cavitation flow properties, we have selected the cavitation cloud length for regression modelling. The results of the cavitation length model are shown in Table 3. Cavitation length dependence on the parameter β<sup>0</sup> is modest at around 6% in favor of existing hydrofoil. The modified hydrofoil's dependence on the parameter β<sup>1</sup> shows that the cavitation length of the modified hydrofoil is lower than that of the existing one. The most important improvement of the modified hydrofoil over the existing one is hidden in the dependence on cavitation number σ (achieving a value at parameter β<sup>2</sup> of −3.4, improving over the original −2.93). The results of the cavitation length model confirm the modified hydrofoil's superior cavitation characteristics over the existing one. parameter β − − 's superior cavitation characteristics parameter β − −

**Table 3.** Cavitation length model parameters for modelled cavitation length L.


The regression diagrams of the measured and modelled cavitation lengths are shown in Figures 14 and 15. The coefficients of determination for both the existing and modified hydrofoil profiles are reasonably good, with values above R <sup>2</sup> > 0.95. The high values of both coefficients of determination confirm the validity of the cavitation length model.

<sup>−</sup> −

<sup>−</sup> −

's superior cavitation characteristics

**Figure 14.** Cavitation model for existing hydrofoil, coefficient of determination R <sup>2</sup> = 0.9586.

**Figure 15.** Cavitation model for the modified hydrofoil, coefficient of determination R <sup>2</sup> = 0.9858.

#### **5. Conclusions**

In future power grids, waterpower flexibility will be an important parameter for ensuring network stability. Water turbines will need to operate at wide flow intervals and in transient regimes, and also allow for fast response times. For this reason, among others, further research on cavitation is required.

The goal of this research was to demonstrate the relation between the hydrofoil shape and cavitation phenomena. Existing and modified hydrofoils were extracted from runner blades on the 95% runner blade radius and modified. The modifications were made on the hydrofoil's leading edge, as well as at the points of maximum thickness and curvature. The leading-edge radius was enlarged to avoid flow separation. Maximum thickness and curvature were moved in the direction of the leading edge to stabilize the flow around the hydrofoil towards the trailing edge. Both hydrofoils were implemented in a bulb turbine and measured under equal energy numbers. Measurements of both turbine runners were performed on the measuring station according to the standard for acceptance tests of water turbines.

The results of this study show improved cavitation characteristics. However, we show here only a single modification of the hydrofoil and the results only demonstrate the efficiency of this modification. Of particular benefit for future design would be a comprehensive study, providing for a hydrofoil and runner shape optimization algorithm, based on selection and later minimization of the objective functions, including the manufacture and measurements of a selection of hydrofoils and runners.

To evaluate the cavitation properties, the flow was visualized in the cavitation tunnel on the suction side. The modified hydrofoil exhibited improved cavitation properties and we chose to evaluate the intensity of the cavitation using the cavitation cloud length. A regression model prediction of the cavitation length shows that modified hydrofoil features improved cavitation characteristics in comparison with the existing hydrofoil.

The primary goal of this study was to provide for improved cavitation characteristics on the modified hydrofoil. Using the modified hydrofoil invariably changes runner characteristics and this resulted in operating points shifting to the left in the hill diagram. This leftward shift resulted in increased efficiency as shown in [6].

We believe that future requirements for water turbines will emphasize flexibility. For this reason, the future research on cavitation in water turbines will focus on the minimization of cavitation occurrence and erosion, while at the same time providing adequate energy output and efficiency.

**Author Contributions:** Conceptualization, A.P., M.H. and M.D.; methodology, A.P., M.H., L.N. and M.D.; software, A.P. and L.N.; validation, A.P., M.H., L.N. and M.D.; formal analysis, A.P., M.H. and L.N.; investigation, A.P., M.H., L.N. and M.D.; resources, A.P., M.H. and M.D.; data curation, A.P., M.H. and L.N.; writing—original draft preparation, A.P. and M.H.; writing—review and editing, A.P., M.H., L.N. and M.D.; visualization, A.P. and M.H.; supervision, M.H. and M.D.; project administration, A.P. and M.H.; funding acquisition, M.H. and M.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was performed with financial support of Research Program P2-0401 Energy Engineering from Slovenian Research Agency ARRS.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Investigation of the Time Resolution Set Up Used to Compute the Full Load Vortex Rope in a Francis Turbine**

**Jean Decaix 1,\* , Andres Müller <sup>2</sup> , Arthur Favrel <sup>3</sup> , François Avellan <sup>2</sup> and Cécile Münch-Allig né <sup>4</sup>**

> 1 Institute of Sustainable Energy, School of Engineering, HES-SO Valais-Wallis, Rawil 47, 1950 Sion, Switzerland


**Abstract:** The flow in a Francis turbine at full load is characterised by the development of an axial vortex rope in the draft tube. The vortex rope often promotes cavitation if the turbine is operated at a sufficiently low Thoma number. Furthermore, the vortex rope can evolve from a stable to an unstable behaviour. For CFD, such a flow is a challenge since it requires solving an unsteady cavitating flow including rotor/stator interfaces. Usually, the numerical investigations focus on the cavitation model or the turbulence model. In the present works, attention is paid to the strategy used for the time integration. The vortex rope considered is an unstable cavitating one that develops downstream the runner. The vortex rope shows a periodic behaviour characterized by the development of the vortex rope followed by a strong collapse leading to the shedding of bubbles from the runner area. Three unsteady RANS simulations are performed using the ANSYS CFX 17.2 software. The turbulence and cavitation models are, respectively, the SST and Zwart models. Regarding the time integration, a second order backward scheme is used excepted for the transport equation for the liquid volume fraction, for which a first order backward scheme is used. The simulations differ by the time step and the number of internal loops per time step. One simulation is carried out with a time step equal to one degree of revolution per time step and five internal loops. A second simulation used the same time step but 15 internal loops. The third simulations used three internal loops and an adaptive time step computed based on a maximum CFL lower than 2. The results show an influence of the time integration strategy on the cavitation volume time history both in the runner and in the draft tube with a risk of divergence of the solution if a standard set up is used.

**Keywords:** vortex rope; Francis turbine; cavitation; CFD; RANS

#### **1. Introduction**

The rapidly growing share of new renewable energy sources in the past few decades leads to new challenges for the electrical network stability. Due to their rapid response time, hydraulic power plants are often used to compensate load fluctuations in the grid. However, the injection of the suitable amount of power requires the hydraulic machines to run outside of their Best Efficiency Point (BEP) for which they were designed. As a consequence, the flow leaving the runner possesses a significant residual swirl, leading to an inhomogeneous pressure distribution in the draft tube and eventually to the inception of cavitation.

For Francis turbines, running at an operating point for which the flow rate is above the nominal one leads to the development of an axial vortex rope. If the pressure level inside in the vortex is sufficiently low, cavitation occurs. For a sufficiently high Thoma number,

**Citation:** Lastname, F.; Lastname, F.; Lastname, F. Investigation of the Time Resolution Set Up Used to Compute the Full Load Vortex Rope in a Francis Turbine. *Appl. Sci.* **2021**, *11*, 1168. https://doi.org/10.3390/app11031168

Academic Editor: Florent Ravelet Received: 2 December 2020 Accepted: 20 January 2021 Published: 27 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the amount of cavitation is still small and the cavitating vortex rope is stable, whereas, by decreasing further the Thoma number, the vortex rope enters in a self-oscillating mode [1].

On one hand, various studies have been performed to better understand the selfoscillating behaviour of the cavitating vortex rope based on mathematical approaches [2] or 1D analyses [3], experimental investigations [4], 3D CFD [5] and coupling between 1D–3D numerical simulations [6,7]. On the other hand, several investigations look at identifying some key parameters (mainly the mass flow gain factor) useful for the calibration and the development of 1D hydro-acoustic models [8,9] using either experimental measurements [10,11] or CFD simulations [12–14].

Among the 3D CFD studies, the set up for the time integration is not often discussed since it is rare that the mean and the maximum Courant–Friedrichs–Lewy (CFL) numbers or the number of internal loops per time step are provided. However, due to the large range of values that can be taken by the speed of sound in a vapour/liquid mixture and the large time density gradient that can occur in a cell, the time integration requires attention. Usually, in order to improve the accuracy of a simulation at each time step, two ways are available: decreasing the time step or increasing the number of internal loops per time step (if such an option is available). The latter choice can provoke some instabilities [15,16]. In addition, by considering the Ansys CFX software, which is a coupled solver, the equation for the volume fraction is treated in a segregated manner and, in case of the simulation of a rotating turbomachinery, only a first order time scheme is available.

Based on these considerations, this paper investigates the influence of the time resolution set up used to compute an unstable cavitating vortex rope at full load using the Ansys CFX software. The sensitivity to the time integration set up is often disregarded due the large CPU effort required and authors focus rather on grid sensitivity. In addition, sensitivity analysis is often carried out at BEP. In the present paper, the authors focus on an unstable operating point of a Francis turbine at full load, for which a standard set up that allows for accurately capturing the stable cavitating vortex rope [17] is shown to not be able to capture the unstable one. The sensitivity of the time integration set up is performed either by increasing the number of internal loops per time step or by imposing a limit on the time step by limiting the maximum CFL number. The three simulation results are described and some points are discussed by comparisons with experimental data.

#### **2. Test Case**

The test case is a 1:16 reduced scale physical model of a Francis turbine with 16 runner blades and 20 guide vanes. The real prototype machine, featuring a rated power of over 444 MW at a specific speed of *ν* = 0.27 (dimensionless parameters are defined in Equation (1)), experiences severe pressure surges at high load conditions. The load corresponds to 130% of the value at the BEP, i.e., the discharge factor is equal to *Q*ED = 0.26 instead of 0.20 at the BEP. The speed factor is equal to the rated value of *n*ED = 0.288. The specific speed of the model runner at these conditions is *ν*<sup>m</sup> = 0.31. For a Thoma number *σ* = 0.11, the film displays an unstable cavitating vortex rope (see Figure 1) characterised by a self-excited breathing motion, i.e., a periodic collapse and reformation of the cavity accompanied by large amplitude pressure fluctuations throughout the system. Additional details regarding the experimental set up and the analysis of the measurements are available in [1,4].

The cited dimensionless parameters are defined as follows:

$$\nu = \frac{\omega}{2^{3/4} \pi^{1/2}} \frac{\mathbb{Q}^{1/2}}{\mathbb{E}^{3/4}}; \quad \mathbb{Q}\_{\rm ED} = \frac{\mathbb{Q}}{\mathbb{D}^2 \mathbb{E}^{1/2}}; \quad \mathbb{\eta}\_{\rm ED} = \frac{nD}{\mathbb{E}^{1/2}}; \quad T\_{\rm ED} = \frac{T\_m}{\rho \mathbb{D}^3 \mathbb{E}}; \quad \sigma = \frac{\text{NPSE}}{\text{E}} \tag{1}$$

with *ω* the rotational speed in rad/s, *Q* the flow rate in m3/s, E the total energy in J/kg equal to 262.83, *D* the outlet runner diameter in *m* equal to 0.35 m, *n* the runner frequency in s−<sup>1</sup> equal to 13.33, and *NPSE* the Net Positive Suction Energy.

Maximum extension of the vortex rope

Vortex rope collapse

**Figure 1.** Experimental visualization of the vortex rope at full load.

#### **3. Numerical Settings**

The flow is modelled based on the homogeneous Unsteady Reynolds-Averaged Navier–Stokes (U-RANS) equations assuming thermodynamic and mechanical equilibrium [18]. Discarding the energy equation, only the mass and momentum conservation equations are solved. The Reynolds stresses are modelled using the Boussinesq's assumption that introduces a turbulent viscosity *µ<sup>t</sup>* computed with the SST turbulence model [19]. At the wall, the turbulence model is coupled with an automatic wall function. The phase change due to cavitation is modelled using the model proposed by Zwart with the standard values for the parameters [20]. Additional details regarding the flow modelling are available in [17]. Other simulations have also been performed by Wack [21].

The computational domain considered includes a part of the inlet pipe, the spiral case, the stay vanes, the guide vanes, the runner and the draft tube (see Figure 2). The mesh is a structured hexahedral one based on a blocking technique and made with ICEM CFD (see Figure 2). The number of nodes is slightly above 11.5 million with 2.6 million of nodes in the runner domain and 4.3 million in the draft tube. The minimum angle is above 20 degrees and only 9% of the cells have an angle below 50 degrees. The *y* <sup>+</sup> value does not exceed 85 and the averaged value over each solid walls is lower than 20. A study of the grid sensitivity has been done and published in [17], showing around 2% of differences at the BEP between the computed *TED*, *QED* and *nED*.

The flow rate is imposed at the inlet, whereas the pressure is set at the outlet. The rotational speed of the runner is imposed on the experimental value. For all the solid walls, a no slip condition is imposed.

The convective fluxes in the momentum equation are computed with a high order scheme [22], whereas the turbulence convective fluxes are computed with an upwind scheme. Regarding the time scheme, the second order backward Euler scheme is applied [23] excepted for the volume fraction equation that is discretized with a first order backward Euler scheme (In Ansys CFX 17.2, for rotating turbomachineries, this option is the only one available.). The time step, the CFL number, the number of internal loops per time step and the time duration of the simulation expressed as the number of runner revolutions achieved are given in Table 1 for each case considered. *Case 1* corresponds to the standard set up in Ansys CFX for a flow in a Francis turbine with a time step set

to less than one degree of revolution per time step and five internal loops per time step. This set up has been validated against experimental data for stable cavitating vortex rope in [17]. *Case 2* differs from *Case 1* only by the number of internal loops, which is set to 15. *Case 3* differs from the other two cases since the time step is not imposed by the user but is computed automatically by imposing a limitation on the CFL number to a value of 2. This value is sufficiently small to deal with the accuracy of the first order backward Euler scheme used for the volume fraction equation keeping advantage of an implicit formulation. In addition, the number of internal loops is decreased to a value of 3 even if two internal loops are enough to reach a second order time accuracy (see [24]) for the equations discretized with a second order scheme. The unsteady simulations are started from a converged steady simulation (see [17]).

**Figure 2.** Computational domain (left top), view of the surface mesh of blade (left bottom) and side view of the structured mesh in a mid-plane (right).

Even if CFX used a coupled solver, the volume fraction equation is treated in a segregated manner. The interfaces between the stationary and rotating domains are handled using a General Grid Interface (GGI) algorithm.

**Table 1.** Parameters set for the time integration. Bold terms refer to the imposed values by the user.


#### **4. Results**

✾✺ ✾✻ ✾✼ ✾✽ ✾✾ ✶✵✵ ✶✵✶ ✶✵✷ ✶✵✸

✶✵✹ ✶✵✺ ✶✵✻ ✶✵✼ ✶✵✽ ✶✵✾

First of all, it has to be mentioned that the simulation *Case 1* diverges when the vortex rope in the draft tube collapses (see the next section for a detailed discussion). However, until this point, all the simulations show the development of a cavitating vortex rope and a cavitation sheet over the trailing edge of the suction side (see Figure 3). Such features are supported by experimental visualisation [4]. Compared to the experiment (Figure 1), the volume of vapour inside the vortex rope is underestimated.

**Figure 3.** Instantaneous iso-surface of the liquid volume fraction *α<sup>L</sup>* = 0.9. Downstream view (top) and side view (bottom).

#### *4.1. Integrated Variables*

The time-history of the efficiency, the dimensionless volume of vapour in the computational domain, the torque factor *TED* and the speed factor *nED* are plotted on Figure 4 (the experimental values for the *TED* and *nED* are added in the caption of the figure). The divergence of the simulation *Case 1* just after the collapse of the cavitation volume is clearly put in evidence. On the contrary, simulations *Case 2* and *Case 3* do not show a complete collapse of the cavitation volume even if a time oscillation of the different quantities is observed. The first cycle seems to be shortened by increasing the requirement on the time integration since it is shorter for *Case 2* compared to *Case 1* and the shortest for *Case 3*.

Regarding the efficiency, the torque factor and the speed factor, the differences between the simulations are around 1.5%, whereas, for the volume of vapour, simulation *Case 1* predicts a volume that is at least four times higher. Compared to the measurements, the *TED* and *nED* are overestimated respectively by approximately 2.5% and 3.5%. Simulations *Case 2* and *Case 3* have a rather similar behaviour even if a time shift is observed due to a faster first cycle for *Case 3*. Furthermore, due to the smaller time step imposed in *Case 3*, fluctuations at higher frequencies are observed in all the signals.

For simulation *Case 2* between the seventh and eleventh runner revolution, the slow variation of the efficiency, the torque factor, and the speed factors disappear, whereas a fluctuation at a low frequency is still observed for the volume of vapour. By splitting the time history of the volume of vapour between the volume in the runner and draft tube domains (see Figure 5), it is noticeable that between the seventh and eleventh runner revolution, the volume of vapour in the draft tube is close to zero and then increases again after the tenth runner revolution. Therefore, the low frequency seems to be related to the volume of vapour in the draft tube, which means with the presence of a cavitating vortex rope. This point seems to be in agreement with the experimental investigations [1].

**Figure 4.** Time-history of the efficiency, the dimensionless volume of vapour, the torque factor *TED* (measured value: 0.127) and the speed factor *nED* (measured value: 0.288).

**Figure 5.** Time-history of the dimensionless volume of vapour in the runner (left) and in the draft tube (right).

A Fourier analysis of the efficiency signal is shown in Figure 6 on the left. The frequency is divided by the runner frequency *fn*. The modes put in evidence are the ones associated with the blade passage frequency *f<sup>b</sup>* = 16 *f<sup>n</sup>* and the first and fifth harmonic. Both *Case 2* and *Case 3* capture the fifth harmonic, which corresponds to the rotor/stator interaction since 5 *f<sup>b</sup>* = 4 *fgv* = 80 *fn*. A Fourier analysis of the volume of vapour has also been carried out even if only the simulation *Case 2* has achieved enough cycles to provide an estimation of the low frequency characteristic of the time evolution of the volume of vapour. By focusing only on *Case 2*, a low frequency around 0.3 *f<sup>n</sup>* is put in evidence, which is closed to 0.25 *fn*, the value observed on the torque during the measurements [1].

**Figure 6.** Magnitude of the Fourier transform of the efficiency (left) and volume of vapour (right) signals.

#### *4.2. Local Pressure*

Several pressure probes have been set in the computational domain as shown in Figure 7. Two probes are located in two different guide vane channels, six are located on a runner blade with three on the pressure side and three others on the suction side, two times four probes are distributed at ninety degrees on the wall in two sections of the draft tube cone (as in the experiment) and finally two probes are set in the draft tube elbow.

**Figure 7.** Positions of pressure probe in the computational domain (left). Zoom on the runner probe locations (right).

For each probe, the fluctuations of the pressure coefficient are computed by subtracting the average pressure at the probe location and then by dividing the result by 0.5 *ρ U*<sup>2</sup> 1¯ (with *U*1¯ the peripheral velocity at the runner outlet). In addition, for the probes in the cone of the draft tube, a spatial averaging is also performed.

The pressure signals on the pressure side of the runner at the leading edge (probe Ru1) and on the suction side of the trailing edge (probe Ru6) and in the two sections of the draft tube cone are displayed in Figure 8. In the draft tube cone, the experimental pressure signals are also added.

At the leading edge of the runner, the pressure fluctuations captured by the three simulations are identical with a high frequency content. At the trailing edge, the fluctuations at high frequency alternate with a flat period corresponding to the period during which the probe is covered by a pure vapour sheet. No clear pattern can be observed. In the draft tube cone, compared to the experiment, the simulations *Case 2* and *Case 3* are smoother without a strong peak of pressure after the partial collapse of the vortex rope. Simulation *Case 1* presents a strong pressure peak at the collapse of the vortex rope, which seems to follow qualitatively the trends of the experimental one. However, as already mentioned, the simulation diverges after the collapse.

**Figure 8.** Time history of the pressure fluctuations in the runner and the draft tube cone. The zero of experimental signal has been set arbitrary.

The Fourier transforms of the pressure signals in the runner and the draft tube are represented in Figure 9. At the leading edge of the runner, the guide vane passage frequency *fgv* = 20 *f<sup>n</sup>* and its harmonic until the fourth are clearly observed. At the trailing edge, only the fundamental is observed. In addition, a frequency around 120 *f<sup>n</sup>* is captured by the simulations *Case 2* and *Case 3*. This frequency is also observed in the cone of the draft tube but with an attenuation. The origin of this frequency is not clear.

In the draft tube, the comparison of the simulation *Case 2* with the experiment makes the presence of the low frequency appear around 0.3 *f<sup>n</sup>* in the simulation instead of 0.25 *f<sup>n</sup>* in the experiment. Obviously, the amplitude is low compared to the experiment due to the

inability of the simulation to capture accurately the collapse of the vortex rope as shown in Figure 8.

**Figure 9.** Magnitude of the Fourier transform of the pressure signal in the runner and the draft tube.

#### **5. Discussion**

The results described in the previous section make a question appear—why the simulation *Case 1* seems to follow more closely the behaviour of the vortex rope observed experimentally than the simulations *Case 2* and *Case 3*, which use a refined time integration set up.

The Root Mean Square (RMS) and the maximum RMS of the pressure equation are compared in Figure 10. The RMS reach values lower than 10 −3 for all the simulations except simulation *Case 1* before the simulation crashes. For *Case 2*, the residuals reach values below 10 −5 , which shows the role plays by the number of internal loops per time step. The maximum RMS are quite high for all simulations with some values above 0.1. Interestingly, the maximum RMS of *Case 1* and *Case 3* are around 0.4, but *Case 3* does not crash.

**Figure 10.** RMS (left) and maximum RMS (right) of the pressure equation.

The time history of the dimensionless mass flow *Q*<sup>∗</sup> = *Q*/*Qinlet* at the inlet of the runner and at the outlet of the computational domain are compared for the three simulations in Figure 11. In the runner inlet section, upstream from the region where cavitation occurred, the mass flow is constant and equal to 1. Downstream, due to the phase change, the mass flow oscillates. The oscillation is stronger for *Case 1* with a strong decrease of the mass flow, which is a consequence of the vapour volume growth (see Figure 5).

**Figure 11.** Time history of the dimensionless mass flow at the runner inlet section (left) and at the outlet of the computational domain (right).

By assuming a density of the vapour negligible compared to the liquid density, the 1D mass conservation equation applied for the whole computational domain writes:

$$Q\_1 - Q\_2 = -\frac{dV}{dt} \tag{2}$$

with *Q*<sup>1</sup> and *Q*2, respectively, the mass flow rate at the inlet and outlet of the computational domain and *V* the volume of vapour. The left- and right-hand terms of Equation (2) are plotted in Figure 12 for the three cases. For *Case 2* and *3*, the mass conservation is respected contrary to *Case 1* for which the time variation of the volume of vapour is smaller than the difference of the flow rate. This mass unbalanced seems to be responsible for the collapse of the vortex rope and the crash of the simulation. As the boundary condition are fixed, the mass flow rate unbalanced cannot be supported infinitely since it is a non-equilibrium state. Therefore, the collapse of the vortex is the only way to return to an equilibrium state.

**Figure 12.** Time history of the left- and right- hand terms of Equation (2) for each case.

#### **6. Conclusions**

The influence of the time resolution set up on the simulation of an unstable cavitating rope has been performed using the Ansys CFX software. A standard set up *Case 1* usually used for the URANS computation of a Francis turbine is not able to capture the flow due to the crash of the simulation after the first collapse of the cavitating vortex rope. By increasing the number of internal loops from 5 to 15, the simulation *Case 2* is able to capture several cycles of the vortex rope with a dampened amplitude compared to the experiment but with a frequency close to the measured one. It has also been noticed that the low frequency observed on the torque or the efficiency seems to be related to the presence of a cavitating vortex rope in the draft tube since, in its absence, the low frequency disappears. The simulation *Case 3* carried out with an imposed maximum CFL number lower than 2, i.e., with a smaller time step, and only three internal loops does not provide a dynamic behaviour very different than *Case 2* even if a longer simulation could show differences.

Surprisingly, the simulation *Case 1* is the only one that gives a duration and an amplitude of the pressure peak in agreement with the experiment. Nevertheless, the simulation stops when the vortex rope collapse due to a solver crashes and no successive cycles are computed. On the contrary, the simulations *Case 2* and *Case 3* do not show any growth or collapse of the cavitating vortex rope. This difference is explained by the fact that the simulation *Case 1* does not respect the mass conservation equation at each time step contrary to the two other simulations.

To compute an unstable vortex rope, it is therefore necessary to take care of the time integration set up by increasing at least the number of internal loops. To improve the accuracy of the simulation, the U-RANS simulation needs, in addition, to be coupled with a 1D hydroacoustic model allowing for adjusting the boundary conditions at the inlet and outlet boundaries of the computational domain [6,7].

**Author Contributions:** J.D. contributed to carrying out the simulations and the post-treatment of the CFD results. He is the main writer of the paper. A.M. and A.F. performed the experimental work and provide the experimental results used in this paper. They also contribute to the redaction of the paper. F.A. and C.M.-A. were responsible for the supervision, project administration, and the funding acquisition. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research leading to the results published in this paper is part of the HYPERBOLE research project, granted by the European Commission (ERC/FP7-ENERGY-2013-1-Grant 608532)

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study could be available on request from the corresponding author. The data are not publicly available due to NDA.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


## *Article* **Numerical Investigation of Methodologies for Cavitation Suppression Inside Globe Valves**

**Jun-ye Li <sup>1</sup> , Zhi-xin Gao 1,2,\*, Hui Wu <sup>1</sup> and Zhi-jiang Jin <sup>2</sup>**


Received: 20 July 2020; Accepted: 7 August 2020; Published: 11 August 2020

**Abstract:** Cavitation inside globe valves, which is a common phenomenon if there is a high-pressure drop, is numerically investigated in this study. Firstly, the cavitation phenomenon in globe valves with a different number of cages is compared. When there is no valve cage, cavitation mainly appears at the valve seat, the bottom of the valve core, and the downstream pipelines. By installing a valve cage, cavitation bubbles can be restricted around the valve cage protecting the valve body from being damaged. Secondly, the effects of the outlet pressure, the working temperature, and the installation angle of two valve cages in a two-cage globe valve are studied to find out the best method to suppress cavitation, and cavitation number is utilized to evaluate cavitation intensity. Results show that cavitation intensity inside globe valves can be reduced by increasing the valve outlet pressure, decreasing the working temperature, or increasing the installation angle. Results suggest that increasing the outlet pressure is the most efficient way to suppress cavitation intensity in a globe valve, and the working temperature has a minimal effect on cavitation intensity.

**Keywords:** cavitation; cavitation number; globe valve; valve cage; computational fluid dynamics

#### **1. Introduction**

For the valves applied in the industrial pipelines, there is often a reduction of the flow channel around the valve seat or the valve plug, where the liquid pressure may fall below the saturated vapor pressure and the cavitation phenomenon may occur. Hence cavitation is one of the reasons that can lead to valve failure if a valve is not designed appropriately [1].

Valves play important roles in many working scenarios, and common applications include flow, pressure, and temperature control; common valves including ball valves [2], control valves [3], and globe valves [4], etc. are summarized and shown in Table 1. Works about valves mainly focus on the flow characteristics and structure optimization. Qian et al. [5] investigated the acoustic power with different thermal conditions in a temperature and pressure regulation valve. Fratini and Pasta [6] investigated the effects of the residual stress on the fatigue crack. Wang et al. [7] designed a new spool structure and various pressure and valve openings were investigated to find out their effects on the stability of the flow field. Lin et al. [8] studied the two-phase flow characteristics of a gate valve and the effects of size and numbers of the particles on erosion were analyzed. Jin et al. [9] and Qian et al. [10–12] applied the Tesla valves in the hydrogen decompression process and the flow characteristics and the aerodynamic noise of a single and a multi-stage Tesla valve was investigated. Fan et al. [13] focused on the contaminated friction of a spool valve. Qian et al. [3] investigated the flow characteristics of a feed-water valve with different structures of the throttling window. Yan et al. [14] focused on the erosion in a servo valve. Pasta et al. [15] and Rinaudo et al. [16] investigated the hemodynamics in

biomedical heart valves and the results could be helpful for the clinical diagnosis. Zhang et al. [17] studied the effects of the damping sleeve on the flow force of an on-off valve. Zhong et al. [18] proposed a two-level fuzzy control algorithm based on the calculated flow rate feedback from the spool displacement to realize high precision flow control.

**Table 1.** Common valves applied in different industries.

<sup>1</sup> Reproduced from reference [2]. Copyright 2017, Elsevier.

Cavitation that appears in the valve can lead to self-excited noise, cavitation erosion, and system jam, etc. [19,20], and it may even cause severe damage to the valve spool as shown in Figure 1 [21]. To date, more and more works focusing on the cavitation problem in various kinds of valves can be found, and an experimentally verified computational fluid dynamics (CFD) method is becoming a popular tool [22,23]. Wen et al. [24] studied the velocity of a needle valve motion on cavitation using the dynamic grid overlapping techniques and a visual experimental platform. Li et al. [25] focused on the cavitation in bileaflet mechanical heart valves and their results showed that it was induced by the high closing velocity and water hammer effect. Jin et al. [26,27] investigated the cavitation in a globe valve and a sleeve regulating valve, and the effects of the valve core shape and the inlet velocity were obtained. Liu et al. [28] studied the effects of inlet pressure on cavitation location and area in a regulating valve, and Qiu et al. [29] focused the effects of valve opening and pressure drop of a sleeve regulating valve. Liu et al. [30] investigated the cavitation erosion in a butterfly valve with various inlet pressure and valve opening angles. Jin et al. [31] numerically analyzed the cavitating flow in microchannels and the effects of the diameter ratio and pressure drop were discussed. Zhang [32] studied the flow characteristics of a cone throttle valve under cavitation, and the cavitation luminescence was used to help relative experiments. Yuan et al. [33] compared the compressible and incompressible method, and their results showed that the compressible method could provide more precise results.

**Figure 1.** Damage of the valve plug results from cavitation erosion [21]-Copyright 2016, Elsevier.

To suppress or prevent the occurrence of cavitation, valve structure optimization, adding outlet pressure, or adding stage/anti-cavitation trims is inevitable. Yaghoubi et al. [4] numerically studied the effects of the numbers of valve trims on cavitation intensity, and their results showed that the numbers of valve trims should be below two. Xu et al. [34] performed structure optimization of the valve port in a relief valve to reduce cavitation intensity. Shi et al. [35] utilized a drainage device to suppress the cavitation in a throttle valve, and cavitation intensity could be reduced with appropriate inlet and outlet pressure. Lee et al. [36] optimized the structure of the bottom plug in a three-way reversing valve to suppress valve cavitation. Chern et al. [37] found that by adding a valve cage, cavitation could be limited to the region close to the valve cage, and the valve body could be protected from damage resulting from valve cavitation.

Globe valves are widely used in nuclear and thermal power stations, or in petrochemical fields, where there is often a high-pressure drop condition resulting in cavitation inside globe valves [37]. In this paper, a globe valve with two valve cages is proposed to suppress valve cavitation, and each cage has 20 evenly distributed slender orifices. Firstly, the flow characteristics of the globe valve with and without valve cage are compared to investigate cavitation distribution and cavitation intensity. Then the effects of the working temperature and the valve outlet pressure on cavitation intensity for the globe valve with two valve cages are studied. At last, the effects of the installation angle of two valve cages are analyzed and the most efficient method to reduce cavitation intensity is found out.

#### **2. Numerical Methods**

The structure of the investigated globe valve is shown in Figure 2; here the globe valve with a different number of cages and the geometry of the cage are displayed. The two-stage globe valve is self-designed to be used in the pipe system with a high-pressure drop, and the model is created using commercial software SOLIDWORKS. The nominal pipe diameter of the globe valve, *D*, is 50 mm, and the entrance pipe length is 3 *D* and the downstream pipe length is 8 *D*. The throttling orifices on each cage are the same and are evenly distributed along the circumference. The throttling orifices on two cages can either be installed aligned or be installed with an installation angle, α, as shown in Figure 2. *α*

**Figure 2.** Structure of the investigated globe valve.

*ε* Fluid flow through the two-cage globe valve with a fully opened state is numerically solved. Liquid water and water vapor under different temperatures are considered as the incompressible working fluid, and the properties of the working fluid are shown in Table 2. The continuity and the momentum equations are solved using commercial software ANSYS Fluent, and the governing equations and the equations of the Realizable *k-*ε turbulence model can be found in the previous work [26]. The cavitation model proposed by Zwart et al. [38] is used and the equations describing cavitation are shown below.

$$\frac{\partial}{\partial t}(\alpha \rho\_{\upsilon}) + \nabla(\alpha \rho\_{\upsilon} u\_{\upsilon}) = R\_{\epsilon} - R\_{\epsilon} \tag{1}$$

*Appl. Sci.* **2020**, *10*, 5541

$$\stackrel{\circ}{\partial}\_{\mathcal{O}}(a\rho \;) + \nabla \begin{pmatrix} a\rho & \\ \end{pmatrix} = \ -$$

$$\mathbb{R}\_{\mathbf{r}} = F\_{\text{vap}} \frac{\partial \alpha\_{n[\mathbf{r}]}(\mathbf{1}\cdot\mathbf{c}\alpha\lambda)\rho\_{\mathbf{v}}}{R\_{B}} \sqrt{\frac{2}{3}\frac{p\_{\mathbf{v}}\cdot\mathbf{p}}{\beta\lambda}} \quad (\mathfrak{p} \le \mathfrak{p}\_{\mathbf{v}}) \quad \mathbf{j} \tag{2}$$

$$R\_{\mathfrak{k}} = F\_{cond} \frac{\mathfrak{Z} \rho\_p \rho\_v}{R\_B} \sqrt{\frac{2}{3} \frac{p = p\_v}{\beta^l}} \quad (p \ge \underline{p}\_v) \quad \tag{3}$$

here α is vapor volume fraction, subscript *v* stands for vapor phase, ρ stands for density, *u* stands for velocity, *R<sup>e</sup>* and *R<sup>c</sup>* are related to the growth and collapse of the vapor bubbles, *RB*, *Fvap*, and *Fcond* are model constant. During the simulation, the pressure-based SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm is utilized, and the second-order upwind discretization scheme is applied. The residuals are set as 0.001 and the area-weighted inlet velocity is monitored, and the solution is considered as converged if the residual is below 0.001 and the inlet velocity is constant. *α ρ*


**Table 2.** Properties of water and water vapor.

Pressure inlet and pressure outlet are applied as the main boundary conditions. The inlet pressure is 1.17 MPa, and the outlet pressure varies from 0.1 to 0.25 MPa. The no-slip wall with the scalable wall function is used. The temperature variation is neglected during the simulation. −3 − −5 −5 −5

A polyhedral mesh with boundary layer is utilized in this study, and the zone near the throttling window of the cage is refined. Figure 3 shows the mesh of the investigated two-cage globe valve. Mesh independence is done with a different number of grid cells when the outlet pressure is 0.25 MPa and the temperature is 80 ◦C. The inlet velocity and pressure drop are compared and shown in Table 3, and it is clear that the results are close when the cell number is 824,206 and 1,550,930, thus 824,206 cells are adopted eventually.

**Figure 3.** Mesh of the two-cage globe valve.



#### **3. Results and Discussion**

The globe valve with and without cages is compared in this study, and the installation angle between two cages, the outlet pressure, and the working temperature are investigated. Flow on the middle plane of the valve body including the velocity and pressure distributions and the streamlines are analyzed to help understanding the cavitation phenomenon. Cavitation number, loss coefficient, and induced water vapor volume are quantitatively studied and compared.

#### *3.1. Comparison of Cavitation in Globe Valves with and without Cages*

When there is a large pressure drop, local pressure is very likely to drop under the saturated vapor pressure. Here, the globe valve without a cage, the globe valve with one cage and the globe valve with two cages are simply compared to find out the effects of the cage. The outlet pressure is chosen as 0.1 MPa and the working temperature is 80 ◦C.

The velocity and pressure distributions in the globe valve with a different number of cages are shown in Figure 4, and an obvious difference in the velocity and pressure distributions can be found between different globe valves. The Reynolds number is around 7 <sup>×</sup> <sup>10</sup><sup>5</sup> , which means the flow regime is turbulent. At the bottom of the valve body, there is a high-pressure and low-velocity zone, which results from the sudden change of the flow direction. When there is no cage, most water flows directly toward the valve downstream, resulting in a small high-pressure zone at the bottom of the valve core closing to the valve outlet. While cages are added, the valve resistance increases, and the pressure distribution at the valve core bottom becomes more even. From Figure 4, it can also be found that the cage number has little influence on the streamline distribution once the cage is added.

**Figure 4.** Velocity and pressure distributions for the globe valve with a different number of cages.

When water vapor volume fraction is above 0.1, the water vapor volume distribution in different globe valves is compared qualitatively [4,37] and shown in Figure 5. It is clear that cavitation caused water vapor mainly congregates around the valve seat and the bottom of the valve core if there is no cage, but the water vapor only appears at the location near the valve outlet, which is because of the high velocity and low pressure as shown in Figure 4. Besides, a relatively large water vapor volume can be found at the location where the valve body and the downstream pipe connect. When a cage is added, the region where water vapor appears varies and is limited close to the cage itself, which protects the valve body from being damaged. When two cages are added, the region where cavitation appears reduces further but is still near the two cages.

**Figure 5.** Water vapor volume distribution for globe valve with a different number of cages: (**a**) no cage; (**b**) 1 cage; (**c**) 2 cages.

*σ* To compare cavitation intensity in a globe valve with a different number of cages, dimensionless cavitation number σ is used. The expression of the cavitation number is shown below.

$$
\sigma \equiv \overline{\sigma} \frac{p\_o \not p \not p \tau}{1/2 \rho v^2} p
\tag{4}
$$

*ρ* here *p<sup>o</sup>* is the pressure at valve outlet, *p<sup>v</sup>* is the saturated vapor pressure, ρ is the density of water, *v* is the velocity at the valve outlet.

*ξ* The loss coefficient ξ is also calculated to analyze the flow characteristics of different globe valves. The expression of the loss coefficient is shown below.

$$\xi \underset{\overline{\xi}}{\rightleftharpoons} \frac{p\_i}{\overline{\xi}} \overline{\text{p}} \frac{p\_a}{2\rho v^2} p$$

here *p<sup>i</sup>* is the pressure at valve inlet.

The comparison of the cavitation number and the loss coefficient between different globe valves is shown in Figure 6. It can be found that the cavitation number and the loss coefficient increase with the increase in the valve cage number, which means that cavitation intensity can be suppressed by adding a valve cage in a globe valve. In Figure 5, it is noticed that the water vapor volume in the globe valve with two cages is higher than other globe valves. The reason is that when the total pressure in the valve inlet is the same, the velocity inside the globe valve with two cages is lower due to high resistance, which leads to a higher static pressure at valve inlet. Thus, the actual water vapor volume is higher in the globe valve with two cages. Together with Figures 5 and 6, it can also be found that the cavitation number is useful when deciding the cavitation intensity in valves.

**Figure 6.** Comparison of cavitation number and loss coefficient for the globe valve with a different number of cages.

#### *3.2. E*ff*ects of the Working Temperature and Outlet Pressure*

When the outlet pressure and the working temperature vary, the total pressure drop or the saturated vapor pressure changes, thus leading to the variation of cavitation intensity. In this section, the effects of the working temperature and the outlet pressure are studied and compared.

Figures 7 and 8 show the water vapor volume distribution in the two-cage globe valve under different outlet pressures and different working temperatures. The working temperature for Figure 7 is 80 ◦C, and the outlet pressure for Figure 8 is 0.1 MPa. The water vapor volume decreases with the increase in the outlet pressure or the decrease in the working temperature. Based on the water vapor volume distributions, one can find that increased outlet pressure has more positive effects on cavitation intensity than the working temperature.

**Figure 7.** Water vapor distribution in two-cage globe valve with different outlet pressure: (**a**) 0.15 MPa; (**b**) 0.20 MPa; (**c**) 0.25 MPa.

**Figure 8.** Water vapor volume distribution in two-cage globe valve with different working temperature: (**a**) 20 ◦C; (**b**) 40 ◦C; (**c**) 60 ◦C.

To investigate the effects of the outlet pressure and the working temperature quantitatively, the cavitation number and the loss coefficient are compared in Figure 9, where the straight line stands for cavitation number, and the dash line stands for the loss coefficient. It can be found that the cavitation number increases with the increase in the outlet pressure and the decrease in the working temperature, which is consistent with the above results. Cavitation number almost has a linear relationship with the outlet pressure. When the outlet pressure increases by 0.05 MPa, the cavitation number increases nearly by 0.4. While the increment of the cavitation number decreases with the decrease in the working temperature. When the working temperature decreases from 80 ◦C to 60 ◦C, the average increment of cavitation number is 0.134, and the average increment of cavitation number is 0.0538 if the working temperature decreases from 40 ◦C to 20 ◦C.

**Figure 9.** Cavitation number and loss coefficient of two-cage globe valve with different outlet pressure and working temperature.

The loss coefficient barely changes with varying outlet pressure which is because the loss coefficient is affected by the pressure difference between the valve inlet and valve outlet shown in Equation (5). However, with the increase in the working temperature, the loss coefficient decreases firstly then increases rapidly, and when the working temperature is 80 ◦C, the loss coefficient has the maximum value. The difference of the loss coefficient may be caused by the cavitation phenomenon inside valves.

In the numerical simulation applied in this study, the finite volume method is used. Thus, the water vapor volume in valves can be calculated based on the water vapor fraction and the discrete volume, and the expression is shown below. <sup>=</sup> ∑

α

$$V = \sum\_{1}^{n} \alpha\_{\upsilon} dV \tag{6}$$

here *V* stands for the total water vapor volume in a valve, *n* stands for the total number of the discrete volume, *dV* stands for the volume of a discrete volume, α*<sup>v</sup>* stands for the water vapor fraction in a discrete volume. The water vapor volume in the investigated two-cage globe valve under different outlet pressures and the different working temperatures is shown in Figure 10. It can be found that when the outlet pressure is larger than 0.2 MPa, the variation of water vapor volume under different working temperatures is small. While the outlet pressure is smaller than 0.2 MPa, the water vapor volume increases with the increase in the working temperature and the outlet pressure, and the higher the working temperature, the greater the influence of the outlet pressure on water vapor volume. Together with Figure 8, it can be concluded that increase in the outlet pressure is more effective in suppressing the cavitation than the decrease in the working temperature. *α*

**Figure 10.** Water vapor volume in the two-cage globe valve with different outlet pressure and working temperature.

#### *3.3. E*ff*ects of the Installation Angle*

When there are two cages inside a globe valve, the installation angle between the two cages influences the flow resistance, thus affecting the cavitation intensity. For the investigated two-cage globe valve, the number of the throttling window on each cage is 20, thus the installation angle from 0 ◦ to 9◦ is investigated. Figure 11 shows the streamline on a cross-section in the two-cage globe valve under different installation angles. It can be seen that when the installation is 0 ◦ , most water flows directly from the inner cage to the outer cage, resulting in small flow resistance. With the increase in the installation angle, water comes from the inner cage and collides with the outer cage, and then changes flow direction two times before finally flowing to the valve downstream, resulting in larger flow resistance.

**Figure 11.** Streamline on a cross-section for the two-cage globe valve with different installation angles.

Figure 12 shows the loss coefficient in two-cage globe valves under different installation angles. It can be found that under different installation angles, the effects of the working temperature on the loss coefficient are the same with the above analysis, and the loss coefficient has the maximum value when the working temperature is 80 ◦C. Besides, with the increase in the installation angle, the loss coefficient becomes dependent on the outlet pressure even under the same working temperature, and the larger the outlet pressure, the smaller the loss coefficient. Together with Figures 11 and 12, it can also be concluded that the larger the installation angle, the larger the flow resistance, which means the smaller the cavitation intensity.

Cavitation number in two-cage globe valves with different installation angles is shown in Figure 13, it can be found that with the increase in the outlet pressure and the decrease in the working temperature, cavitation number increases, which is consistent with the above results. With the increase in the installation angle, the cavitation number also increases, which means cavitation intensity in two-cage globe valves with a large installation angle is small. The increment of cavitation number decreases with the decrease in the outlet pressure but barely changes with the variation of the working temperature. When the outlet pressure is 0.25 MPa, the increment of cavitation number is about 0.39 (the installation angle varies from 0◦ to 9◦ ), but the increment of cavitation number is only 0.144 when the outlet pressure is 0.1 MPa. Based on the results analyzed in Section 3.2, the most effective method to suppress cavitation intensity in the two-cage globe valve is by increasing the outlet pressure.

**Figure 12.** Loss coefficient in the two-cage globe valve with different installation angles.

**Figure 13.** Cavitation number in the two-cage globe valve with different installation angles.

#### **4. Conclusions**

Cavitation phenomenon and cavitation intensity in a two-cage globe valve are investigated in this study using the computational fluid dynamics method. A comparison of the water vapor volume distribution caused by cavitation for a globe valve with a different number of cages is made, and cavitation can be restricted around the valve cage and the valve body damage can be suppressed by adding a cage in the valve body, and cavitation intensity decreases with the increase in valve cage number. The effects of the valve cage with 20 slender orifices, the outlet pressure, and the working temperature on cavitation intensity and loss coefficient inside the investigated two-cage globe valve

are studied, and cavitation number and water vapor volume are chosen to stand for the cavitation intensity. With the increasing outlet pressure and the installation angle, or the decreasing working temperature, cavitation number increases which means cavitation intensity decreases. When cavitation occurs, the loss coefficient firstly decreases with the increase in the working temperature, and then rapidly increases, while the loss coefficient barely changes with varying outlet pressure but increases with the increasing installation angle. Besides, among three focused parameters, reducing the outlet pressure can remarkably decrease the cavitation intensity compared to the others and the working temperature has minimal effects on cavitation intensity.

**Author Contributions:** Conceptualization, J.-y.L. and Z.-x.G.; methodology, J.-y.L.; software, Z.-x.G.; validation, J.-y.L.; formal analysis, Z.-x.G.; resources, H.W. and Z.-j.J.; writing—original draft preparation, J.-y.L. and Z.-x.G.; writing—review and editing, Z.-x.G.; supervision, J.-y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Applied Sciences* Editorial Office E-mail: applsci@mdpi.com www.mdpi.com/journal/applsci

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-1545-8