Numerical Investigation of Tip-Vortex Cavitation Noise of an Elliptic Wing Using Coupled Eulerian-Lagrangian Approaches

Abstract

In this study, a numerical methodology is developed to investigate the tip-vortex cavitation of NACA16-020 wings and their flow noise. The numerical method consists of a sequential one-way coupled application of Eulerian and Lagrangian approaches. First, the Eulerian method based on Reynolds-averaged Navier–Stokes equation is applied to predict the single-phase flow field around the wing, with particular emphasis on capturing high-resolution tip-vortex flow structures. Subsequently, the tip-vortex flow field is regenerated by applying the Scully vortex model. Secondly, the Lagrangian approach is applied to predict the tip-vortex cavitation inception and noise of the wing. The initial nuclei are distributed upstream of the wing. The subsequent time-varying size and position of each nucleus are traced by solving spherically symmetric bubble dynamics equations for the nuclei in combination with the flow field predicted from the Eulerian approach. The acoustic pressure at the observer position is computed by modelling each bubble as a point source. The numerical results of the acoustic pressure spectrum are best matched to the measured results when the nuclei number density of freshwater is used. Finally, the current numerical method is applied to the flows of various cavitation numbers. The results reveal that the cavitation inception determined by the predicted acoustic pressure spectrum well matched the experimental result.

1. Introduction

Global economic growth has led to high demand for worldwide shipping. However, since the increase in maritime traffic requires many ships, the noise caused by the ships is rising in coastal areas and underwater ecosystems.

Murphy and King [1] reported and assessed the extent of noise exposure for residents within the vicinity of Dublin Port by long-term measurements for 45 days. Paschalidou et al. [2] also focused on the noise level in the Mediterranean port city and suggested the Strategic Noise Map and Action Plans to reduce the noise levels and exposure of the population to the noise by performing 24-h noise measurements. Bernardini et al. [3] and Fredianelli et al. [4] assessed an acoustical characterization by using short- and long-term measurements for five different small vessels that move around at various speeds in every type of port and produced noise map in Livorno’s canals. As part of the INTERREG Marittimo-Maritime projects, Bolognese et al. [5] reported the complaints by the citizens living near the ports by investigating the port noise in the North Tyrrhenian Sea. Underwater-radiated noise from ships has also adversely affected marine ecosystems by disturbing natural acoustic signals, such as that of whales, in the ocean. Bae et al. [6] reported that artificial underwater noise could harm fish hearing and their body. Sohn et al. [7] reported the adverse impact of underwater noise on marine mammals like dolphins and proposed the underwater noise management system. In this regard, developed countries, including the European Union and the United States, have already adopted specific guidelines for regulating underwater noise from ships. Underwater-radiated sound is also crucial from a military aspect, and specialized weapon systems, such as submarines, should be able to simultaneously operate at high speed and secure stealth via low noise design.

Underwater-radiated noise generated by surface and underwater vehicles, such as ships and submarines, can be subdivided into two parts based on its generation mechanism. The first part involves mechanical noises due to the structural vibration, and the second part involves flow noise generated via the flow dynamic mechanism without structural vibrations. Engines and generators are primary sources of structural vibration in vehicles, and underwater propellers are dominant sources of flow noise. At low speeds, mechanical noise typically contributes more to underwater-radiated sound than the flow noise. However, as the vehicle speed increases, the flow noise becomes more dominant. Specifically, when the speed reaches the so-called cavitation inception speed (CIS), the cavitation noise provides the highest contribution to the overall underwater noise. Thus, propeller cavitation noise is one of the most critical issues, which is considered in developing high-speed and low-noise surface or underwater vehicles [8].

Cavitation in an underwater propeller can be classified into sheet cavitation and tip-vortex cavitation, based on the structure of cavitation formation. The volumetric acceleration of cavitation produces cavitation noise as per the acoustic analogy. Given that the sheet cavitation forms and collapses periodically in association with propeller rotation, peak acoustic pressure signal due to sheet cavitation is observed at the blade passing frequency of the propeller, and thus it shows tonal characteristics [9,10]. Conversely, tip-vortex cavitation develops along the tip-vortex line, which forms helical curves that start from the propeller blade tip and extend towards the downstream direction. The underwater-radiated noise due to tip-vortex cavitation exhibits the characteristic of broadband noise due to the random behavior of the cavitation bubble [11,12]. Given that the noise levels are broadly increased in the high-frequency range when cavitation occurs, the high-speed operation of the ship and submarine is restricted. Therefore, it is important to understand the onset condition of cavitation in the propeller design process for developing low-noise propellers at high-speed operations. Primarily, tip-vortex cavitation is known to occur earlier than others such as sheet cavitation as the operation speed increases, and the accurate prediction of the tip-vortex cavitation noise is crucial for designing low-noise propellers [13,14].

Conversely, the prediction of CIS for full-scale vehicles is of considerable significance to developers and manufacturers. In practice, a small-scale laboratory test is typically performed for the analysis of tip-vortex cavitation in a design stage. Subsequently, the scaling law is employed to extrapolate the small-scale result for the full-scale application. However, there are some cases where the classical scaling law based on the McCormick formula should be corrected [15,16]. The correction is required due to real flow effects, such as nuclei distribution and dynamics and bubble/flow interaction and unsteadiness, which are generally termed as scale-effects. Based on a study by McCormick [16], several experimental studies have been conducted to extrapolate the model-scale results in full-scale results by varying the McCormick exponent [17,18].

Furthermore, a proper criterion for determining the onset of cavitation inception is not clear. Two representative measures correspond to acoustic and optical criteria. The former is the variation of sound pressure level, and the latter involves a bubble size larger than a threshold. In any case, developing a high-fidelity numerical method that can accurately and efficiently analyze cavitation initiation and noise is critical for evaluating CIS at the design stage.

The numerical methods for predicting cavitation flows can be classified into two: one is based on Eulerian approaches that the other is based on Eulerian–Lagrangian approaches. The Eulerian approaches are further classified into two: one is based on the Reynolds-averaged Navier–Stokes (RANS) solvers, and the other is based on the high-order solvers such as Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS) methods, both of which are combined with cavitation models. Kim et al. [19] investigated the effects of viscous flux in using the preconditioning and dual-time stepping methods with a homogeneous mixture model on numerical prediction of the sheet cavitation and noise of hydrofoils. Ha et al. [20,21] performed two-phase flow simulations using the Chimera grid method on massage passing interface (MPI) parallelization. Park et al. [22] used the Sauer cavitation model and volume of fluid (VOF) technique to predict tip-vortex cavitation of a NACA 66₂-415 wing. In general, the RANS simulation results [19,20,21,22] were in reasonable agreement with those that were experimentally observed when the cavitation is strong and large, such as sheet and cloud cavitation. However, it is challenging to simulate the tip vortex cavitation using the RANS solver alone because the tip-vortex core center where the cavitation bubble starts to develop belongs to the micro-scale phenomenon. It is reported that the much higher grid resolution is required to capture the vortex core [23]. Besides, a turbulence model contaminates the numerical solution of the vortex core due to excessive turbulence kinetic energy and dissipation [24]. Recently, to overcome the difficulty in resolving the tip-vortex flow, the LES and DNS simulations are employed [25,26,27], which require much numerical cost and still have difficulty for the prediction of cavitation noise. Furthermore, in the case of the Eulerian approaches, cavitation inception is determined by only using the vapor pressure. Here, it is assumed that different phases are homogeneously distributed in the given control volume. In this respect, Eulerian methods exhibit two limitations in simulating the cavitation inception. Typically, vapor pressure is assumed to be a function of only temperature, based on the assumption that the distance between vapor particles is sufficiently large. Hence, the molecular force can be ignored, which is valid for sheet cavitation. However, because the nuclei that are tiny seed particles causing tip-vortex cavitation are at a micro-scale, molecular forces cannot be ignored [28]. Given that the Eulerian method uses the volume fraction to quantify the cavitation flow, it is difficult to distinguish the interface between liquid and vapor precisely. Thus, it is difficult to determine the exact location and size of the cavitation.

A Lagrangian approach based on bubble dynamics can be applied to evaluate the cavitation occurrence with more direct physical quantities. The bubble dynamics model aids in determining the time-varying position and size of the bubble that develops from each nucleus, and thus enables accurate prediction of the noise due to the bubble. Hsiao et al. [29] predicted the cavitation inception in a line vortex flow using a Lagrangian approach combined with the single-phase RANS solvers. It was shown that that the size of nuclei plays an essential role in extrapolating the small-scale results to full-scale results.

In this study, a systematic numerical methodology is developed to investigate the tip-vortex cavitation inception and noise of an elliptic NACA16-020 wing. The numerical method consists of sequential one-way coupled applications of Eulerian and Lagrangian approaches. First, the Eulerian method is applied to predict the hydrodynamic flow field around the wing with a focus on accurate simulation of the tip-vortex flow. The RANS solver is used to predict the single-phase hydrodynamic flow around the wing. Subsequently, the tip-vortex flow field is regenerated by employing the Scully vortex model. The model parameters, such as the circulation and core-radius, are determined using the information from the predicted hydrodynamic flow. Secondly, the Lagrangian approach is applied to predict the tip-vortex cavitation inception and noise of the wing. The initial nuclei are distributed upstream. The subsequent time-varying volume and location of the nucleus are computed by solving bubble dynamics equations in combination with the hydrodynamic flow field of the wing, especially the regenerated tip-vortex flow. Finally, the compressible acoustic pressure at a selected monitoring position is computed by modelling each bubble as a point monopole source. The validity of the present numerical methodology is confirmed by comparing its prediction results with the measured results.

The overall procedure of the current numerical method is similar to that reported by Hsiao et al. [29]. However, in the present study, the numerical methods used in each numerical step are investigated in more detail. Thus, the four main contributions of the present study are as follows. First, the effects of turbulence models and grid resolutions on the prediction of the tip-vortex flow are quantitatively compared, which can provide a guideline for performing an accurate numerical simulation of tip-vortex flow. Second, the systematic application of the vortex model is presented for regenerating high-resolved tip-vortex flow from the background flow obtained from the RANS solver. It includes the methods for tracing the tip-vortex core location and determining the vortex strength and core-radius of the vortex model. Third, the bubble dynamics equations are investigated in detail with a focus on the effects of each term of the equations on prediction accuracy. Fourth, the critical physical issues in predicting acoustic pressure due to tip-vortex cavitation are examined in detail. Notably, the effects of the number density of initial nuclei on the tip-vortex cavitation inception and noise are quantitatively analyzed. The results indicate the merit and constraints of the current numerical methods and their suitability as an analysis tool for investigating tip-vortex cavitation.

The structure of the present paper is as follows. In Section 2, the numerical methods are presented to describe the RANS simulation, vortex model, initial nuclei distribution, bubble dynamics equations, and acoustic analogy that are used in the subsequent numerical simulation. In Section 3, the experimental setup is introduced to visualize the cavitating tip-vortex of the wing and to measure the hydro-acoustic pressure due to cavitation. In Section 4, the numerical results are presented. First, the results of the grid-refinement study are presented, and then the effects of the turbulence model on the resolution of predicted tip-vortex flow are shown. Subsequently, the reconstruction of tip-vortex flow with the vortex model is described in detail. Finally, the simulation of tip-vortex cavitation via the bubble dynamics equations and prediction of acoustic pressure due to the cavitation are presented. The prediction results of cavitation tip-vortex formation and that of the hydro-acoustic pressure due to the cavitation are compared with the experimental results. The effects of the number density of initial nuclei on the prediction results are also investigated.

2. Numerical Methods

In this section, the numerical methods presented in this study are described in detail. The sequential one-way coupled Eulerian–Lagrangian method is applied to simulate the tip-vortex cavitation inception and noise of a NACA16-020 blade. It consists of six steps, as shown in Figure 1.

In the first stage, the hydrodynamic steady flow field around the wing is obtained using a single-phase RANS solver. However, the accuracy of the numerical solution in resolving tip-vortex flow depends on the dissipation of numerical schemes, grid-resolution, and turbulence models. Therefore, a grid-refinement study is performed using four types of grids. Subsequently, the influence of turbulence models is quantified by comparing the tip-vortex structure predicted using each turbulence model.

The vortex formed at the wing-tip is physically dependent on the Reynolds number and is convected in a long line-shape along with the downstream. Typically, the vortex core is located on the lowest pressure region, such that the wing-tip vortex cavitation occurs along the vortex core line, starting from the wing-tip and stretching to downstream. Therefore, the center of the vortex core is identified in the flow field obtained from the RANS result and used as the central axis of the vortex regenerated using the vortex model.

From a macroscopic point of view, cavitation is a phase change phenomenon from liquid to vapor that occurs at a pressure lower than vapor pressure. However, from a microscopic view, nuclei grow and form cavitation bubbles when the ambient pressure is lower than the critical pressure. Nuclei distribution depends on water quality, which varies according to temperature, seawater or freshwater, seasons, and regions.

Experimental studies were performed to find the relationship between water quality and nuclei distribution. Vanning et al. [30] investigated water susceptibility and background nuclei content in a water tunnel using a cavitation susceptibility meter and showed that the cumulative histogram of nuclei concentration against critical pressure followed a power-law dependence over a large range of concentrations and pressures. O’ Hern et al. [31] measured the distribution of the number concentration density of nuclei in the Pacific Ocean and suggested the number density as a function of the initial nuclei radius. In this study, initial nuclei distributions are determined using this data as the standard seawater condition.

The initial nuclei are distributed upstream of the wing. The subsequent behavior of nuclei (or cavitation bubble), such as bubble radius variation and bubble trajectory in the flow field, is predicted by solving the bubble dynamics equations. A single bubble is assumed as spherically symmetric and is treated as a point monopole source. Subsequently, the flow noise due to cavitation is predicted using the acceleration in bubble volume, which is a function of the radius of a spherical bubble. A detailed description of each step of numerical methods is provided in the following sections.

2.1. Background Hydrodynamic Flow Simulations Using RANS Solver

In this study, an elliptic wing with NACA 16-020 cross-section [32] is utilized for the investigation of wing-tip vortex cavitation and noise. The wing of the same geometry is also experimentally investigated. The detailed shape of the targeted wing, dimensions of the test section of the water tunnel, and applied boundary conditions are depicted in Figure 2. The chord and height of the wing are 0.5 m, and the angle of attack is set at 8°. The dimensions of the entire computational domain are determined based on the actual size of the test section of the large cavitation tunnel (LCT) in Korea Research Institute of Ships and Ocean Engineering (KRISO). The related boundary conditions are set as identical to those used in a previous study [18]. The boundary conditions of inlet velocity and outlet pressure are set at 9 m/s and 49.885 kPa, respectively. The vapor pressure is 1.304 kPa, which corresponds to a cavitation number of 1.2.

The hydrodynamic flow field around the wing is predicted by numerically solving the three-dimensional steady incompressible RANS equations. The semi implicit method for pressure linked equations (SIMPLE) scheme, where the continuity and momentum equations are sequentially solved, is used to satisfy the mass conservation. The gravity effect is applied along the z-coordinate direction. Fluent v19.2 tool is used to perform the numerical RANS simulations described above.

A wing-tip vortex begins to form in the lowest pressure region around the wing-tip and develops in the downstream direction. It becomes partially laminar inside its vortex core due to the extremely low turbulent kinetic energy when it is fully established [22]. However, it is difficult to numerically resolve the vortex structure due to the excessive turbulent viscosity of the turbulence model. In this study, five turbulence models, namely, standard k-ε, realizable k-ε, re-normalization group (RNG) k-ε, k-ω SST, and the Reynolds stress model (RSM), are considered. The second-order accuracy is retained for numerical differencing schemes, and thus a grid-refinement study is performed to find an appropriate grid size near the vortex structure. Four types of grids are used, and their detailed information is listed in

Table 1. Grid refinement in vortex domain.

Grid Refinement
Domain	Maximum Cell Size (m)	Number of Cells (Vortex Domain Only)	Total Number of Cells (Tunnel + Wing-Near + Vortex)
Tunnel	~0.1	-	-
Wing near	~0.0028	-	-
Vortex	~0.0052 (Test)	8 M *	21.7 M
	~0.0014 (Coarse)	3 M	22.1 M
	~0.0007 (Middle)	14.5 M	47.8 M
	~0.00035 (Dense)	83.2 M	116.5 M

* M: Million.

. The computational domains are divided into three regions: tunnel, wing-near, and vortex zones, as shown in Figure 3. Furthermore, the identical grid-resolution is retained in the tunnel and wing-near zones for all the four grid types. However, the resolution in the vortex zone is varied. The dependence of numerical solutions on grid-resolution is analyzed in terms of the variation of vortex structure along the downstream direction, especially the tangential velocity distribution around the vortex core and vortex core-radius.

To identify the location of the vortex core line, the preliminary RANS simulation is performed with the test grids. The computational area for the simulation is shown on the left side of Figure 3. After the vortex core line is identified from this simulation result, the vortex zone is newly modelled to reduce the numerical cost, as shown on the right side of Figure 3. The computational domain, shown on the right side of Figure 3, is adopted for the grid-refinement study using the so-called coarse, middle, and dense grids. Furthermore, detailed information on the grids is provided in Table 1. The y⁺ distribution on the wing surface is shown in Figure 4a. The grids of y⁺ ≈ 20~30 are used in the vicinity of the leading and trailing edges, and the grids of y⁺ < 110 are used on the other surface regions. Park [33] estimated the free surface shape and hydrodynamic drag of a ship using the grids of the resolutions between 50 < y⁺ < 150. The predicted results showed good agreements with the measured ones. The unstructured grid distribution is presented in Figure 4b. The tetrahedron grids are used in the Tunnel domain, and the hexahedron grids are used in the wing-near and vortex domains.

2.2. Tip-Vortex Regeneration Using Vortex Model

It is, generally, difficult to capture the vortex structure using the RANS simulation where the vortex tends to be over-dissipated due to numerical viscosity of numerical difference schemes and excessive turbulence viscosity of the turbulence model. To resolve this issue, the tip-vortex structure is regenerated by applying a vortex model to the background flow field obtained from the RANS simulation.

The position of the vortex core should be identified in the RANS result before the vortex model is applied. The location of the vortex core is determined using the λ₂ criterion, which is the second-largest eigenvalue of the matrix. It is the sum of the strain tensor (S) and vorticity tensor (Ω) of the flow, as shown in the following equations for the y-z plane.

$$\nabla \overrightarrow{u}=\left(\begin{array}{cc}\frac{\partial v}{\partial y}& \frac{\partial v}{\partial z}\\ \frac{\partial w}{\partial y}& \frac{\partial w}{\partial z}\end{array}\right)=S+\Omega$$

(1)

$${\lambda }_2=\frac{{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial z}\right)}^2}{2}+\frac{\partial v}{\partial z}\frac{\partial w}{\partial y}$$

(2)

Here, y and z denote Cartesian coordinates on the cross-sectional plane perpendicular to the mean flow direction, and v and w denote velocity components in the y- and z-directions, respectively. The λ₂ criterion is widely used to identify coherent vortex structure because it can accurately find the coherent vortex core for various flow fields. The region of negative λ₂ is defined as the vortex core, and the minimum corresponds to the central axis of the vortex core. It is possible to define a vortex core on a cross-section plane perpendicular to the tip-vortex vector [34]. In this study, λ₂ is defined in the two-dimensional plane (y-z plane) in Equation (2) because the vortex is formed on a cross-sectional plane normal to the mean flow direction (x-dir.).

The vortex model is applied on each plane normal to the vortex core line determined from the background flow field. Hsiao et al. [35] generated vortex flow fields using the Rankine vortex model. They determined the vortex properties of the circulation, Г, and vortex core-radius, a_c, using the reference length and freestream velocity based on the classical thin wing theory and assuming a fully turbulent boundary layer. However, this method can overestimate the vortex strength. Therefore, in this study, the circulation is estimated using the vortex core-radius and maximum tangential velocity, $u_{\theta }$, in the region near downstream of the wing-tip, where the vortex begins to form. The vortex flow field is regenerated using the Scully vortex model. The Rankine vortex model, which assumes a solid body rotation due to the viscous effect inside the vortex core and a free (potential) vortex outside, exhibits a discontinuity at the boundary of the vortex core from the fixed vortex to free vortex. The Scully vortex model used an algebraic velocity distribution to overcome this singularity problem. The tangential velocity in the Scully vortex model is defined as a function of r as follows:

$$V_{\theta }\left(r\right)=\frac{\Gamma }{2\pi }\left(\frac{r}{a_c^2+r^2}\right)$$

(3)

Equation (3) is used to determine the circulation from the tangential velocity and the vortex core-radius. Subsequently, the pressure distribution is determined from the three-dimensional steady-state Navier–Stokes equation for the radial direction as follows:

$$\begin{array}{ll}{u}_r\frac{\partial u_r}{\partial r}+\frac{u_{\theta }}{r}\frac{\partial u_r}{\partial \theta }& +w\frac{\partial u_r}{\partial z}-\frac{u_{\theta }^2}{r}\\ & =-\frac{1}{\rho }\frac{\partial p}{\partial r}+\nu \left(\frac{{\partial }^2{u}_r}{\partial r^2}+\frac{1}{r^2}\frac{{\partial }^2{u}_r}{\partial {\theta }^2}+\frac{{\partial }^2{u}_r}{\partial z^2}+\frac{1}{r}\frac{\partial u_r}{\partial r}-\frac{2}{r^2}\frac{\partial u_{\theta }}{\partial \theta }-\frac{u_r}{r^2}\right)\end{array}$$

(4)

Typically, the radial velocity, u_r, is negligible when compared to the tangential velocity $u_{\theta }$ of vortex core [36]. Therefore, by ignoring the radial direction component and assuming no velocity variation in the tangential direction, i.e., $\frac{\partial u_{\theta }}{\partial \theta }=0$, a more straightforward relationship between the tangential velocity and pressure [37,38] is obtained from Equation (4) as follows:

$$\frac{u_{\theta }^2}{r}=\frac{1}{\rho }\frac{\partial p}{\partial r}$$

(5)

where, subscript r, θ, and z denote the radial, tangential, and axial directions, respectively, with its origin on the center of the vortex core. From Equation (5), pressure distribution, p (r), in the Scully vortex model can be derived as follows:

$$p\left(r\right)=p_{\infty }-\frac{\rho {\mathsf{\Gamma }}^2}{8{\pi }^2}\frac{1}{r^2+a_c^2}$$

(6)

The pressure field around the vortex core is reconstructed by combining Equation (6) with the background RANS field. Azimuthal velocity components in the RANS solution are replaced with that of the Scully model, while the axial velocity components in the RANS solutions are retained.

The vortex structure generated on the wing-tip is physically smeared with the thickening of the core-radius via viscous diffusion, as it is convected downstream. Moore and Saffman [39] derived the formula for predicting the variation in the tangential and axial velocities of the tip-vortex during its convection in the downstream direction using the lifting-line theory under the assumption of inviscid flow. The singularity that occurs on the vortex center in the inviscid flow field is treated by considering the influence of viscous stress inside the vortex core. The vortex core-radius also increases due to the viscous effects, as the vortex moves toward the downstream direction. The variation of the vortex core-radius due to the viscous effects is considered using the following equation [39].

$$\frac{\mathsf{\Delta }{a}_c}{c}=2.92{\left(\frac{x}{c}\right)}^{0.5}\frac{1}{R{e}_c^{0.5}}$$

(7)

where, $\mathsf{\Delta }{a}_c$ denotes the variation in vortex core-radius, c denotes the chord length, x denotes the distance between wing-tip and the concerned position downstream, and Re_c denotes Reynolds number dependent on the chord length.

2.3. Initial Nuclei Distribution and Critical Pressure

Nuclei distributed in an undissolved state of water grow into a bubble, i.e., cavitation occurs when the pressure around nuclei is lower than the critical pressure of nuclei [29,31]. The pressure value below which nuclei begin to grow into bubbles is termed as the critical pressure. The nuclei distribution and critical pressure vary with nucleus size and are the main properties of the water quality [17].

Nuclei distribution in the inlet flow should be determined for the prediction of cavitation inception. Nuclei distribution varies with water type (seawater or freshwater), sea area, water temperature, and seasons. Cavitation is more likely to occur in seawater than in freshwater because seawater includes more impurities, such as plankton, gas bubble, and dust, and has more nuclei than freshwater. O ‘Hern et al. [31] measured the nuclei distribution in different areas of Pacific seawater. H. Kamiirisa [40] compared the nuclei distribution between seawater and freshwater and reported that the seawater environment could be approximated by changing the air content rate in the freshwater. It was reported that the number of nuclei in seawater is about ten times higher than that in freshwater. The nuclei distribution of the freshwater used in the experiment is reproduced by multiplying the reported seawater condition by 0.1. The number density of initial nuclei used for the simulations is shown in Figure 5, which is a modified version of the number concentration density distribution of nuclei reported by O’Hern et al. [31]. O’Hern et al. showed that the holographic detection method was a more reliable technique. Besides, the nuclei distribution for the range of 10~25 µm was measured by using a sampled water of small volume, 1 mL, while the nuclei larger than 25 µm were measured by using a sampled water of large volume, 150 mL. Therefore, the nuclei larger than 25 µm are used in the present study. A similar range of initial nuclei sizes was used in the previous studies of Park et al. [37], Hasio et al. [35,41], and Kamiirisa [40]. It should be noted that a parametric study with the varying number density of nuclei is also conducted to investigate the effects of the initial nuclei distribution on the predicted cavitation inception and noise.

When the pressure around nuclei is low enough to reach critical pressure, the nuclei grow into bubbles, and cavitation occurs. Brennen [28] described this behavior of nuclei in terms of intermolecular forces from a microscopic point of view.

Figure 6 shows the variation in potential Φ of intermolecular force based on the distance between molecules, x. Furthermore, slope $\frac{\partial \mathsf{\Phi }}{\partial \mathrm{x}}$ is the intermolecular force. Intermolecular force is attractive at distances where the slope is positive, while it is repulsive at distances where the slope is negative. The so-called surface tension originates from the intermolecular force. Subsequently, position $x_0$ represents an equilibrium state where there is no intermolecular force, i.e., no surface tension. The initial nuclei are assumed to be in an equilibrium where external pressure is equal to the internal pressure of the nuclei. When the external pressure decreases, the volume of the nuclei increases, and intermolecular distance on the nuclei surface increases. Hence, the attraction force is applied. Therefore, for nuclei to grow into cavitation bubbles, the external pressure must be sufficiently lowered based on the difference between vapor pressure and additional force due to surface tension. Typically, the equilibrium distance $x_0$ of the liquid is extremely small, approximately 10⁻¹⁰ m. If the cavitation shape is relatively large, similar to that of sheet cavitation, then the effect of surface tension can be neglected. In this case, only the vapor pressure should be considered because the distance between molecules is large enough, and thus it is in the region of the straight line shown in Figure 6. However, to predict the tip-vortex cavitation, the behavior of nuclei in the micro-scale is of significant importance, and thus the critical pressure must be considered [42]. The critical pressure varies according to nuclei size, and the degree of cavitation is closely related to the nuclei distribution. The critical pressure is defined as follows:

$$p_{cr}=p_v-\frac{4S}{3R_0}{\left[3\left(1+\frac{p_0-p_v}{2S/R_0}\right)\right]}^{-1/2}$$

(8)

It can be inferred from Equation (8) that when the initial nuclear radius is reduced, the critical pressure is lowered. Hence, small nuclei are less likely to become cavitation bubbles than large nuclei.

Table 2. Critical pressure according to initial nuclei size.

R0 (µm)	30	40	50	60	70	80	90	100
pcr (Pa)	725	924	1030	1094	1137	1167	1189	1206

listed the critical pressure of nuclei according to its initial size.

2.4. Spherical Bubble Dynamics

Among various types of cavitations observed in an underwater propeller, wing-tip vortex cavitation occurs first as the operating speed increases. Therefore, to design low-noise and high-performance underwater propellers, the process of creation of wing-tip vortex cavitation should be closely examined. As described in the introduction, the RANS solvers with a homogeneous mixture model have been frequently used for the cavitation simulations of hydrofoils and propellers. Specifically, many studies have been conducted on sheet cavitation and tip-vortex cavitation, and the predicted results exhibit qualitative agreements with the measured data in terms of the overall shape of cavitation [19,20,21,22]. To define the CIS, Gaggero, et al. [43] computed the minimum pressure and volume of tip-vortex cavitation and visually compared the cavitation tendency with the experiment. However, it was reported that the cavitation tendency differs between the experimental and numerical results, and the difference increases when the body shape is complicated, such as ducted propellers. Furthermore, Raju et al. [44] suggested that although the homogeneous mixture model accurately predicts the macroscopic behavior of bubbles, local oscillations within macroscopic movements can be accurately predicted using the bubble dynamics.

Although there are no generally accepted criteria to define the CIS, two approaches are frequently used. The first approach involves a visual approach where the size of bubbles is visually observed, and the second approach involves an acoustic approach where the sudden increase in noise level is used as an indication of the cavitation inception [29]. Therefore, in this study, wing-tip vortex cavitation is simulated using the Lagrangian approach, where the bubble dynamics model is used under the assumption of the spherically symmetric bubble. The advantage of using the bubble dynamics model is that it can effectively track the volume variation of bubbles and the noise due to the volumetric bubble acceleration, which can be used as a visual and acoustic index of cavitation inception, respectively.

The volume variation of the bubble can be predicted using the Keller–Herring equation, where the bubble radius is determined with the given external information such as external pressure outside the bubble. The location of the bubble can be traced by solving the trajectory equation with the given flow information around the bubble. The Keller–Herring equation includes terms that reflect liquid compressibility and consider energy dissipation during a volume change, while the Rayleigh–Plesset equation maintains the energy as constant under the assumption of incompressible flow [28,45]. The Rayleigh–Plesset equation with the incompressibility assumption is as follows:

$$R\ddot{R}+\frac{3}{2}{\dot{R}}^2=\frac{1}{\rho }\left[p_v+p_{g0}{\left(\frac{R_0}{R}\right)}^{3k}-p-\frac{2\gamma }{R}-\frac{4\mu }{R}\dot{R}\right]+\frac{{\left(\overline{U}-{\overline{U}}_b\right)}^2}{4}$$

(9)

The Keller–Herring equation with the compressibility effects is as follows:

$$\left(1-\frac{\dot{R}}{c}\right)R\ddot{R}+\frac{3}{2}\left(1-\frac{1}{3}\frac{\dot{R}}{c}\right){\dot{R}}^2=\left(1+\frac{\dot{R}}{c}+\frac{R}{c}\frac{d}{dt}\right)\left[h-\frac{p_v}{\rho }\right]+\frac{{\left(\overline{U}-{\overline{U}}_b\right)}^2}{4}$$

(10)

where $h=\left(p_B-p\right)/\rho , p_B=p_g-\frac{1}{R}\left(2\gamma +4\mu \dot{R}\right)$, R denotes the bubble radius, p_v denotes vapor pressure, p_g0 denotes the initial internal gas pressure in the nuclei, R₀ denotes initial nuclei radius, γ denotes surface tension, μ denotes liquid dynamic viscosity, p denotes fluid pressure, and $\overline{U}$ denotes the velocity on the bubble surface. Furthermore, ${\overline{U}}_b$ denotes the moving velocity of a bubble, and the superscript ‘∙’ denotes the time derivative term. Additionally, c, h, and k denote liquid sound speed, enthalpy, and polytropic gas constant, respectively. If the speed of sound becomes infinite, then Equation (10) is reduced to Equation (9). In this study, the Keller–Herring equation is used for predicting the variation in the volume of the bubble. However, in classical bubble dynamics, the external pressure on the bubble surface is extrapolated from the pressure at the position, which is the bubble center if there is no bubble. This approximation neglects the variation of the bubble surface pressure due to the change in bubble radius. Therefore, the external pressure in Equation (10) is determined using the surface-averaged pressure method [29,41], where the external properties are computed by averaging the pressure and velocities in the flow field on the bubble surface.

Gas in the nuclei or bubble is assumed to be an ideal gas. However, there is a significant difference in the time-scale between the growth and decrease phases of the bubble: the time-scale in the bubble growth phase is much higher than that in the bubble decrease phase. Therefore, the bubble growth process can be assumed as isothermal, and the bubble decrease process can be assumed as adiabatic. The polytropic exponent for isothermal and adiabatic processes is set as [32] follows:

$$polytropicexponent,k=\{\begin{array}{cc}1& Isothermalprocess\left(\dot{R}>0\right)\\ 1.4& Adiabaticprocess\left(\dot{R}<0\right)\end{array}$$

Nuclei containing non-condensable gas and liquid–vapor are in a dynamic equilibrium state with the surrounding fluid. The pressure due to non-condensable gas, p_g, and the pressure due to vapor, p_v, equilibrates with the surface tension on the nuclei surface and the pressure due to external liquid, p_L, as shown in Figure 7.

This equilibrium relationship can be written as follows:

$$p_L=p_v+p_g-\frac{2\gamma }{R}\hspace{1em}\hspace{1em}where{p}_g=p_{g0}{\left(\frac{R_0}{R}\right)}^{3k}$$

(11)

Here, the initial gas pressure p_g0 is determined from initial flow field information and initial nuclei radius [28]. Given that the bubble surface variation occurs very fast when compared with the time-scale of bubble dynamics simulation, the vapor can be assumed to be in equilibrium at each moment, and the vapor pressure remains constant. Conversely, the expansion of the gas occurs more slowly when compared with the time-scale of the bubble dynamics simulation. Therefore, it is necessary to consider the change in the internal gas pressure by changing the bubble radius at every moment. The gas pressure in Equation (11) is varied by changing the bubble radius.

The bubble location can be traced using the bubble trajectory equation, which is as [42,46] follows:

$$\frac{d{\overline{U}}_b}{dt}=-\frac{3}{{\rho }_l}\nabla p-2g+\frac{3}{4}\frac{C_D}{R}\left(\overline{U}-{\overline{U}}_b\right)\left|\overline{U}-{\overline{U}}_b\right|+\frac{3}{2\pi }\frac{C_L}{R}\sqrt{\frac{{\mu }_l}{{\rho }_l}}\frac{\left(\overline{U}-{\overline{U}}_b\right)\times \overrightarrow{\omega }}{\sqrt{\left|\overrightarrow{\omega }\right|}}+\frac{3}{R}\left(\overline{U}-{\overline{U}}_b\right)\dot{R}$$

(12)

where subscripts l and b denote the liquid and bubble property, respectively, and the terms without the subscript denote the external features on the bubble surface. The first term in the right-hand side of Equation (12) presents the effects of the pressure gradient on the bubble surface, and the second term denotes the buoyancy effect. The third and fourth terms denote the drag and lift acting on the bubble, respectively. The final term is the Bjerknes force term, which reflects the coupling effects between bubble volume change and bubble motion. Furthermore, C_D, which can be derived from the empirical formula of Haberman and Morton [47], is the drag coefficient of the bubble due to the surrounding fluid. Specifically, C_L is the bubble lift coefficient given in [48,49].

2.5. Acoustic Pressure Due to Cavitation Bubble

The bubble is assumed to be spherically symmetric during its life, and each bubble is considered as a monopole point source. Fitzpatrick and Strasberg [50] proposed a bubble noise model where the noise emitted from a bubble is computed using the volume acceleration of the bubble. This model is consistent with that derived using the Lighthill’s acoustic analogy [51]. Flow noise sources are classified into monopole, dipole, and quadrupole sources in the acoustic analogy. The strength of the monopole source is proportional to volumetric acceleration. An acoustic pressure signal at the receiver is computed using the source information at the retarded time as follows:

$$p_a\left(x, t\right)=\frac{R\rho }{l}\left[R\ddot{R}\left(y, \tau \right)+2{\dot{R}}^2\left(y, \tau \right)\right],\tau =t-\frac{\left|x-y\right|}{c}$$

(13)

where, x and y denote the position vector of the receiver point and source location, respectively, and τ and t denote the retarded time (or source time) and observer time, respectively.

3. Experimental Methods

To validate the numerical results, the corresponding experiments are performed. The experiments on cavitation inception and noise of the wings are conducted in the LCT at KRISO to provide the comparison data for the prediction results. The overall dimensions of the LCT are 60.0 m (L) $\times$ 6.5 m (B) $\times$ 19.8 m (H), and those of the test section are 12.5 m (L) $\times$ 2.8 m (B) $\times$ 1.8 m (H), as shown in Figure 8. More detailed specifications of the LCT can be found in [18,52]. The experimental setup is shown in Figure 9. A wing with the chord length and span both of which are 500 mm and a maximum thickness of 100 mm is installed at the centerline of the ceiling of the test section, which is the same as the numerical setup shown in Figure 2. A charge-coupled device (CCD) video camera with the accompanying lighting system is set outside the test section to visualize the tip-vortex cavitation. The inflow velocity in the test section is set to 9.0 m/s, and the cavitation number is controlled by changing the tunnel pressure. Acoustic pressure is measured using a B&K Type 8105 hydrophone located under the floor of the test section with a sampling rate of 256 kHz.

4. Numerical Results

As described in Section 2, the current numerical method consists of the sequential applications of the Eulerian and Lagrangian approaches. In Section 4.1 and Section 4.2, the effects of grid resolution and turbulence models on the tip-vortex flow predicted using the RANS solver are presented. In Section 4.3, the reconstruction of tip-vortex flow is detailed. In Section 4.4, the initial nuclei distribution used for the simulations is explained. Finally, the tip-vortex cavitation inception and noise are predicted and compared with the measured data in Section 4.5.

4.1. Effects of Grid Resolutions

First, a grid-refinement study is conducted. Figure 10 shows the numerical results obtained from the RANS simulations using the four types of grids listed in Table 1. The Reynolds stress model is used as the turbulence model. The results are presented in terms of turbulent kinetic energy and vorticity magnitude.

Given that the roll-up process of vortex formation involves a small-scale vortex structure and high gradient change near the vortex core, the grid-resolution is of critical importance to the accuracy of the numerical solutions.

Table 3. The number of grid points across the vortex core to predict tip-vortex flow.

Number of Grid Points across the Vortex Core
Test Grids	Coarse Grids	Middle Grids	Dense Grids
3.6	13.5	27	54

compares the grid-resolution of different grid types in terms of the number of grids that are used to resolve the vortex core.

The vortex core diameter, 0.0189 m, is computed using the distance between peak-to-peak of the tangential velocity profile obtained from the simulation result using the coarse grids. In the RANS simulation, the wing-tip vortex structure is dissipated rapidly along the downstream direction as the grid size increases. Hence, as the number of grids in the vortex region increases, the sharpness of the resolved vortex structure increases. Specifically, the vortex generated at the wing-tip is known to exhibit a partial laminar flow region in the vortex core due to the low rate of turbulence energy production. This fact can be confirmed in Figure 10 by identifying the vortex core region, where turbulent kinetic energy is low.

Figure 11 compares the non-dimensional profiles of tangential velocity and pressure across the vortex core region at the position 0.5c downstream of the wing-tip. It is inferred that as the grid-resolution becomes dense, the predicted profiles of tangential velocity and pressure become steeper and converge to a single curve.

4.2. Effects of Turbulence Models

Figure 12 compares the distribution of turbulence kinetic energy in the tip-vortex obtained from the RANS results using five different turbulence models on the coarse grids. Evidently, the vortex structure cannot be captured from the numerical result obtained using the standard k-ε model. It is important to note that the numerical result does not converge for the standard k-ε model. The detached eddy on the wing tip surface is critical to the accurate simulation of the tip vortex. However, the standard k-ε model, which was initially developed for the eddy simulation in the free field, tends to over-estimate the turbulence kinetic energy near the wall surface [24]. This limitation of the standard k-ε model seems to degrade the convergence of the RANS simulation.

In the cases of k-ω SST and realizable k-ε models, the vortex structure can be captured. However, it is contaminated due to excessive turbulent energy during its convection. Given that the RNG k-ε model has been widely used for simulating rotating flows, it resolves the vortex structure in a better manner than the preceding three models. However, it still overestimates the turbulent kinetic energy. The best result is obtained from the simulation using the Reynolds stress model. The vortex structure is resolved extremely well, and the partial laminar region that occurs in the vortex core is also reproduced well [22].

Figure 13 compares the distributions of non-dimensional tangential velocity and pressure coefficient on the line passing through the vortex core in the plane 0.5c downstream from the wing-tip. All the cases except for the standard k-ε model show similar results. A closer examination of the detailed distributions reveals that the result obtained using the Reynolds stress model captures the vortex core at the highest sharpness. This is due to the lower turbulent kinetic energy along the downstream direction without noticeable dissipation, as shown in Figure 12.

4.3. Vortex Core Identification and Vortex Modelling

Although the simulation result using the Reynolds stress model shows the maximum resolved flow field of the tip-vortex, the predicted vortex structure becomes weaker during its convection along the downstream direction due to numerical damping. To deal with this issue, the vortex model is used to reproduce the tip-vortex with high-resolution. First, the λ₂ criteria are used to identify the vortex core in the flow field. Figure 14 shows the line constructed by tracing the minimum value of λ₂ on each plane normal to the mean flow. The constructed line matches well with the vortex core line. Figure 15 compares the vortex core trajectories predicted in this manner for the cases with different grids and turbulence models except for the standard k-ε models. Hence, there is no noticeable variation in the identified vortex core lines among considered the cases, and the numerical results well matched the experimental result.

Figure 16 shows the pressure coefficient variation along the vortex core trajectory predicted using different grids and turbulence models. As the grids become coarser, the pressure increases. The pressure difference between the different grids also increases as the vortex is convected downstream. These results imply higher vortex dissipation as the grid size increases. It is observed that the pressure predicted using the turbulence models, except for the Reynolds stress model, increases as the vortex is convected in the downstream direction. Furthermore, the pressure increment is maximum and comparable for the cases of the realizable k-ε and the k-ω SST models. Furthermore, it is relatively lower for the RNG k-ε model. Given that the standard k-ε model does not accurately simulate the vortex structure and even cannot converge, as shown in Figure 13, it is not considered in this comparison. The pressure distribution predicted using the Reynolds stress model is minimum and appears to stay stable in the downstream direction. It can be inferred from the results in Figure 15 and Figure 16 that the trajectory of predicted vortex core is not sensitive to the grid-resolution and turbulence models. However, the strength of the vortex represented by the minimum pressure is susceptible to grid-resolution and turbulence models. Therefore, in this study, all the following simulations are conducted using the Reynolds stress model with a coarse grid. This ensures the accuracy of numerical solutions at a moderate numerical cost.

Figure 17 shows the variation in pressure and turbulent kinetic energy along the vortex core trajectory. Theoretically, the vortex is a rotating flow about the vortex center on which the rotation speed is zero. Therefore, within the vortex core, there is a region where the flow is partially laminar. As shown in Figure 17, this region begins to form between sections (e) to (f) and develops into a complete vortex structure. The origin of tip-vortex begins to develop from the wing-tip. Although the region with the minimum pressure occurs at the wing-tip, the pressure gradually increases as it moves from (a) to (c) and then decreases from (c). The pressure increase between these planes seems to be due to the adverse pressure gradient along with the chord-direction of wing airfoil sections near the wing-tip. Typically, the wing-tip vortex cavitation begins to form in the location slightly downstream from the wing-tip and gradually approaches the tip from the back in the early stage of cavitation. This phenomenon is related to the distribution of lower pressure region formed downstream of the wing-tip along the vortex core trajectory.

Given that the wing-tip vortex is dissipated by the numerical damping, as it is convected downstream, the tip-vortex is reproduced by applying the Scully vortex model on the plane perpendicular to the vortex core trajectory. To regenerate the vortex structure, the initial vortex core-radius and tangential velocity of the Scully vortex model should be determined. The initial vortex core-radius is defined as half the distance between the points of maximum and minimum velocities on the plane perpendicular to the vortex core trajectory. The distance can be easily estimated using the tangential velocity distribution, as shown in Figure 11 and Figure 13. The initial tangential velocity can also be calculated using the same approach. Along the downstream direction, the physical vortex dissipation due to flow viscosity is predicted using Equation (7), which provides an estimation for the increase in the vortex core-radius.

Figure 18 compares the variation of pressure coefficient along the vortex core trajectory and the distribution of z-directional velocity in the tip-vortex flow between the RANS and vortex-regenerated flow fields. The strength of the vortex becomes weaker as it is convected downstream in the RANS results. Conversely, the vortex structure is maintained as almost constant in the flow field regenerated using the vortex model, which is more similar to the experimental observation [32].

4.4. Initial Nuclei Distribution

Figure 19 illustrates the initial distribution of nuclei of various sizes in an arbitrary cuboid volume. The initial nuclei are homogeneously distributed in an arbitrary volume upstream of the wing-tip [37]. The number density of nuclei based on size was determined from the number density shown in Figure 5. The total number of initial nuclei is 1721. As shown on the right side of Figure 19, nuclei corresponding to each size are evenly distributed within an arbitrary volume under the assumption of homogeneous distribution in the flow field. Figure 20 shows the initial distribution of nuclei in the computational domain. Furthermore, the nuclei are located upstream of the wing, and their subsequent volume and location are computed by solving Equations (10) and (12) during their convection via the regenerated vortex flow field embedded in the background RANS result.

4.5. Bubble Dynamics and Noise

Figure 21 compares the predicted tip-vortex cavitation with the measured tip-vortex cavitation. Figure 21a shows the snapshot of tip-vortex cavitation, which is recorded using the CCD video camera, as explained in Section 3. Figure 21b–i show the predicted time-history of bubble behavior in the flow field. The black dashed line denotes the vortex core trajectory that is extracted from the experiment, and the blue sphere denotes the bubble. The total number of nuclei considered in this study is 10,326 and the time interval for bubble simulation is 0.5 µs. Among the nuclei convected through the regenerated flow field, a few grow into cavitation bubbles as they pass through the vortex core, which is one of the lower pressure regions. The bubbles move downstream along the vortex core trajectory and repeat the compression and expansion of the volume in response to changes in the surrounding pressure.

Figure 22 shows the time-history of the bubble radius for the nuclei of the initial radius 30 µm and 40 µm which contribute dominantly to the radiated acoustic pressure. The corresponding acoustic pressure is computed using Equation (13) for the nuclei. Given that a single bubble is assumed to be a point monopole sound source, acoustic pressure due to the acceleration in bubble volume can be predicted using the time-history of bubble radius, as expressed in Equation (13). The observer position is set at the measured point, as shown in Figure 9. The next generation of nuclei follows the same trajectory as the current generation of nuclei because nuclei are repeatedly generated upstream of the wing in the steady-state flow field. Thus, the subsequent variations in the volume and location of the bubbles are repeated generation by generation. The power spectral density (PSD) of acoustic pressure is calculated using the signal processing parameters shown in

Table 4. Signal processing setup for numerical and experiment results.

Signal Processing Setup
	Numerical	Experiment
Time resolution (Δt, s)	0.5 × 10−6	3.90625 × 10−6
Frequency resolution (Δf, Hz)	11.1	1
Acquisition time (T, s)	0.1	60

.

Figure 23 compares the predicted PSD and the one-third octave band spectrum of sound pressure level (SPL) with the measured. The experimental setup was presented in Section 3. Given that the cavitation noise is impulsive, it contributes mainly to broadband noise in a high-frequency range. There is a good agreement between the two results in the frequency range exceeding 5 kHz. The difference between the two results in the frequency range higher than 5 kHz is approximately within 10 dB in terms of one-third octave band levels.

Figure 24 shows the time-variation in radius and radial velocity of a particular nucleus of initial radius, 40 µm, as it moves in the downstream direction. The lifetime of this bubble can be divided into three phases according to its dynamical characteristics. The first phase is a growing phase, the second phase is the shrinking and rebounding phase, and the final phase is a free vibration phase. Brennen separated these phases into two parts. One part corresponds to the stable acoustic cavitation, which consists of bubble oscillation with the production of subharmonics, and the other part corresponds to transient acoustic cavitation, which is composed of explosive cavitation growth and violent collapse [28]. Typically, the time-scale is shorter when bubbles shrink than that when they grow, and thus bubbles release louder acoustic pressure when they shrink as compared to that when they grow. As shown in Figure 24a,b, the radial velocity changes more rapidly in the decaying and rebounding phases, especially at the moment of transition from decaying to rebounding. Furthermore, it is evident in Figure 24c that the acoustic pressure generated by the transient acoustic cavitation step is higher than the acoustic pressure generated by a stable acoustic cavitation step, and the acoustic pressure in the decaying and rebounding phase is higher than the acoustic pressure in the growing phase. Generally, the most substantial radius reduction and acoustic pressure occur in the first decaying stage because the maximum energy is released during the first rebound process.

The Lagrangian computations are repeated to examine the effects of the number density of initial nuclei distribution on the cavitation inception pattern and noise. The number of nuclei in the water is an important parameter for determining the water quality in terms of cavitation inception. Figure 25 shows the distributions of the number concentration density of the initial nuclei used for the computations. The number of total nuclei is varied, but the relative distribution of the number density of nuclei according to their sizes are kept the same as that of seawater shown in Figure 5. The distributed pattern and location are set the same as those shown in Figure 19 and Figure 20, respectively.

Figure 26 shows the predicted tip-vortex cavitation shapes and the PSDs of the corresponding acoustic pressure for all the considered cases along with the measured spectrum. As the number density of the initial nuclei increases, more nuclei evolve into cavitation bubbles, and the corresponding acoustic pressure also becomes louder. However, the cavitation shape is more similar to the picture taken via camera, as shown in Figure 21a, as the number density of the nuclei increases to the number density of seawater. Simultaneously, the predicted acoustic spectrum shows a closer agreement with the measured spectrum as the number density of the nuclei approaches that of freshwater. These results are due to the constraints of the bubble dynamics model, which assumes the spherical symmetry for the bubble. It can be inferred from the camera-recorded cavitation shape that the assumed symmetry is broken when the nucleus evolves into a bubble. Bubbles expand into the direction of tip-cavitation and merge with adjacent bubbles. Hence, the predicted cavitation pattern differs from the experimentally visualized pattern when fewer initial nuclei are used. Some studies such as Pennings et al. [53] tackled the problem of wing-tip vortex cavitation noise by assuming stable, attached cavity and then a monopolar line source. In terms of the noise generation mechanism, fundamental physics is not different between the two approaches. The present study uses a point monopole source, while the prescribed study uses a line monopole source. Point monopole sources gathering on a line is equivalent to a line monopole source since only volume acceleration of the cavity determines the magnitude of sound source irrespective of its shape.

Figure 27 illustrates the situation where the two kinds of sources generate the same sound, although one is intermittent bubble activation, and the other is stable and attached cavity taken from Pennins et al. [53]. Since a point monopole source is an elementary source, a specific combination of point monopole sources produces any vortex cavity deformation such as a monopole, dipole, and quadrupole sources, as schematically described in Figure 27. In this regard, there is no fundamental difference in predicting the cavitation noise by using the point source and the line source. It is noted that point sources can form a line source, but the reverse is impossible. However, the goal of the present study is to predict the CIS of underwater propellers. At the beginning of the tip vortex, intermittent bubble activation occurs. The result shown in Figure 26 reveals the transition from intermittent cavitation to attached stable cavitation, as the number density of nuclei increases.

To highlight the ability of the present numerical method for predicting the CIS of the wing, additional computations are carried out by varying the cavitation number. The variation of the SPL spectral density spectrum predicted according to the cavitation numbers are shown and compared with the measured ones in Figure 28. There are good agreements between the numerical and experimental results in terms of variation of the spectrum due to the cavitation number. To find the cavitation inception number (or speed), the spectral densities shown in Figure 28 are redrawn as the octave band spectrum and regrouped into the experimental and numerical results, respectively, in Figure 29. From these results, it is seen that the broadband noise increases significantly when the cavitation number reaches $\sigma =1.4$. Since it is well known that the high-frequency broadband noise suddenly increases when the tip-vortex cavitation occurs, the cavitation inception number is estimated 1.4, which is the same for the numerical and experimental results.

5. Conclusions

To overcome the limitations of the Eulerian approach based on the homogeneous mixture model, numerical methods consisting of sequential one-way coupled applications of Eulerian and Lagrangian approaches were proposed and applied for investigating the wing-tip vortex noise of a NACA16-020 wing. A single-phase flow field around the wing was predicted using the Eulerian method consisting of the Reynolds-averaged Navier–Stokes (RANS) solver and vortex model. The cavitating tip-vortex and its noise were predicted using the Lagrangian approach based on the bubble dynamics model combined with the acoustic analogy. First, the single-phase flow field around the wing was predicted using the RANS solver. The effects of grid resolutions and turbulence models on the numerical results are analyzed in terms of location and strength of the tip-vortex core along the mean flow direction. It was observed that the predicted trajectory of the tip-vortex core center is not sensitive to the grid resolutions and the choice of turbulence models except for the standard k-ε model. However, the predicted strength of tip-vortex was sensitive to the grid-resolution and turbulence model used. The numerical solutions used more than 27 grids over the vortex core and showed convergent results in terms of minimum pressure and tangential velocity distribution across the vortex core. The numerical solution using the Reynolds stress model showed the sharpest distributions of pressure and tangential velocity across the vortex core. Furthermore, the strength of tip-vortex predicted using the Reynolds stress model maintained its value during its convection in the downstream direction. Secondly, the vortex model was applied to deal with the decay of the vortex structure along the downstream direction. This is inevitable when using the RANS solver with a moderate numerical cost. The λ₂ criteria were used to trace the vortex core line from the RANS results. It was shown that the trajectory of the vortex core center detected using the λ₂ criteria was consistent irrespective of grid resolutions and the turbulence models. Subsequently, the tip-vortex flow was regenerated by applying the Scully vortex model to the predicted background flow. The physical decay of tip-vortex, due to viscosity along with its convection in the downstream direction, was computed using the theoretical vortex dissipation model. Finally, the tip-vortex cavitation was predicted by applying the bubble dynamics equations to the regenerated tip-vortex flow field. Initial nuclei were distributed upstream of the wing. Its number density versus nuclei size was determined to match the experimentally measured property of the freshwater. The wing-tip vortex cavitation was simulated by computing the bubble radius variation and bubble position obtained from the bubble dynamics equations. Although the predicted tip-vortex cavitation shape does not exactly match the line-shape observed in the experiment, the bubbles growing from the nuclei appeared in the tip-vortex line and showed reasonable agreement with the camera-recorded image. The cavitation noise was predicted by modelling each bubble as a point monopole source and its strength was determined by its volume acceleration. The predicted PSD of sound pressure was compared with the measured PSD of sound pressure. There was good agreement between the two results. By comparing the time-history of radius and radial velocity of a single bubble with the corresponding time-variation of sound pressure, it was confirmed that the strong acoustic pressure is generated when a bubble is in the first decay phase.

To further investigate the reason for the difference between the predicted and measured cavitation shapes, the Lagrangian computations were repeated by varying the number density of the initial nuclei from freshwater to seawater conditions. The corresponding tip-vortex cavitation shape and radiated acoustic pressure were compared. It was observed that as the number density of nuclei increases from freshwater to seawater, more bubbles develop from the nuclei. Thus, the corresponding acoustic pressure also increases. As more bubbles developed, the cavitation pattern resembled a line-shape observed in the experiment. However, the predicted acoustic spectrum became higher than the measured acoustic spectrum. Given that the strength of point monopole sources assumed in the acoustic computation was proportional to its volumetric acceleration and the total volume of cavitation predicted using the lower nuclei number density was similar to that of the measured tip-vortex cavitation, it was inferred that the simulation using the freshwater condition is best matched to the actual experimental condition. Therefore, the difference between the numerical and experimental results for the tip-vortex cavitation shape was due to the limit of the spherically symmetric bubble dynamics model. Based on the camera-recorded image of the tip-vortex cavitation, it was inferred that the symmetry of the bubbles is broken during the evolution of nuclei into bubbles. Although the current numerical method exhibited difficulty in reproducing the exact shape of tip-vortex cavitation, the predicted cavitation noise matched well with the measured cavitation noise.

The additional computations were carried out for the same wing by varying the cavitation number to highlight the ability of the present numerical method to predict the cavitation inception of the wing. The predicted cavitation noise spectrums were compared with the experimental results. There were good agreements between the numerical and experimental results in terms of variation of the spectrum due to the cavitation numbers. The cavitation inception number estimated from the predicted spectrum was the same as that determined from the measured spectrum. Based on these results, the present numerical procedure, utilizing a combination of Eulerian and Lagrangian approaches, can be used as an important tool for analyzing tip-vortex cavitation inception and noise.

In future studies, the current numerical methods will be applied to investigate the scale-effects on the CIS of the wings. Furthermore, they will be used to investigate tip-vortex cavitation inception and the noise of underwater propellers.