1   Introduction

The construction of an increasing number of high-speed civil vessels and naval equipment has increased the intensity of underwater noise, which poses a harmful impact on the marine environment and ecosystem. Among the various noise sources, propeller noise is regarded as a dominant source because it not only contributes a significant noise level but also exhibits unique acoustic characteristics—a broadband, random signal superimposed with several discrete peaks of strong target strength—which threatens the stealth capability of submarines (Stewart et al., 2015; Yao et al., 2020; Wang et al., 2020b). Recent developments have resulted in noncavitating propeller designs, which shift the focus to the noncavitating characteristics of propellers, such as hydrodynamic noise and the hull excitation it induces. Discrete tones at the blade passing frequency and its harmonics are caused by periodic force fluctuations when operating in the wake field. Conversely, broadband noise is generated due to irregular force pulsations caused by turbulent ingestion (Nathan et al., 2017). Considering that the decay of acoustic power observed in the low-frequency range is relatively slow, the prediction and optimization of low-frequency broadband noise have consistently been a research focus. Moreover, the noise source—low-frequency broadband thrust—induces hull pressure pulses and potential resonance (Ge, 2021; Xiong et al., 2021; Xiong et al., 2022).

The earliest and most widely used approach for hydroacoustics is Lighthill's equation, which is originally rooted in aeroacoustics; this equation was later extended to the FW-H equation, which includes the effects of moving surfaces (Testa et al., 2018; Savas et al., 2021; Si et al., 2023; Wu et al., 2018). The right-hand side can be interpreted as three sound source terms: a monopole due to the moving object, a dipole due to the unsteady loading on bodies, and a quadrupole generated by fluctuating Reynolds stresses in a volume; these terms can be obtained by experimental measurements and computational fluid dynamics methods (Wang et al., 2006). As early as the 2070s, Sevik (Sevik, 1973; Paul and Uhlman, 2013; Wu and Huang, 2021) experimentally investigated the broadband-spectrum force of a blade model propeller in a homogeneous turbulent flow field with "humps" or "haystacks" at multiples of the blade rate. This work has since served as a benchmark for assessing numerical prediction methods. However, experimentally detecting low-amplitude broadband force sources or noise is extremely challenging. Moreover, the scale effects (Jukka and Antti 2022) between the model scale and the full scale are not clearly interpreted. With advancements in computational fluid mechanics, comprehensive sound sources can be acquired for input into acoustic analogies. Jiang et al. (2020) used a modified strip method to predict low-frequency discrete and broadband-spectrum forces with the assumption of inviscid and irrotational flow. Wang et al. (2020a) examined the sound field properties of noncavitating marine propellers based on the boundary element method approach, which suggests the dominance of the loading term in nonuniform inflow conditions. However, the thickness term prevails in open-water scenarios.

Another broad category involves the direct approach to simultaneously compute the hydrodynamic and hydroacoustic fields by solving the compressible flow equations through direct numerical simulation (DNS) or large-eddy simulation (LES) (Al-Am et al., 2022). The immersed boundary method (IBM), as a popular tool for imposing boundary conditions on boundaries that do not align with the computational grid, reduces the effort of grid generation and shows natural potential for coupling the hydrodynamic and hydroacoustic fields (Tran and Plourde et al., 2014). Given that the wavelength scale is much larger than that of the turbulent structure, the simulation domain must be sufficiently large to encompass all the sound sources of interest and at least part of the acoustic near field. On this basis, extension to the acoustic far field can then be achieved using various analytical and numerical means (Bodony and Lele, 2005). Initially, the direct prediction of flow noise concentrates on pure turbulent noise. Bogey et al. (2003) computed the sound field of a Mach 0.9 jet at a Reynolds number of 6.5 × 104 using LES and validated their flow and sound statistics with experimental data. Seo and Mittal (2011) developed a sharp interface IBM for the linearized perturbed compressible equations and validated the method with flow-induced noise from stationary objects and acoustic scattering problems. Accurately including the effects of nonlinear flow and reflections at the domain boundaries remains challenging when dealing with noise generation around moving bodies immersed in flow. Subsequently, Schlanderer et al. (2017) introduced a virtual boundary method for compressible viscous fluid flow, which can precisely represent noise radiation from two-dimensional stationary and oscillating cylinders without any penalty in the allowable time step. Insights into the near-field noise propagation illuminate noise-control strategies.

A highly-resolved ILES simulation for compressible flow based on the BDIM is used to perform a flow-induced noise prediction and explore the generation mechanism in the source region of a model-scale underwater rotating propeller, which avoids modeling approximations from turbulence models. The noise directivity and low-frequency broadband spectra of the primary dipole noise source—integral thrust fluctuation due to turbulent injection—are accurately captured, with an accurate representation of the effect of turbulence injection.

2   Methodology

2.1   Governing equations

To derive the BDIM meta equation, a domain that includes a solid body subdomain Ω_b and a compressible fluid subdomain Ω_f is considered, as illustrated in Figure 1. A smooth transition in the normal direction of the boundary is achieved by convolving both subdomains over a region with a half-width of ε, which results in the general formulation an arbitrary field variable as shown in Eq. (1).

Figure  1  Subdomains for deriving the BDIM equations and computational domain

Download: Full-Size Img

where b_f and b_s represent the control equations for an arbitrary physical variable φ in the fluid and solid subdomains, respectively. μ^ε₀ (Schlanderer et al., 2017) is a dimensionless smooth transition function, with values ranging from 0 to 1, as represented by Eq.(2).

For fluids, the compressible Navier–Stokes equations are numerically resolved to simultaneously investigate the acoustic and hydrodynamic behaviors as follows.

where the shear stress is denoted by τ_ij, the fluid velocity by u_i, density by ρ, and static pressure by p. The following relationship holds:

with the sound velocity c_T set as 1 500 m/s in the present underwater simulation. γ is the specific heat capacity coefficient of water (γ = 6.1), and B denotes the stiffness coefficient of water, which shares the same dimension as pressure p (B = 345 MPa). By substituting Eq. (5) into Eq. (3), an equation concerning the pressure term can be deduced from the continuity equation and the state equation.

For solids, velocity and pressure are naturally decoupled, where the velocity u_b can be solved directly by preconditions. The governing equation of pressure in the solid body can be represented as follows:

2.2   Numerical methods for solving the compressible Naiver–Stokes equations

A finite difference method is used to discretize the fluid equations. As shown in Figure 2, field variables are defined at the nodes of the grid cell. The convection term on the right-hand side of Eq. (4) is computed using a secondorder upstream difference scheme, as shown by streamwise velocity in Eq. (8). This scheme considers the influence of two nodes on the adjacent side of the control node through Eq. (9). Thus, setting values for the discrete points in the outermost two layers during initialization is necessary, with the secondary outer grid nodes obtained by the first-order upwind scheme. The diffusion term utilizes the second-order central difference scheme, and a first-order upwind scheme is applied for time advancement.

Figure  2  Sketch of uniform Cartesian grid cell

Download: Full-Size Img

Figure 3 depicts a chain decomposition method for solving the compressible equations. The pressure field is first explicitly predicted by Eq. (10). Then, the velocity field is updated using the pressure field obtained from the differential equation outlined in Eq. (11).

Figure  3  Procedure for solving the compressible Naiver–Stokes equations

Download: Full-Size Img

where the superscript n represents the current time of the physical quantity, and n + 1 represents the next time. The forces acting on the solid boundary are calculated from the pressure values on the surface. Given that the values within the smoothing region are not physical, the surface quantities were evaluated with a distance of ɛ from the actual surface, which is expected to remain unaffected.

2.3   Computational geometry and mesh

A computational domain was established around a 10-blade rotor to numerically investigate the hydrodynamic and acoustic fields of a high-speed rotating propeller. The rotor has a diameter of 20.32 cm and a constant chord length of 2.54 cm. This setup aligns with Sevik's research (Sevik, 1973). The focus is on the ability of the program code to interpret incoming turbulence on the broadband unsteady force. A turbulence-generating grid with a spacing size of 10.16 cm was used to create inflow turbulence. The velocity inlet was positioned 198 cm upstream of the grid, with the inflow speed set at 4.57 m/s and the rotation speed at 18.43 rps due to limited computing resources. No-slip walls were applied to the top, bottom, and sides of the domain.

Figure  4  Computational domain and boundary conditions

Download: Full-Size Img

The total number of the most refined resolution level is approximately 1.6 billion (1 600 × 320 × 320) in the domain with more than 30 nodes covering the diameter of the propeller, as shown in Figure 5. The smoothing region halfwidth was set to ɛ = 2∆x = 2∆y = 2∆z. Notably, the number of grids cannot meet the requirements of DNS. Thus, it belongs to ILES (Len and Fenando, 2007). No explicit function was applied for the SGS tensor B_ij, and the numerical dissipation was considered sufficient to mimic its action. The near-wall normalized spatial spacing was y⁺ = u_τΔy/ν ≤ 1.25, where u_τ is the friction velocity on the wall, Δy is the normal distance from the first mesh point to the solid surface, and ν is the kinematic viscosity of the liquid. The simulation was performed with a time step of 1×10⁻⁵ s and less than 20 iterations per step, which leads to a maximum Courant number, Co = (U · n)Δt/Δx ≈ 0.01.

Figure  5  Computational domain mesh of axial section

Download: Full-Size Img

3   Results and discussion

3.1   Validation of the numerical results

According to classical theory for broadband extraction algorithms, the thrust spectrum depends on the turbulence level, integral scale of turbulence, and geometrical dimensions of the rotor. The turbulence characteristics downstream of the turbulence-generating grid were first investigated. These characteristics include the turbulence level $\varepsilon=\sqrt{1 / n \cdot\left(\sum\left(u_{x i}-\bar{u}_x\right)\right)}$ and turbulent length scale (L = $\int\left(R\left(u^{\prime}, \tau\right) / R\left(u^{\prime}, 0\right)\right) \mathrm{d} t$) based on the velocity series extracted from Figure 6(b) of the ILES and the autocorrelation functions illustrated in Figure 6(c), which are associated with Taylor's frozen turbulence assumption. The computed turbulence level was 5.4%, and the turbulent length scale was 2.08×10⁻⁹.

Figure  6  Turbulence level and turbulent length scale estimate extracted from the simulation

Download: Full-Size Img

Considering that pressure fluctuations on the rotating blades and their surface integral are the dominant sources of radiated noise, the predicted thrust spectrum was compared with test results for validation, as depicted in Figure 7 (Hu et al., 2023). Notably, the broadband thrust spectrum was processed using a Fourier transform and the Savitzky– Golay (Zimmermann and Kohler, 2013) approach based on the principle of least squares. This method accurately reflects the fluctuating amplitude of the unsteady force and eliminates random components. The amplitudes were converted into dB units using Eq. (12) due to the different orders of magnitude for various frequency components.

Figure  7  Comparison of predicted and experimental broadband force spectral characteristics

Download: Full-Size Img

where F_x(N) represents thrust in Newtons and F_x(dB) represents thrust in Decibel (dB). Apparently, it shows a similar trend to the experimental thrust spectrum in the low-frequency region, with the "humps" predicted to be occurring near the blade-passage frequency, which is 184.3 HZ in the present simulation. Therefore, the accuracy of the unsteady simulation results over the concerned frequency range is essentially validated.

To assess the grid independence of the simulation, three additional sets of coarser resolution mesh, as presented in Table 1, were adopted, and the formal order of convergence was calculated. The Courant number was kept constant, and the mesh refinement factor was 1.2. Considering the lack of an analytical solution for this problem, the data from the highly-resolved case presented in Figure 6(b) were used as a reference.

The mean velocity data obtained at the probe point in Figure 6, when the vortex structure fully developed, was adopted for the analysis. Figure 8 presents its L₂ norm error over the different grid resolutions, which was considered a measure of the global error made by the boundary scheme. Approximately second-order convergence can be found.

Figure  8  Convergence of the velocity field for the flow of the turbulence immersed propeller

Download: Full-Size Img

3.2   Unsteady flow structures and noise sources

Figure 9 illustrates the vortex structures identified by the Q criterion (Q = 100), pressure, and streamwise velocity near the rotating blades when the turbulence reaching the propeller is fully developed. According to the views of the incoming flow of the propeller at planes A and B in Figure 9(a), a higher ingesting turbulence level is expected due to pronounced vortex shedding from the grille. These slender vortices rejoin and dissipate until they encounter the rotating propeller. The vorticity contour downstream of the propeller disk is subsequently depicted from plane C to plane D. The tip vortex with its elongated shape can be observed in Figure 9(b). It interacts with the boundary layer of the model hull and serves as a noise source. In addition, significant pressure differences inside and outside the vortical region with concentrated lower velocity are shown in Figures 9(c) and (d). Figure 10 shows the power spectral density (PSD) of the pressure fluctuations at the same probe point of the velocity monitor used to investigate the inflow characteristics. Given that the probe point is far upstream of the hub, no obvious discrete tones appear. The slope is close to the theoretical f^−5/3 scaling in the spectra, which suggests that the resolved scales in the present case are asymptotic to the inertial subrange (Ji et al., 2015). Moreover, the mesh is sufficiently fine to solve the major proportion of the energy.

Figure  9  Unsteady structures around the propeller

Download: Full-Size Img

Figure  10  PSD analysis of pressure fluctuations at monitor positions

Download: Full-Size Img

An additional advantage of the compressible BDIM simulation is identifying the location of the noise source and tracking the propagation of acoustic waves. Figure 11 presents the instantaneous dilatation contours to provide an initial insight into the acoustic field. Upstream of the turbulence grille, the noise level is quite low due to the uniform inflow condition. As turbulence develops in the approach stream, a substantial increase in noise becomes evident. As expected, the acoustic field is dominated by an axial dipole radiation pattern from the propeller, with higher dilatation levels in the x-direction for z/L = 0.5 and y/L = 0.5 sections. Conversely, noise radiated along the axis perpendicular to this direction diminishes significantly.

Figure  11  Instantaneous dilatation contours at an arbitrary instant of time

Download: Full-Size Img

3.3   Effect of turbulence characteristics on the flow field

Another turbulence-generating grid upstream of the rotating propeller with a larger interval size of 15.24 cm, as shown in Figure 12, was simulated to further investigate the effect of inflow turbulence characteristics on the broadband force spectrum. The computed turbulence level decreases to 4.1%, and the turbulent length scale increases to 3.07 × 10⁻⁸ based on the monitored velocity. The difference in scale can be observed by tip vortices extracted by the Q iso-surface in Figure 13. Figure 13 also illustrates the linear distribution of streamwise velocity and mean velocity in the x–y plane through the probe point. A larger axial velocity near the rotor blade disk under the disturbance downstream of Grid 2 is observed, which is expected to show a lower-level thrust magnitude due to the increased dominance of the axial velocity component on the local blade angle of attack, denoted by θ-β, where θ is the fixed angle of pitch, and β is the increased advance angle.

Figure  12  Computational domain for the propeller downstream of Grid 2

Download: Full-Size Img

Figure  13  Velocity analysis of the flow upstream of the propeller

Download: Full-Size Img

The time histories of the pressure at the aforementioned probe point were extracted to compute the PSD for a better understanding of the injected structures in Figure 14. The results show that the lower turbulence level and larger length of scale under Grid 2 lead to a lower amplitude PSD, which means less energy input into the rotating blade. Furthermore, the broadband spectra were compared through the same processing in Section 3.1, as shown in Figure 15. They exhibit the same trend, with the main frequency located near the 1BPF, which results in a lower amplitude and, thus, a lower noise level.

Figure  14  Comparison of PSD analysis of pressure fluctuation under different grids

Download: Full-Size Img

Figure  15  Comparison of thrust spectra under different grids

Download: Full-Size Img

4   Conclusions

The unsteady compressible flow around a rotating propeller was simulated with the BDIM framework, with simultaneous visualization of the near-field flow field and sound propagation through velocity dilatation. The main conclusions can be summarized as follows:

1) The algorithm developed based on the BDIM effectively captures the unsteady thrust spectrum of the rotating propeller under approach turbulence, which displays a broadband hump near 1BPF.

2) The study offers insights into the interaction between inflow turbulent structures and the formation of intense tip vortices around the high-speed rotating propeller. In addition, the dipole formation wave front along the axial direction demonstrates that the compressible BDIM accurately represents noise radiation from flow-induced noise generation due to the scattering of the rotating blade.

3) The increased energy attenuation from more widely spaced turbulence-generating grids corresponds with a reduction in the amplitude of pressure fluctuation and an increase in velocity profiles near the rotor disk. Moreover, the interaction between the upstream, which is characterized by a lower turbulence level and a higher turbulent length of scale, and the rotating propeller results in lower broadband components in force spectra and a slight increase in the scale of tip vortices.

The attempts reported in this paper have important implications for noise mitigation in underwater vehicle design. Future research in this field could focus on establishing advanced noise reduction strategies and further refining computational approaches to gain comprehensive insights into underwater acoustics. These advancements could lead to the creation of quieter and more environmentally friendly underwater vehicles for addressing crucial issues in marine acoustics and ecology.