Numerical Investigation of Tip Vortex Cavitation Inception and Noise of Underwater Propellers of Submarine Using Sequential Eulerian–Lagrangian Approaches

: In this study, the high-ﬁdelity numerical methods are developed to investigate the tip vortex cavitation (TVC) inception and noise of underwater propellers, namely, Model-A and Model-B, which are designed to investigate the e ﬀ ects of sweep angle on cavitation inception and noise. In addition, the entire body of the DARPA Subo ﬀ submarine is included to consider the e ﬀ ects of the inﬂow distortion originating from the boundary layer ﬂow of the submarine body on the cavitating ﬂow of the propellers. The Eulerian approach consisting of Reynolds-averaged Navier–Stokes (RANS) solver and the vortex model is coupled with the Lagrangian approach using the bubble dynamics equations and the acoustic analogy for nuclei initially distributed in inlet ﬂow. First, three-dimensional incompressible unsteady RANS simulations are performed to predict the hydrodynamic ﬂow ﬁeld driven by underwater propellers installed on a DARPA Subo ﬀ submarine body. The Scully vortex model and dissipation vortex model (DVM) are used to regenerate the tip vortex dissipated by artiﬁcial numerical damping and low grid resolution around the vortex core center, which is identiﬁed by using minimum λ 2 -criterion in the swirling ﬂow ﬁeld originating from the propeller blade tip. Then, tip vortex cavitation inception is simulated by applying the bubble dynamics equations to nuclei initially distributed in the inﬂow region. The volume and location of each nucleus are obtained by solving the bubble dynamics equations on the ﬂow ﬁeld obtained using the Eulerian method. Finally, the cavitation noise is predicted by modeling each bubble with a point monopole source whose strength is proportional to its volume acceleration. The validity of the present numerical methods is conﬁrmed by comparing the predicted acoustic pressure spectrum with the measured ones.


Introduction
Recently, there has been an increase in demand for propulsion systems with higher hydrodynamic performance and lower underwater-radiated noise, as environmental issues are gaining more attention in addition to the traditional military necessity. The underwater-radiated noise from underwater propellers can be grossly categorized into two parts: one is non-cavitation noise, and the other is cavitation noise. It is well known that, when the cavitation occurs on the propeller, the tonal and broadband noise increase rapidly in the wide frequency range. The sheet cavitation is of periodic characteristics at the rotating frequency, and the tip vortex cavitation (TVC) is of statistically ergodic characteristics in a broadband frequency range [1]. Therefore, the reliable and efficient numerical methods which can be used to predict and analyze the cavitation inception and noise are essential for developing low-noise and high-performance underwater propellers, in practice by increasing the so-called cavitation inception speed (CIS).
There have been a lot of previous experimental studies to predict the CIS of underwater propellers. Lim et al. [2] experimentally showed that cavitation might be caused by protrusions such as bolts in a full-scale ship and increased the cavitation inception speed by streamlining the propeller surface, including bolts. Park et al. [3] tried to delay the TVC inception using the additional structure of twisted thread on the propeller blade tip. Ahn et al. [4] evaluated the propulsion efficiency and cavitation performance of propellers of a container ship and improved the cavitation performance by decreasing the expanded area ratio and increasing the average pitch. The thru-holes near the leading edge of model propellers were recently applied to delay the cavitation inception by flowing the fluid into holes in a uniform flow condition [5]. Pereira et al. [6] tried to investigate the quantitative correlations between cavitation patterns and pressure field variation near a propeller. The cavitation patterns recorded using a high-speed camera in a cavitation tunnel were compared with the measured pressure in the flow field near the propeller when the non-uniform wake was incident on the propeller. The empirical formula using volume variation of cavitation of a propeller was suggested to predict pressure fluctuation at a location apart from the propeller. Han et al. [7] and Kim et al. [8] evaluated the cavitation inception speed using the detection of envelope modulation on noise (DEMON) analysis. The acceleration signal on the hull near the propeller was used because it is difficult to directly measure the acoustic signal. DEMON analysis, however, can identify specific components only in a particular frequency band. To resolve this limitation and quantitatively determine CIS, Lee et al. [9] extracted the cyclo-stationarity from the acceleration signal on the hull. However, the experimental method is challenging to apply in a design stage, and the scale effect for extending the experimental results obtained for the scalded models to the full-scale product is not straightforward.
The numerical methods widely used to predict the cavitation inception of wings and propellers are the Eulerian approaches consisting of the Reynolds-averaged Navier-Stokes (RANS) solvers and the large-eddy simulation (LES) techniques. The potential theory was frequently used in a design stage for underwater-propellers due to its cost-effectiveness. Seol et al. [10,11] predicted the tonal noise caused by sheet cavitation using cavitation volume variation obtained from the potential-based flow simulation. These results were validated through the comparison with the measured data. Park et al. [12] and Byeon et al. [13] verified the effectiveness of numerical methods by comparing the RANS results obtained from the commercial code, STAR-CCM+, with those obtained from the experiments in the large cavitation tunnel (LCT) of Korea Research Institute of Ships and Ocean Engineering (KRISO). The effects of the wall boundary layer profile on the contraction part of the LCT, located before the test section, were analyzed by applying the same flow distribution as the measured one on the inlet flow of the test section. Then, the wake distribution, propulsion performance, and cavitation phenomena of several propellers were estimated to focus on the hull effects. In addition, the influence of turbulence models on the TVC of a NACA66 2 -415 elliptical wing was also investigated. It was found that the use of the realizable k-ε model resulted in the overestimation of the turbulence viscosity inside the TVC by over-predicting the turbulence production. However, the Reynolds stress model (RSM) well described the vortex structure and variation in flow momentum by accurately predicting the laminar region inside the vortex core [14]. Asnaghi et al. [15] performed the LES simulation using open-source code, OpenFOAM, to simulate the TVC of the elliptical wing of NACA66 2 -415 cross-section and compared the simulation results with the measure ones. The Rayleigh-Plesset equation and the minimum pressure criterion based on the single bubble assumption were used to predict the cavitation inception, which revealed the importance of the water quality effect affecting the nuclei distribution. Cheng et al. [16] simulated the tip leakage cavitation of a NACA0009 wing using the LES and analyzed the interaction between tip leakage vortex and tip-separated vortex structures in terms of Lagrangian coherent structures. It was shown that the combination of two vortex structures generated the induced vortex, and the stretching effect of vorticity was dominant in the tip leakage cavitation flow. Bhatt et al. [17] simulated the cavitation flow on the sharp wedge using the compressible LES and homogeneous mixture assumption. The periodic pattern from cavitation generation to collapse was described. Re-entrant jet and bubbly shock wave during the transition from the sheet cavitation to the cloud cavitation were observed. Notably, it was confirmed that the supersonic flow was instantaneously formed when bubbles were collapsed, and the bubble shock wave was triggered by a compression wave generated through the disappearance of the detached cloud cavitation. Tao et al. [18] investigated the designs to delay cavitation inception in the impeller using the three-dimensional incompressible RANS. The impeller wing was simplified to a NACA0006, of which the cavitation inception was compared by varying the design parameters such as the long-short axis ratio of elliptical leading edge shape, the thickness diffusion angle, and the exponent of thickness diffusion angle. It was found that, as the leading edge shape is closer to the thick ellipse, not the arc, cavitation inception was delayed at the cost of loss of less than 5% in hydrodynamic performance.
Because LES requires a lot of computational resources and time, it is difficult to carry out the simulation effectively at the design stage. This is why RANS simulations are widely used. However, it is challenging to accurately reproduce the vortex structure convected in the downstream direction without enormous numerical cost by using the RANS solvers. It suffers from being numerically attenuated during its convection, which is caused by two reasons: one is the numerical damping intrinsic to the difference equations obtained in the discretization process of the original governing partial differential equations, and the other is excessive turbulence intensity caused by the use of the turbulence model. This problem can be mitigated by using high order schemes, enough grids, and a high-accuracy turbulence model, but it requires as much computational cost as the LES [15].
The hybrid methods where the Eulerian methods are coupled with the Lagrangian approach have been developed to simulate the cavitation flow as an alternative. These methods estimate cavitating flow by tracking the variation of nuclei distributed in the inflow region of the associated hydrodynamic flow field. Compared to the conventional Eulerian approaches, the Lagrangian approach has the advantage that it can give a simple and straightforward but reasonably comprehensive account of the cavitation inception from the bubble dynamics at a much reduced numerical cost. In the conventional computational fluid dynamics (CFD) methods such as the RANS solver and LES techniques, the cavitation is described using vapor volume fraction in the fluid under the assumption of a homogeneous mixture. However, there is intrinsic ambiguity in determining the cavitation inception by using the vapor volume fraction. If a low volume fraction value, such as vapor volume fraction α v = 0.1, is used to identify the cavitation flow at the interface between the cavitation and the surrounding fluid, the cavitation may be excessively over-predicted. If a high volume fraction value, such as α v = 0.9, is applied, the cavitation may be underestimated. There is no generally accepted criterion about which value of volume fraction is physically proper for determining the cavitation inception, especially at the beginning of cavitation. Besides, because the cavitation generation criterion depends solely on the vapor pressure, which is a function of temperature, water quality effects, such as the number density and the detailed distribution of nuclei, are not considered. To resolve this issue, some studies used modified cavitation models considering water quality effects by defining the number of nuclei per unit volume [15]. However, it is assumed that all nuclei have the same size. This means that it does not reflect the distribution of nuclei with various sizes. Regarding the nuclei distribution, O' Hern et al. [19] characterized the quality of seawater in the specific region of the Pacific Ocean in terms of the number concentration density distribution as a function of radius. Kamiirisa [20] compared the cavitation trend in terms of cavitation inception due to the nuclei distribution by applying the different air content rate in the freshwater in the cavitation tunnel. It was confirmed from this result that nuclei distribution in the freshwater could be approximately 1/10 times that of seawater. Because the Eulerian-Lagrangian coupling method can consider nuclei distribution directly in the flow field, cavitation inception can be defined quantitatively by using bubble properties, such as the bubble radius, the number of cavitation bubbles, and acoustic signal by the bubble. Additionally, water quality effects can be directly reflected by varying the initial distribution of the nuclei of different sizes.
Park et al. [21] investigated the TVC behavior of the wing having the NACA0015 cross-section using the Eulerian-Lagrangian approaches. The vortex model was used to reconstruct the tip vortex. The dissipation vortex model (DVM) was introduced to consider the vortex dissipation effect during its convection. The cavitation phenomena were simulated using the modified Rayleigh-Plesset equation for homogeneously distributed nuclei. Cavitation inception was analyzed using the cavitation bubble radius and the number of peak values in the time-history of acoustic pressure. Hsiao et al. [22][23][24] and Ma et al. [25] have developed the 3DYNAFS solver, which is an Eulerian-Lagrangian two-way coupling solver that can simulate both micro and macro variations in cavitation flow. This method can simulate bubble deformations, such as growth, elongation, separation, and collapse, during the transition process between cavitation patterns. Based on this two-way coupling method, cavitation bubble elongation along the vortex core trajectory and interaction between cavitation bubbles and re-entrant jet on the suction side of NACA0015 were simulated [26]. However, most of these methods were applied to problems involving relatively simple geometries.
In this study, the cavitation inception and noise of underwater propellers are numerically and experimentally investigated with a consideration of the full body of a submarine of a DARPF Suboff. The effects of inflow due to the submarine body on the cavitation of the underwater propeller are fully included in all of the presented numerical and experimental results so that the results in the current study are of great practical importance to the actual developers. Reliable and efficient numerical methodology is developed to predict cavitation inception and cavitation noise of submarine propellers. The proposed method is based on the one-way coupled Eulerian-Lagrangian method, of which the validity was confirmed through its application to the prediction of the tip vortex cavitation and cavitation noise of a NACA16-020 wing in the previous study [27]. Since the cavitation pattern that first occurs on the propeller is generally the TVC, the current numerical methods are focused on the prediction and analysis of the TVC and cavitation noise of underwater propellers. Two propellers of the same specifications except for the sweep angle are designed and manufactured to investigate the sweep angle effects on cavitation inception and noise. The sweep angle of one propeller is 17°, and that of the other is 34°, which are named Model-A and Model-B, respectively.
The main contributions of the present study are four-fold. First, the systematic numerical methodology is presented, which can provide reliable predictions of the tip vortex cavitation inception and noise of underwater propellers with full consideration of the entire hydrodynamic flow field of a submarine body. Second, the effects of sweep angle of the underwater propellers on cavitation inception and noise are analyzed quantitatively. It is found in both numerical and experimental results that the propeller with a higher sweep angle induces the tip vortex of less strength, generates weaker cavitation bubbles, and produces less cavitation noise. Third, the effects of incident wake from the boundary layer flow of the submarine body on the cavitating flow of underwater propellers are quantitively analyzed, with a focus on the effects of the wake flow of the rudders and the sail. The wake flow is found to reduce the incident axial flow velocity, increase the local angle of attack, and induce the cavitation on the propeller surface. Fourth, the effects of initial nuclei distributions on the cavitation flow and noise are quantitatively examined by varying the nuclei's number densities according to their size. The result reveals that the water quality described by the nuclei distribution is critical to the reliable assessment of cavitation flow and noise of underwater propellers.
The present paper is organized in the following way: In Section 2, the overall numerical procedures are presented. The Eulerian approach consisting of the RANS solver and the vortex model are described. Then, the Lagrangian approach based on the bubble dynamics model and the acoustic analogy is explained. In Section 3, the experimental method is described, which is used to capture the propeller's cavitating tip vortex and to measure hydro-acoustic pressure caused by cavitation. The numerical results are presented with the experimental results in Section 4. The effects of the inflow distortion due to submarine body, especially, the rudders and the sail on the tip vortex flow of underwater propellers are investigated in detail. Then, the tip vortex flow reconstructed by using the vortex model is presented. Finally, the tip vortex cavitation predicted by using the bubble dynamics equations and the acoustic pressure spectrum predicted by applying the acoustic analogy are presented. The predicted hydro-acoustic pressure is compared with the measured spectrum. The effects of the sweep angle of the propeller blade and initial nuclei distributions are also examined in detail.

Numerical Methods
The one-way coupled Eulerian-Lagrangian method is used in the present study. Firstly, the Eulerian method is used to predict the hydrodynamic flow field of the full submarine body with the propeller. The background flow field around the submarine body equipped with the propellers is calculated by using the RANS solver. However, the tip vortex structure predicted using the RANS solver is subject to dissipation due to the numerical damping during its convection in a downstream direction. To resolve this issue, the tip vortex flow field is reconstructed by using the Scully vortex model [27]. However, the physical dissipation of the tip vortex during its convection in a downstream direction is considered using the DVM. After the same number of tip vortex flow fields as the blade number is regenerated, the Lagrangian method is applied, which consists of the bubble dynamics model and the acoustic analogy to predict the TVC inception and cavitation noise, respectively. The initial nuclei are distributed homogeneously over a broad inlet region upstream of the propeller blade tip. The subsequent variation in the location and volume of each nucleus are traced by solving the bubble dynamics equations on the flow field obtained from the Eulerian method. For the comparison with the experiment carried out in the LCT of the KRISO, the initial size distribution of nuclei is determined by using the data for freshwater available in the literature. The time-varying data of each bubble's radius and position are fed to the acoustic solver as the input data. Assuming that each cavitation bubble is a monopole source, cavitation noise emitted by the bubble is predicted. Numerical methods used in each step are described in detail in the following subsections.

Reynolds Averaged Navier-Stokes Solver
The test propellers named Model-A and Model-B are designed to have the same specifications except for the sweep angle, which is often employed as a design parameter to delay CIS. Because the Model-A propeller has a lower sweep angle than the Model-B propeller, the latter is expected to have less cavitation and thus less cavitation noise than the former. One of the benchmarking points is whether the current numerical method can adequately describe the effects of the sweep angle on the cavitation inception and noise. Moreover, the propeller located at the stern of the submarine is exposed to wake-like disturbance produced by the boundary layer flow of the hull body. In this respect, the cavitation inception and noise of the propeller may be affected by flow disturbance caused by the submarine body. This is the second benchmarking point. To account for this, the entire body of the submarine is considered along with the propellers. Because the inflow disturbance to the propellers is caused by the upstream flow from the submarine body and appendages, the DARPA Suboff 5470 is used to model a submarine body. Geometry information for submarine body and propellers are presented in Figure 1. The dimensions of the submarine body and the propeller are listed in Table 1.
The entire computation domain and boundary conditions for the computation of the hydrodynamic flow of the submarine body with the propeller are shown in Figure 2. The dimensions of the computation domain are set the same as those of the test section of LCT in KRISO. The applied boundary conditions are determined to reproduce the experimental conditions, which is described in detail in Section 3.
caused by the submarine body. This is the second benchmarking point. To account for this, the entire body of the submarine is considered along with the propellers. Because the inflow disturbance to the propellers is caused by the upstream flow from the submarine body and appendages, the DARPA Suboff 5470 is used to model a submarine body. Geometry information for submarine body and propellers are presented in Figure 1. The dimensions of the submarine body and the propeller are listed in Table 1.

Modeled Propellers
Propeller Diameter  id density, ρ, and kinematic viscosity, ν, are equal to 997.05 kg/m 3 and 1.004 × 10 ely. The test-section sizes (L × W × H) of LCT are equal to 12.5 × 2.8 × 1.8 m. Inflow v t the inlet, gauge pressure is 0 Pa at the outlet, and operating pressure is set to be 51 e experimental background pressure. The other four side planes are set to be a sol e no-slip boundary condition is applied. The sliding mesh method is used to reprod of the propeller. The rotating speed of the propellers is set to be 25.05 rps. The cav Liquid density, ρ, and kinematic viscosity, ν, are equal to 997.05 kg/m 3 and 1.004 × 10 −6 m 2 /s, respectively. The test-section sizes (L × W × H) of LCT are equal to 12.5 × 2.8 × 1.8 m. Inflow velocity is 5 m/s at the inlet, gauge pressure is 0 Pa at the outlet, and operating pressure is set to be 51 kPa to realize the experimental background pressure. The other four side planes are set to be a solid wall where the no-slip boundary condition is applied. The sliding mesh method is used to reproduce the rotation of the propeller. The rotating speed of the propellers is set to be 25.05 rps. The cavitation number for the rotating propeller is set to be 1.78, which is given in the form, where p ∞ and p v denote the freestream pressure and vapor pressure, respectively. ρ, n and D denote liquid density, round per second, and propeller diameter. The computational domain is divided into two parts: one is the rotating part for the flow driven by the propeller, and the other is the flow field elsewhere around the submarine body. The entire grid distribution with the zoomed grids around the propellers is presented in Figure 3. The value of y + is set to 1 to accurately simulate the blockage effect and wake flow caused by the wall boundary layer on the submarine body and cavitation tunnel as well as the tip vortex flow of the propeller.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 27 The numbers of grids in the rotating part for the propeller and the non-rotating part for Suboff body are 5.8 and 4.6M, respectively. The conservation equations of mass and momentum for unsteady incompressible RANS can be written as Here, ρ, p, ⃗ and ⃗ are liquid density, pressure, velocity vector and gravitational acceleration, respectively. ⊗ denotes the tensor product and ̿ means the viscous stress tensor. The commercial code STAR CCM+ is used to numerically solve Equations (2) and (3) [12,14]. The numerical discretization is based on a finite volume-based solver. The second-order accuracy schemes are used for time and spatial derivatives. The SIMPLE method for coupling the pressure and velocity is used to satisfy the mass conservation. Reynolds stress model (RSM) is employed as the turbulent model because it accurately simulates the boundary layer and vortex structure [12]. The Schnerr-Sauer cavitation model [28] is used to consider the mass transfer between liquid and vapor. The volume of fluid (VOF) method [29] is applied to identify the cavitating flow region. The number concentration density of nuclei is considered using the Rayleigh-Plesset equation [30]. However, the detailed distribution of nuclei according to various nuclei size cannot be accounted for by the Eulerian model.

Vortex Core Identification and Vortex Model
Although partial differential equations are non-dispersive, the difference equations obtained by discretizing the partial differential equations are dispersive and dissipative. This fact implies that the numerical scheme of low order requires more numbers of grids to resolve the complex flow structure more accurately [15]. It is challenging to apply high-order numerical schemes or lots of grids in practical engineering fields because of enormous numerical cost and time, though these methods have been used to investigate relatively simple problems for academic purposes. Due to this difficulty, a lot of previous studies using the CFD techniques such as the RANS solvers have used the grids with local high resolution in the vicinity of the target object or wall to accurately simulate the boundary layer. However, the TVC of underwater propellers typically forms along a long spiral vortex trajectory in the downstream direction away from the propeller surface. This implies that the enormous numerical cost and time are required to simulate the TVC accurately using the RANS solver only. In this study, to deal with this issue, the tip vortex is reconstructed by applying the vortex The numbers of grids in the rotating part for the propeller and the non-rotating part for Suboff body are 5.8 and 4.6M, respectively. The conservation equations of mass and momentum for unsteady incompressible RANS can be written as Here, ρ, p, → u and → g are liquid density, pressure, velocity vector and gravitational acceleration, respectively. ⊗ denotes the tensor product and = τ means the viscous stress tensor. The commercial code STAR CCM+ is used to numerically solve Equations (2) and (3) [12,14]. The numerical discretization is based on a finite volume-based solver. The second-order accuracy schemes are used for time and spatial derivatives. The SIMPLE method for coupling the pressure and velocity is used to satisfy the mass conservation. Reynolds stress model (RSM) is employed as the turbulent model because it accurately simulates the boundary layer and vortex structure [12]. The Schnerr-Sauer cavitation model [28] is used to consider the mass transfer between liquid and vapor. The volume of fluid (VOF) method [29] is applied to identify the cavitating flow region. The number concentration density of nuclei is considered using the Rayleigh-Plesset equation [30]. However, the detailed distribution of nuclei according to various nuclei size cannot be accounted for by the Eulerian model.

Vortex Core Identification and Vortex Model
Although partial differential equations are non-dispersive, the difference equations obtained by discretizing the partial differential equations are dispersive and dissipative. This fact implies that the numerical scheme of low order requires more numbers of grids to resolve the complex flow structure more accurately [15]. It is challenging to apply high-order numerical schemes or lots of grids in practical engineering fields because of enormous numerical cost and time, though these methods have been used to investigate relatively simple problems for academic purposes. Due to this difficulty, a lot of previous studies using the CFD techniques such as the RANS solvers have used the grids with local high resolution in the vicinity of the target object or wall to accurately simulate the boundary layer.
However, the TVC of underwater propellers typically forms along a long spiral vortex trajectory in the downstream direction away from the propeller surface. This implies that the enormous numerical cost and time are required to simulate the TVC accurately using the RANS solver only. In this study, to deal with this issue, the tip vortex is reconstructed by applying the vortex model to the flow field obtained from the RANS solver. The unknown parameters of the vortex model, the minimum vortex core radius and maximum tangential velocity (or azimuthal velocity), are determined from the pressure and velocity fields of the RANS results. Because the vortex structure of the vortex model is assumed to be axisymmetric with an origin on the vortex core center, the vortex core identification needs to be made before the vortex reconstruction. Several criteria are available to identify the vortex core, such as minimum pressure, maximum vorticity, maximum Q-criterion, and minimum λ 2 -criterion. In particular, the λ 2 -criterion has been successfully used for identifying the coherent vortex structure for various types of problems [31]. λ 2 is defined as the second largest eigenvalue of S 2 + Ω 2 , where S and Ω are the symmetric and antisymmetric parts, respectively, of the flow velocity gradient, J ≡ ∇u. The symmetric part presents a strain rate tensor of velocity and describes the deformation rate, while the anti-symmetric part is directly related to the rate of rotation of a fluid, but not its deformation at all.
Because S 2 + Ω 2 has a symmetric matrix, it only has real eigenvalues. The symmetric and antisymmetric parts are defined as follows: Here, u, v, and w denote the velocity components in the x-, yand z-directions, respectively. The vortex core region is defined as the region where λ 2 is negative (λ 2 < 0), and the vortex core center is defined as the location where λ 2 is a minimum. The tangential velocity and pressure fields of vortex structure regenerated by using the Scully vortex model are written in the forms, respectively, [27] By separating the tangential velocity (V θ ) to x-, yand zvelocity components, the velocity field in the Cartesian coordinate can be obtained. Pressure and velocity fields in the regenerated vortex are a function of distance (r) from the vortex center with parameters of initial vortex core radius (a c ) and circulation (Γ). The vortex flow field can be divided into two areas according to its flow characteristics: one is the inner region of the vortex core, where a forced rotating flow forms, and the other is the outer region of the vortex core where a free vortex flow occurs. Some vortex models are available to define the vortex structure, such as the Rankine vortex, Scully vortex, Vatistas vortex, and Lamb-Oseen vortex. The Rankine vortex is a simple mathematical model of a vortex in a viscous fluid. However, this model has a singularity presented by a sharp velocity profile at the interface between two regions. This singularity can be overcome by using the algebraic velocity profile as Equations (7) and (8) of the Scully vortex model. The Vatistas vortex is also of the algebraic swirl-velocity profile but is based on the stationary vortex assumption. The Lamb-Oseen vortex model is an analytic solution of the one-dimensional Navier-Stokes equation, but with the assumption that the flow is laminar and all velocity components except the azimuthal velocity are zero [32]. Therefore, in this study, the Scully vortex model is used to define the vortex structure due to its merits over the others. However, since the vortex structure is physically dissipated by fluid viscosity during its convection along the downstream, the physical dissipation effects need to be considered in the reconstructed vortex field. Moore and Saffman [33] predicted the variation of tip vortex core radius along with the downstream by considering the vortex dissipation due to viscous effect, and Park et al. [21] applied this method in the vortex modeling process, named the dissipation vortex model (DVM). The formula of DVM is written as Here, a * c denotes a non-dimensional vortex core radius, and a c and ∆a c are the vortex core radius and its variation along with the downstream distance, x, respectively. D and Re D denote the chord length and the chord-based Reynolds number, respectively. The vortex core radius at the location downstream is then the sum of the initial vortex core radius (a c,init ) and the change in vortex core radius as Equation (10). The initial vortex core radius can be determined from the RANS simulation result.

Lagrangian Approach Based on Discrete Properties
The Lagrangian method is used to investigate the TVC inception and noise [23,24,30]. First, nuclei of seeds that can grow into the cavitation bubbles are initially spread out upstream of the vortex-regenerated flow field. The subsequent location and volume of each nucleus are tracked by solving the bubble dynamics equations in association with the background flow field. Finally, the acoustic pressure due to cavitation is predicted by modeling each nucleus with a point monopole source. The detailed process in each step of the particle-based approach is described in detail in the following subsections.

Initial Nuclei Distribution
Nuclei distribution varies according to water temperature, seasons, sea areas, and water types (seawater or freshwater). Experimentally, it is also dependent on the air content rate in the cavitation tunnel. Therefore, it is difficult to establish quantitative standards, but it is possible to grossly characterize the distribution of nuclei between seawater and freshwater from the available literature. O'Hern et al. [19] measured the number concentration density of nuclei as a function of particle radius in the two sea regions of Santa Catalina Island and southwest of Los Angeles by using holographic and Coulter counter measurement techniques. Figure 4 presents the number concentration density of nuclei obtained from the data for Santa Catalina Island seawater in the reference [19]. Because the experiment for cavitation inception of the scaled body is performed in the cavitation tunnel using freshwater, the nuclei distribution needs to be adjusted to represent freshwater quality. Kamiirisa [20] showed that the number of nuclei in freshwater is about ten times less than that in the seawater. In this study, the distribution of number concentration density in Santa Catalina Island seawater is used for the nuclei of radius between 30 to 100 µm, as shown in Figure 4b. However, the number of nuclei is adjusted to reflect the freshwater condition.
scaled body is performed in the cavitation tunnel using freshwater, the nuclei distribution needs to be adjusted to represent freshwater quality. Kamiirisa [20] showed that the number of nuclei in freshwater is about ten times less than that in the seawater. In this study, the distribution of number concentration density in Santa Catalina Island seawater is used for the nuclei of radius between 30 to 100 μm, as shown in Figure 4b. However, the number of nuclei is adjusted to reflect the freshwater condition.   In a flow field where the large-scale cavitation structure occurs due to a decrease in the cavitation number, the molecules on the bubble surface are sufficiently far away for one another; thus, the molecular force is negligible. Then, the critical pressure being a function of the only temperature equals the vapor pressure, which is why the cavitation flows, such as sheet and cloud cavitation, have been successfully predicted with the Eulerian methods based on the RANS solver combined with the homogeneous mixture model [34][35][36][37][38][39][40]. However, the Eulerian method has difficulty in accounting for the water quality effect because the critical pressure varies with the initial nuclei size. The scale of nuclei is generally in a range of micrometers; thus, the molecular force cannot be ignored. The molecular force between particles can be described in terms of potential energy that is a function of the distance between the particles: the molecular force is proportional to the gradient of potential energy. If the distance between molecules is too close, the potential energy increases in the direction of pushing each other. In contrast, if the distance is far, the potential energy increases in the direction of pulling each other. There is an equilibrium position between close and far, where the repulsion and attraction forces are in equilibrium. When the distance is too far away, the potential energy is constant; thus, the molecular force becomes negligible [30]. The critical pressure via which the initial nuclei size affects the cavitation inception is written in the following form [41]: where γ and p 0 denote the surface tension and the freestream pressure, respectively. It can be found from Equation (11) that the smaller the initial nuclei radius (R 0 ), the lower the critical pressure (p cr ). The fact that the ambient pressure around a nucleus must be lower than the critical pressure for the nucleus to grow into the cavitation bubble implies that the smaller the nuclei radius is, the more difficult it is for it to become a bubble. The critical pressure versus the radius of a nucleus is presented in Table 2.

Bubble Dynamics and Noise Models
The variations in the location and volume of a nucleus embedded in the background flow field are tracked by solving the bubble dynamics and trajectory equations. The Keller-Herring equation [42] is used as the bubble dynamics equation in the form: where R denotes the bubble radius, c is the speed of sound, U is the flow velocity, and U b is the bubble moving velocity. The superscript ' . ' and ' .. ' denote the single-and double-time derivatives, respectively, and the subscript 'b' and '∞' indicate the bubble and freestream flow properties, respectively. h, p, and ρ denote enthalpy, pressure, and density, respectively. The bubble radius is computed by solving Equation (12) in combination with the properties of the flow field surrounding the nucleus. The location of the nuclei (or bubble) in the flow field is computed by solving the bubble trajectory equation [43,44] in the following form: where µ l is the liquid dynamic viscosity; ρ l is the liquid density; ρ b and γ represent the bubble density and surface tension, which are set to be 1.0 kg/m 3 and 0.0742 N/m, respectively. The speed of sound in water is 1450 m/s. The first term in the right side of Equation (13) denotes the pressure gradient on the bubble surface, and the second term, gravitational acceleration ( → g ), represents the buoyancy effect. The third term is the drag force, which is a function of bubble radius(R), bubble drag coefficient (C D ) and the slip velocity. The fourth term is the lifting force proportional to the bubble lift coefficient (C L ) and the slip velocity. The last term is the Bjerknes force accounting for the bubble radius variation. The bubble drag coefficient is calculated from the empirical equation by Haberman and Morton [45], and the bubble lift coefficient is set to be 6.44 [44]. Equations (12) and (13) are numerically solved by using the 4th order Runge-Kutta method with the time step size of 0.5 µs.
In the classical Rayleigh-Plesset equation, the liquid properties surrounding the bubble are extrapolated from the point properties of liquid at the position, which is the bubble center unless there is a bubble. However, the equation sometimes allows for the bubble volume to grow infinitely because it neglects the effects of the bubble volume. Hsiao et al. [22] proposed to use the averaged liquid pressure on the bubble surface to avoid this problem, called the surface-averaged pressure (SAP) scheme. The same scheme is adopted in this study to identify the pressure (p) and velocity (U) of liquid surrounding the bubble.
To predict the flow noise due to cavitation, every single bubble is modeled as a point monopole noise source [46]. Consequently, the far-field acoustic pressure due to the TVC bubble is expressed in the following form: where ..
V b denotes the bubble volume acceleration, r is the distance between the bubble center and an observer location, and t and t-r/c denote the observer and source time (or retarded time), respectively. M r means the relative Mach number in the radiated direction.

Experimental Methods
Experiments on cavitating flow and noise of the underwater propellers are performed in the large cavitation tunnel (LCT) in the KRISO. The overall dimensions of the LCT and the test section are 60.0 (L) × 6.5 (W) × 19.8 m (H) and 12.5 (L) × 2.8 (W) × 1.8 m (H), respectively. More details on the LCT of the KRISO can be found in [27]. Figure 5a shows the SUBOFF submarine [47,48] with the Model-A propeller installed in the test section of LCT. Figure 5b schematically presents the experimental setup with the noise measurement system. The submarine is supported with two struts to secure structural stability. However, the effects of these struts on the inflow at the propeller plane were found to be negligible [49]. A B&K Type 8105 hydrophone is used to measure the radiated acoustic pressure, being located under the floor of the test section. The hydrophone is connected to the DAQ system with a frequency range of up to 102.4 kHz. The test conditions for the measurements are listed in Table 3. Here, is the inflow velocity in the test section, is the rotational speed of the propeller shaft, ; is the cavitation number at the center of the propeller, and / is the oxygen saturation that represents the air content in the water. Each measurement is conducted for one minute with the sampling rate of 262,144 samples per second, which ensures that the effect of instantaneous behavior is minimized. The measured time series data are converted to the power spectral density (PSD), using the Hanning window with a bandwidth of 1 sec. and 75% overlap.

Verification of Bubble Dynamics Solver
The bubble dynamics solver is validated by using the measured data, which are obtained from the experimental study by Muller et al. [50] to investigate single bubble behavior in detail. The single bubble was triggered in the undisturbed condition by laser. The subsequent bubble behaviors of expansion and collapse were carefully examined. The same problem is simulated by solving the bubble dynamics equation with the conditions given in Table 4.   Table 4. Figure 6 shows the simulation result in terms of the time variation of the bubble radius. The test conditions for the measurements are listed in Table 3. Here, V is the inflow velocity in the test section, n is the rotational speed of the propeller shaft, σ n;center is the cavitation number at the center of the propeller, and α/α s is the oxygen saturation that represents the air content in the water. Each measurement is conducted for one minute with the sampling rate of 262,144 samples per second, which ensures that the effect of instantaneous behavior is minimized. The measured time series data are converted to the power spectral density (PSD), using the Hanning window with a bandwidth of 1 s and 75% overlap.

Verification of Bubble Dynamics Solver
The bubble dynamics solver is validated by using the measured data, which are obtained from the experimental study by Muller et al. [50] to investigate single bubble behavior in detail. The single bubble was triggered in the undisturbed condition by laser. The subsequent bubble behaviors of expansion and collapse were carefully examined. The same problem is simulated by solving the bubble dynamics equation with the conditions given in Table 4. The bubble simulation is started from the initial conditions, R b (t 0 ) = R max and . R b (t 0 ) =0 at the time when the bubble radius is maximum, that is t 0 = t max = 70.7 µs. The maximum bubble radius, R max , is equal to 747 µm, and the time interval is 5 µs. The detailed computation conditions are listed in Table 4. Figure 6 shows the simulation result in terms of the time variation of the bubble radius.
ci. 2020, 10, x FOR PEER REVIEW 13 Figure 6. Comparison of numerical bubble radius with experiment. t is seen from the experimental result that a single bubble shrinks and rebound after it expa ubble behavior is well followed by the numerical result obtained by solving the Keller-He e dynamics equation. The measured minimum radius is less than 12 μm, and the calcu um radius is 8.8 μm at the end of the collapse. These results confirm the validity of the bu ics model used in the present study.

umerical RANS Solutions
igure 7 compares the RANS simulation results with the experimental results in term ting flow shapes. The snapshots of cavitating flow shape recorded by using the high-s ra in the test section of the LCT are compared with the iso-surface of vapor volume frac ed from the numerical simulation results. Although the TVC and sheet cavitation on n sides of propeller blades are well captured, the cavitating tip vortex trajectory obse y in the experiments cannot be found in the corresponding numerical results. It is challen ture the cavitating tip vortex trajectory by using the RANS solver without enormous nume This difficulty will be overcome effectively by applying the vortex model and the bu ics and trajectory equations. It is seen from the experimental result that a single bubble shrinks and rebound after it expands. This bubble behavior is well followed by the numerical result obtained by solving the Keller-Herring bubble dynamics equation. The measured minimum radius is less than 12 µm, and the calculated minimum radius is 8.8 µm at the end of the collapse. These results confirm the validity of the bubble dynamics model used in the present study. Figure 7 compares the RANS simulation results with the experimental results in terms of cavitating flow shapes. The snapshots of cavitating flow shape recorded by using the high-speed camera in the test section of the LCT are compared with the iso-surface of vapor volume fractions obtained from the numerical simulation results. Although the TVC and sheet cavitation on the suction sides of propeller blades are well captured, the cavitating tip vortex trajectory observed clearly in the experiments cannot be found in the corresponding numerical results. It is challenging to capture the cavitating tip vortex trajectory by using the RANS solver without enormous numerical cost. This difficulty will be overcome effectively by applying the vortex model and the bubble dynamics and trajectory equations. camera in the test section of the LCT are compared with the iso-surface of vapor volume fractions obtained from the numerical simulation results. Although the TVC and sheet cavitation on the suction sides of propeller blades are well captured, the cavitating tip vortex trajectory observed clearly in the experiments cannot be found in the corresponding numerical results. It is challenging to capture the cavitating tip vortex trajectory by using the RANS solver without enormous numerical cost. This difficulty will be overcome effectively by applying the vortex model and the bubble dynamics and trajectory equations.  Further analysis is carried out to confirm the validity of numerical results in terms of hydrodynamic performances. Table 5 compares the predicted thrust with the measured ones. The thrust is compared in terms of non-dimensional values: the thrust coefficient of K T (= T ρn 2 D 3 ). It can be seen that the maximum relative percentage error is less than 2%.  Figure 8 shows the distribution of axial inflow velocity in the region 0.14D upstream from the propeller tip, where D denotes the propeller diameter. The axial velocity is averaged during one rotation on the circle of the propeller tip radius, and its distribution is presented according to the azimuthal angle defined on the left side of Figure 8. The axial velocity distribution shows the deep narrow valleys at the circular positions of 0°, 90°, 180°, and 270°, which coincide with the locations of the rudders. Moreover, the deepest valley is located at 90°, which manifests the additional effects of the sail. Because the reduced axial velocity at these locations locally increases the angle of attack of the propeller blade, the onset of cavitation can occur more likely on the surface of the propeller blade passing through these locations, mostly at 90°. This phenomenon affects the so-called sheet cavitation inception, as shown in Figure 7. The sheet cavitation inception was not expected because the targeted propellers were originally designed to investigate the tip vortex cavitation that is known to occur firstly among the other types of cavitation. However, it is confirmed that the hydrodynamic flow field obtained from the RANS simulation includes the effects of wake flow originating from the boundary layers of the submarined body, especially, the rudders and the sail. Further analysis is carried out to confirm the validity of numerical results in terms of hydrodynamic performances. Table 5 compares the predicted thrust with the measured ones. The thrust is compared in terms of non-dimensional values: the thrust coefficient of KT (= ). It can be seen that the maximum relative percentage error is less than 2%.  Figure 8 shows the distribution of axial inflow velocity in the region 0.14D upstream from the propeller tip, where D denotes the propeller diameter. The axial velocity is averaged during one rotation on the circle of the propeller tip radius, and its distribution is presented according to the azimuthal angle defined on the left side of Figure 8. The axial velocity distribution shows the deep narrow valleys at the circular positions of 0˚, 90˚, 180˚, and 270˚, which coincide with the locations of the rudders. Moreover, the deepest valley is located at 90˚, which manifests the additional effects of the sail. Because the reduced axial velocity at these locations locally increases the angle of attack of the propeller blade, the onset of cavitation can occur more likely on the surface of the propeller blade passing through these locations, mostly at 90˚. This phenomenon affects the so-called sheet cavitation inception, as shown in Figure 7. The sheet cavitation inception was not expected because the targeted propellers were originally designed to investigate the tip vortex cavitation that is known to occur firstly among the other types of cavitation. However, it is confirmed that the hydrodynamic flow field obtained from the RANS simulation includes the effects of wake flow originating from the boundary layers of the submarined body, especially, the rudders and the sail.

Vortex Core Trajectory
Before the tip vortex is regenerated by using the vortex model, the vortex core center needs to be identified for the swirling flow from propeller blades. Coherent vortex structure is identified using the λ2-criterion of the RANS simulation results. The vortex core center is defined as the location of the minimum value of λ2. Figure 9 presents the distribution of the λ2-criterion on the circular plane just behind the propeller. Figures 9a shows the overall distribution of the λ2-criterion in a range

Vortex Core Trajectory
Before the tip vortex is regenerated by using the vortex model, the vortex core center needs to be identified for the swirling flow from propeller blades. Coherent vortex structure is identified using the λ 2 -criterion of the RANS simulation results. The vortex core center is defined as the location of the minimum value of λ 2 . Figure 9 presents the distribution of the λ 2 -criterion on the circular plane just behind the propeller. Figure 9a shows the overall distribution of the λ 2 -criterion in a range between −2400 and 2400. It is seen that the seven similar regions are formed in the flow field, as can be expected from the fact that the propeller consists of seven blades. Under the understanding that the large negative value of the λ 2 indicates the existence of a strong coherent vortex structure, the regions of high negative values are detailed by adjusting the contour levels, together with the projected view of the propeller blades, as shown in Figure 9b. In addition to the contours of the λ 2 -criterion, the circular band contours in the center region denoting the circular angle are added to specify the vortex core location. The vortex core centers are identified by using the minimum value of the λ 2 for each of the blades, and their exact locations are presented with the azimuthal angles, as shown in Figure 9b. Figure 10 presents the distribution of the λ 2 on the black-colored circle shown in Figure 9b. The seven sharp peaks can be identified in the curve, which implies that the vortex center can be identified with high resolution. The trajectories of the tip vortex core are formed by tracing the locations of the minimum λ 2 values identified in this way along with the downstream direction. Figure 11 shows the formed trajectory curves of the tip vortex for two types of propellers.  Figure 9b. Figure 10 presents the distribution of the λ2 on the black-colored circle shown in Figure 9b. The seven sharp peaks can be identified in the curve, which implies that the vortex center can be identified with high resolution. The trajectories of the tip vortex core are formed by tracing the locations of the minimum λ2 values identified in this way along with the downstream direction. Figure 11 shows the formed trajectory curves of the tip vortex for two types of propellers.
(a) wide contour range (b) contour range less than 0    Figure 9b. Figure 10 presents the distribution of the λ2 on the black-colored circle shown in Figure 9b. The seven sharp peaks can be identified in the curve, which implies that the vortex center can be identified with high resolution. The trajectories of the tip vortex core are formed by tracing the locations of the minimum λ2 values identified in this way along with the downstream direction. Figure 11 shows the formed trajectory curves of the tip vortex for two types of propellers.
(a) wide contour range (b) contour range less than 0    The tip vortex core trajectory is known to be of a spiral curve along the downstream directio is seen in Figure 11 that the vortex core trajectory predicted using the minimum λ2-criterion of th NS results well follows this curve. Moreover, the predicted trajectories of the tip vortex core see be similar to the experimental ones shown in Figure 7. However, it noted that the vortex co jectory predicted for the Model-A is attached to the blade tip while detached from the Model- The tip vortex core trajectory is known to be of a spiral curve along the downstream direction. It is seen in Figure 11 that the vortex core trajectory predicted using the minimum λ 2 -criterion of the RANS results well follows this curve. Moreover, the predicted trajectories of the tip vortex core seem to be similar to the experimental ones shown in Figure 7. However, it noted that the vortex core trajectory predicted for the Model-A is attached to the blade tip while detached from the Model-B blade tip. This difference is one of the points to characterize the flow driven by the propellers with different sweep angles.
For further analysis, Figure 12 presents the distribution of the minimum values of λ 2 according to the blade location for two propellers. It is seen that the minimum values of λ 2 for all blades of the Model-B propeller are lower than those of the Model-A propeller. This result reveals that the tip vortex of Model-B with a higher sweep angle is weaker than that of Model-A with a lower sweep angle, and the tip vortex cavitation is also possibly delayed for Model-B in comparison with Model-A. Figure 13 shows the zoomed plot of the starting point of the vortex core trajectory for each propeller shown in Figure 11. It is confirmed that the vortex core trajectory of the Model-B propeller is formed more apart from the blade tip. These results reveal that the underwater propeller of a higher sweep angle induces a weaker tip vortex than a lower sweep angle.
ci. 2020, 10, x FOR PEER REVIEW 16 tip. This difference is one of the points to characterize the flow driven by the propellers ent sweep angles. or further analysis, Figure 12 presents the distribution of the minimum values of λ2 accor blade location for two propellers. It is seen that the minimum values of λ2 for all blades o l-B propeller are lower than those of the Model-A propeller. This result reveals that th x of Model-B with a higher sweep angle is weaker than that of Model-A with a lower sw , and the tip vortex cavitation is also possibly delayed for Model-B in comparison with M gure 13 shows the zoomed plot of the starting point of the vortex core trajectory for ller shown in Figure 11. It is confirmed that the vortex core trajectory of the Model-B prop ed more apart from the blade tip. These results reveal that the underwater propeller of a hi angle induces a weaker tip vortex than a lower sweep angle.    Figure 14 shows the iso-contours of flow velocity components at the plane of x = −3.15 m, which is located just behind the propeller. The steepest gradient is observed in the distribution of the x-directional velocity component. Since the x-direction is parallel to the freestream velocity and almost normal to the propeller rotational plane, the vortex structure is better captured in terms of the x-directional velocity component than the others. Therefore, vortex modeling is made by using the velocity component in the x-direction on the plane normal to the mean flow direction.

Tip Vortex Reconstruction
Appl. Sci. 2020, 10, x FOR PEER REVIEW 17 of 27 Figure 14 shows the iso-contours of flow velocity components at the plane of x = −3.15m, which is located just behind the propeller. The steepest gradient is observed in the distribution of the xdirectional velocity component. Since the x-direction is parallel to the freestream velocity and almost normal to the propeller rotational plane, the vortex structure is better captured in terms of the xdirectional velocity component than the others. Therefore, vortex modeling is made by using the velocity component in the x-direction on the plane normal to the mean flow direction.

Tip Vortex Reconstruction
(a) x-direction velocity (b) y-direction velocity (c) z-direction velocity The vortex core radius and maximum tangential velocity of the vortex model need to be determined to regenerate tip vortex structure by using the vortex model. The procedure to define these values is schematically presented in Figure 15. The distribution of the axial velocity component, which corresponds to the tangential velocity of the vortex model, is extracted from the black-dashed lines between the red-solid lines shown in Figure 15a. The extracted distributions for seven cores are shown in Figure 15b. Here, the vortex core diameter is defined as the distance between the locations of maximum and minimum velocities, and then naturally, the initial vortex core radius is set to be The vortex core radius and maximum tangential velocity of the vortex model need to be determined to regenerate tip vortex structure by using the vortex model. The procedure to define these values is schematically presented in Figure 15. The distribution of the axial velocity component, which corresponds to the tangential velocity of the vortex model, is extracted from the black-dashed lines between the red-solid lines shown in Figure 15a. The extracted distributions for seven cores are shown in Figure 15b. Here, the vortex core diameter is defined as the distance between the locations of maximum and minimum velocities, and then naturally, the initial vortex core radius is set to be half of the diameter. The initial tangential velocity is also defined as the difference between the maximum and minimum velocities, as denoted by the blue text and arrow-line in Figure 15b. Because the velocity is zero at the vortex center, the mean value of the maximum and minimum velocities, as denoted by the green dashed line in Figure 15b, is considered the mean flow velocity in the x-direction. The vortex structure described by Equations (7) and (8) can be determined using the initial vortex core radius and tangential velocity determined in this way. After the vortex core trajectory and the properties of the vortex model are determined, the domains where the tip vortex is reconstructed need to be defined. Figure 16 shows the vortex domain constructed for each blade, together with its construction process. The vortex domain is constructed by using the cuboid of the square cross-section, whose center coincides with the vortex core line, as shown in Figure 16a. This cuboid is bent to follow the spiral curve of the vortex core trajectory, as shown in Figure 16c. The cross-section of the bent cuboid at the plane of the constant azimuthal angle is shown in Figure 16b. The tip vortex for each blade is reconstructed by extending the vortex model in the bent cuboid.  After the vortex core trajectory and the properties of the vortex model are determined, the domains where the tip vortex is reconstructed need to be defined. Figure 16 shows the vortex domain constructed for each blade, together with its construction process. The vortex domain is constructed by using the cuboid of the square cross-section, whose center coincides with the vortex core line, as shown in Figure 16a. This cuboid is bent to follow the spiral curve of the vortex core trajectory, as shown in Figure 16c. The cross-section of the bent cuboid at the plane of the constant azimuthal angle is shown in Figure 16b. The tip vortex for each blade is reconstructed by extending the vortex model in the bent cuboid. Figure 17 compares the reconstructed tip vortex flow with the RANS result in terms of the axial velocity distribution on the six cross-sections of the vortex domain. It is seen from Figure 17a that the axial velocity distributions on the plane nearest the blade tip are almost the same between the two results. However, it can be found in Figure 17b-f that, as the fluid flows in the downstream direction, the vortex structure in the RANS result disappears quickly, while the vortex structure reconstructed using the vortex model keeps almost its original one. Note that the marginal vortex dissipation occurs in the regenerated tip vortex due to the influence of DVM. by using the cuboid of the square cross-section, whose center coincides with the vortex core line, as shown in Figure 16a. This cuboid is bent to follow the spiral curve of the vortex core trajectory, as shown in Figure 16c. The cross-section of the bent cuboid at the plane of the constant azimuthal angle is shown in Figure 16b. The tip vortex for each blade is reconstructed by extending the vortex model in the bent cuboid.

Nuclei Population
In this study, the effects of water quality are accounted for by adjusting the number concentration density (NCD) distribution of nuclei versus its size, along with which the critical pressure varies. The NCD of nuclei is determined to reflect the freshwater quality used in the experiment in the LCT. O'Hern et al. [19] reported the measured NCD of nuclei in seawater, and Kamiirisa [20] compared the NCD of nuclei between seawater and freshwater. The NCD of nuclei for the present water is determined by using these data, which is 10% of the NCD in seawater, as shown in Figure 4. Nuclei are initially distributed homogeneously in the region just upstream of the propeller tip, as shown in Figure 18a. The region is constituted by the donut-shaped domain to ensure the entrance of the nuclei into the tip vortex flow region shown in Figure 18b. Figure 19 illustrates the detailed distribution of nuclei according to their size.

Nuclei Population
In this study, the effects of water quality are accounted for by adjusting the number concentration density (NCD) distribution of nuclei versus its size, along with which the critical pressure varies. The NCD of nuclei is determined to reflect the freshwater quality used in the experiment in the LCT. O'Hern et al. [19] reported the measured NCD of nuclei in seawater, and Kamiirisa [20] compared the NCD of nuclei between seawater and freshwater. The NCD of nuclei for the present water is determined by using these data, which is 10% of the NCD in seawater, as shown in Figure 4. Nuclei are initially distributed homogeneously in the region just upstream of the propeller tip, as shown in Figure 18a. The region is constituted by the donut-shaped domain to ensure the entrance of the nuclei into the tip vortex flow region shown in Figure 18b. Figure 19 illustrates the detailed distribution of nuclei according to their size. igure 4. Nuclei are initially distributed homogeneously in the region just upstream o eller tip, as shown in Figure 18a. The region is constituted by the donut-shaped domain to e ntrance of the nuclei into the tip vortex flow region shown in Figure 18b. Figure 19 illus etailed distribution of nuclei according to their size.

Bubble Dynamics
The tip vortex cavitation of the underwater propeller is predicted by solving the equations of bubble dynamics and trajectory for each nucleus in combination with the information of the background hydrodynamic flow field in which the regenerated tip vortex flow is embedded. Figure  20 shows the predicted cavitation inception of nuclei at the three-time instants for two propellers. The population number of nuclei for a single generation is 991 ea, three generations are used for the simulation, and thus the total number of nuclei is 2973ea. Because the cavitation bubbles are very small when compared to the propeller's dimension, the bubbles shown in Figure 20 are enlarged to twenty times the actual size. The cavitation bubbles from the Model-A propeller appear to be larger than those from the Model-B propeller. This fact is ensured in the results shown in Figure 21.

Bubble Dynamics
The tip vortex cavitation of the underwater propeller is predicted by solving the equations of bubble dynamics and trajectory for each nucleus in combination with the information of the background hydrodynamic flow field in which the regenerated tip vortex flow is embedded. Figure 20 shows the predicted cavitation inception of nuclei at the three-time instants for two propellers. The population number of nuclei for a single generation is 991 ea, three generations are used for the simulation, and thus the total number of nuclei is 2973ea. Because the cavitation bubbles are very small when compared to the propeller's dimension, the bubbles shown in Figure 20 are enlarged to twenty times the actual size. The cavitation bubbles from the Model-A propeller appear to be larger than those from the Model-B propeller. This fact is ensured in the results shown in Figure 21.
The population number of nuclei for a single generation is 991 ea, three generations are used for the simulation, and thus the total number of nuclei is 2973ea. Because the cavitation bubbles are very small when compared to the propeller's dimension, the bubbles shown in Figure 20 are enlarged to twenty times the actual size. The cavitation bubbles from the Model-A propeller appear to be larger than those from the Model-B propeller. This fact is ensured in the results shown in Figure 21.  However, in both cases, the production patterns of TVC predicted using the current numerical method are different from those observed in the experiment, which forms the same spiral lines as the tip vortex core lines, shown in Figures 7 and 11. The difference is due to the assumed spherical symmetry of the bubble dynamics model. When an actual TVC occurs, the cavitation bubbles are known to be stretched along the vortex core trajectory, and merge with the surrounding bubbles and grow, not only in a radial variation [26]. A more complex bubble dynamics model may be needed to reflect the asymmetric deformation of and merging between the bubbles. However, the main goal of the current study is to develop an efficient and reliable numerical methodology for the prediction of the tip vortex cavitation inception and noise. The cavitation inception and the magnitude of the cavitation noise source, which are monopolar characteristics, are not affected by the cavitation's However, in both cases, the production patterns of TVC predicted using the current numerical method are different from those observed in the experiment, which forms the same spiral lines as the tip vortex core lines, shown in Figures 7 and 11. The difference is due to the assumed spherical symmetry of the bubble dynamics model. When an actual TVC occurs, the cavitation bubbles are known to be stretched along the vortex core trajectory, and merge with the surrounding bubbles and grow, not only in a radial variation [26]. A more complex bubble dynamics model may be needed to reflect the asymmetric deformation of and merging between the bubbles. However, the main goal of the current study is to develop an efficient and reliable numerical methodology for the prediction of the tip vortex cavitation inception and noise. The cavitation inception and the magnitude of the cavitation noise source, which are monopolar characteristics, are not affected by the cavitation's detailed shape. This fact justifies the current approach, which will be confirmed by the following results.
The time trace of bubble radius according to its initial size is compared between the Model-A and Model-B propellers in Figure 21. The number of nuclei developing into cavitation bubbles is most for the initial nuclei size of 30 and 40 µm in both cases. Because the tip vortex of the Model-B is delayed and of lower strength, as shown in Figures 11-13, the maximum radius of the cavitation bubble induced by the Model-B propeller is smaller than that of the Model-A propeller. However, the total number of the cavitation bubble is almost the same for both cases, as shown in Table 6. The acoustic pressure due to tip vortex cavitation is computed using the information of cavitation bubbles predicted for the Model-A and Model-B propellers. Each cavitation bubble is assumed to be a point monopole source. The observer location is set to be the same as that of the corresponding experiment shown in Figure 5. The predicted acoustic pressures in the time and frequency domains are shown in Figure 22. Furthermore, the power spectral densities of acoustic pressures are compared with the measured data in the LCT of KRISO. The sampling frequency and the frequency resolution of the predicted results are 2000 kHz and 40 Hz, respectively. The Hanning window is applied to reduce leakage and aliasing. Because the noise due to TVC contributes to the spectral components in a broad frequency range over 1 kHz, the acoustic pressure is compared in the frequency range of 1 to 100 kHz. Although the predicted spectral levels are slightly less than the measured ones, there are fair agreements between the two results in terms of spectral slope over 5 kHz. Moreover, the spectral levels of Model-A are higher than those of Model-B in both numerical and experimental results.
The actual number density of nuclei is known to be dependent on the air content rate in the cavitation tunnel. However, because there is no clear relationship between air content rate in the freshwater and the number density of nuclei, additional computations are carried out to investigate the effects of initial distribution of nuclei on cavitation noise of the underwater propellers. The considered cases are for the NCD of initial nuclei: 33.3% and 100% of seawater condition. It is noted that the case presented in Figure 22 is 10% of the seawater condition. The numerical results predicted using these distributions of nuclei are presented and compared with the experimental ones in Figure 23. In all the cases, the slope of spectrum is almost the same, and the acoustic level due to the cavitation flow of the Model-B is less than that of the Model-A. In the case using 33.3% of the seawater condition, there are good agreements in the frequency ranges over 3.5 kHz between the numerical and experimental results for both propellers.
For more quantitative comparison, the overall sound pressure levels (OASPL) are computed in the given frequency range for the reference pressure of 10 −6 Pa. Table 7 compares the SPLs between the computed and measured data for Model-A and Model-B. It is found that the overall SPL predicted for the Model-B is less by 6.3 to 7.1 dB than that of Model-A, while the measured one for Model-B is less by 2.7 dB than that of Model-A. These results reveal that the current numerical method can adequately appreciate the effects of sweep angle on the cavitation inception and noise of the underwater propeller of the submarine. spectral components in a broad frequency range over 1 kHz, the acoustic pressure is compared in the frequency range of 1 to 100 kHz. Although the predicted spectral levels are slightly less than the measured ones, there are fair agreements between the two results in terms of spectral slope over 5 kHz. Moreover, the spectral levels of Model-A are higher than those of Model-B in both numerical and experimental results. The actual number density of nuclei is known to be dependent on the air content rate in the cavitation tunnel. However, because there is no clear relationship between air content rate in the freshwater and the number density of nuclei, additional computations are carried out to investigate the effects of initial distribution of nuclei on cavitation noise of the underwater propellers. The considered cases are for the NCD of initial nuclei: 33.3% and 100% of seawater condition. It is noted that the case presented in Figure 22 is 10% of the seawater condition. The numerical results predicted using these distributions of nuclei are presented and compared with the experimental ones in Figure  23. In all the cases, the slope of spectrum is almost the same, and the acoustic level due to the cavitation flow of the Model-B is less than that of the Model-A. In the case using 33.3% of the seawater condition, there are good agreements in the frequency ranges over 3.5 kHz between the numerical and experimental results for both propellers.  The actual number density of nuclei is known to be dependent on the air content rate in the cavitation tunnel. However, because there is no clear relationship between air content rate in the freshwater and the number density of nuclei, additional computations are carried out to investigate the effects of initial distribution of nuclei on cavitation noise of the underwater propellers. The considered cases are for the NCD of initial nuclei: 33.3% and 100% of seawater condition. It is noted that the case presented in Figure 22 is 10% of the seawater condition. The numerical results predicted using these distributions of nuclei are presented and compared with the experimental ones in Figure  23. In all the cases, the slope of spectrum is almost the same, and the acoustic level due to the cavitation flow of the Model-B is less than that of the Model-A. In the case using 33.3% of the seawater condition, there are good agreements in the frequency ranges over 3.5 kHz between the numerical and experimental results for both propellers. For more quantitative comparison, the overall sound pressure levels (OASPL) are computed in the given frequency range for the reference pressure of 10 −6 Pa. Table 7 compares the SPLs between the computed and measured data for Model-A and Model-B. It is found that the overall SPL predicted for the Model-B is less by 6.3 to 7.1 dB than that of Model-A, while the measured one for Model-B is less by 2.7 dB than that of Model-A. These results reveal that the current numerical method can adequately appreciate the effects of sweep angle on the cavitation inception and noise of the underwater propeller of the submarine.

Conclusions
A one-way coupled method between Eulerian and Lagrangian approaches was developed to efficiently and accurately predict the tip vortex cavitation (TVC) inception and noise of submarine propellers. Notably, the two propellers of the same specifications except for sweep angle were considered to investigate its effects on the TVC inception and noise. Moreover, the realistic incident hydrodynamic flow field to the propeller was realized by including the entire submarine body to investigate the effects of incident disturbance on the TVC inception and noise.
The Eulerian approach based on the sequential application of the RANS solver and the vortex model was employed to predict the highly resolved tip vortex flow of the submarine propellers. The RANS results were validated by comparing the predicted hydrodynamic performances with the measured ones. The computed axial velocity distribution revealed the structure of incident wake flow to the propeller, which matched the layout of the rudders and the sail location. This incident wake flow affected the cavitating flow of the propeller by increasing the local angle of attack. However, the tip vortex structure obtained from the RANS results showed inadequate agreements with the measured ones due to excessive numerical damping. The vortex model is applied to improve the predicted resolution of the tip vortex flow field. The vortex core trajectory was first identified by tracing the λ 2 -distribution of the RANS results and defining the core center as the minimum λ 2 . Then, the Scully vortex model with the DVM was applied to reproduce a well-organized vortex structure. The tip vortex flow regenerated by using the vortex model also showed good agreement with the recorded snapshot of the tip vortex. Additionally, through the comparison of the minimum λ 2 values for all blades of the Model-A and Model-B propellers, it was found that the larger sweep angle leads to tip vortex of less strength.
The Lagrangian method based on the sequential application of bubble dynamics and trajectory equations for nuclei and of the acoustic equation for the bubbles was applied to predict the TVC inception and noise. The tip vortex cavitation inception of the propellers was predicted by numerically solving the bubble dynamics and bubble trajectory equations for initial nuclei in association with the regenerated tip vortex flow fields. The freshwater quality was reflected by adjusting the number density concentration (NDC) of nuclei to the 1/10 of the NDC reported for Santa Catalina Island. Finally, the acoustic pressure due to tip vortex cavitation was predicted by assuming each bubble as a point monopole source. The predicted spectral density of the acoustic pressure was compared with the measured one. There is good agreement between the predicted and measured ones in terms of spectral slopes in the high-frequency range. Furthermore, the predictions reveal that the TVC noise of Model-B was less than that of Model-A, which also matched the measured result. It was seen that the underwater propeller with a higher sweep angle reduces the TVC noise by decreasing tip vortex strength. Further computations were carried out to investigate the effects of the number concentration density of nuclei. Irrespective of the NCD, the predicted spectral levels for Model-A were higher than those for Model-B.
It was found that there was a noticeable difference in the large structure of TVC between the predicted and measured results. This seems to be due to the restriction of bubble dynamic equations assuming the spherically symmetric shape of the bubble. To reproduce the camera-recoded TVC structure, more complicated models for bubble dynamics may be needed but at the increased numerical cost. However, since the TVC inception and noise are not significantly affected by the detailed shape of TVC, the TVC inception and noise can be successfully predicted by using the current numerical methods. Consequently, the present numerical procedure can be used as an efficient and reliable tool to investigate the TVC inception and noise of the underwater propellers in their early development. A future study aims to develop more advanced cavitation models that reproduce more realistic tip vortex cavitation behavior.
Author Contributions: C.C. provided the fundamental idea for this study. G.K. and I.P. performed the numerical simulations and analyzed the numerical results. H.S. carried out the experiments. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Defense Acquisition Program Administration and Agency for Defense Development.