Sound Field Properties of Non-Cavitating Marine Propellers

The sound field properties of non-cavitating marine propellers are investigated using a hybrid method, in which the FWH (Ffowcs William-Hawkings) analogy is coupled with the BEM (Boundary Element Method) approach. The investigations include both the uniform and non-uniform inflow conditions. For both conditions, the dominant sound source terms and the decay rate of the noise with regard to the distance to propeller centre are investigated. The influence of the permeable surface dimensions in the permeable FWH approach on the hydroacoustic result is also investigated. To carry out the investigations, the formulation to calculate acoustic pressure generated by the propeller wake sheet is proposed for the first time. The issues associated to coupling permeable FWH approach and BEM are also discussed, including the fictitious volume flux problem and the consideration of the ship wake field. It was found that the influence of the permeable surface dimension is little for the 1st BPF (Blade Passage Frequency), but cannot be ignored for the 3rd BPF. In the uniform inflow situation the thickness terms are found to be dominant, while in the non-uniform inflow situation the loading terms are dominant.


Introduction
Human behaviour induced underwater noise has attracted public attention recently. It raises the underwater background noise level by a significant value [1] and very likely causes damage to marine life [2]. Shipping, or more specifically, the operating marine propeller, is recognized as one of the main sources of underwater radiated noise. Thus, different research projects aim to understand and reduce the radiated noise of marine propellers, for example, the European Union (EU) level projects AQUO, SONIC, and ProNoVi [3][4][5].
Numerical methods have already achieved popularity for the analysis of marine propeller hydrodynamics and for the design procedure. Also, for problems related to underwater radiated noise, numerical simulations are expected to help with understanding the mechanism and improving the design of marine propellers. While many experiences have been obtained in the computational aeroacoustics [6], we still need to address the following two questions in the maritime community: (1) How suitable are the existing approaches and best practices for marine propeller applications? (2) What is the optimal set-up to acquire reliable hydroacoustic predictions with feasible computing effort?

State of the Art
Most of the existing studies in the literature are devoted to addressing the two questions mentioned above, that is, they are focused on the validation of the adopted approaches and the finding of best practices for the marine propeller problem, at different levels and aspects. The most popular calculation approach is the hybrid method, combining an incompressible hydrodynamic solver and the FWH (Ffowcs William-Hawkings) acoustic analogy [7]. The P-FWH (permeable FWH) approach proposed by Di Francescantonio [8] is especially frequently used. The direct FWH approach, in which the Farassat 1A formulation [9] is applied to body surfaces and the non-linear quadrupole term, calculated by the direct volume integration, was also adopted in some work for comparison purposes.
The existing numerical investigations can be divided into three categories according to the object being calculated, that is, ship-propeller system, single propeller, and more fundamental sound sources. The simulations of ship-propeller systems mainly aim to gain confidence in acoustic analogy among the maritime community and the results are mostly validated against experimental measurements. Ianniello et al. [10] calculated the radiated noise of a complete scaled ship model by coupling an incompressible RANS solver and the FWH analogy. The FWH approach with direct volume integration and the P-FWH approach were both investigated. The permeable surface was set to envelop the whole ship. The obtained acoustic pressures were compared to the RANS pressures and good correlations were observed. Bensow and Liefvendahl [11] computed the radiated noise of a cavitating propeller behind the research vessel the Princess Royal. The incompressible implicit LES and P-FWH approaches were combined, in which an open permeable surface around the propeller was used for data sampling. When compared to model-scale measurements, the numerical results showed a great underestimation for frequencies less than 400Hz and a good agreement above 400Hz. Li et al. [12] used DDES (Delayed Detached Eddy Simulation) and the P-FWH approach to compute the radiated noise of a full-scale ship and compared the results with sea trial measurements. The permeable surface was placed around the ship. The obtained pressure pulses agreed well with measurements at the first six orders of BPF (Blade Pass Frequency). Significant underestimation of the broadband noise was reported for frequencies between 50 and 100 Hz and also about 200 Hz. Göttsche et al. [13] computed the Princess Royal case by combining the BEM (Boundary Element Method) and the FWH analogy. The Farassat 1A formulation was adopted for the FWH analogy. A general good agreement between the numerical results and the full-scale measurements was achieved. These works showed that reasonable results for the marine propellers' noise can be obtained with acoustic analogies. However, further investigations are necessary to improve the overall prediction accuracy.
The simulations of single propellers aim to find best practices for hydroacoustic calculation and gain a physical understanding of the main sound sources. The open water cases are mostly investigated. Due to the lack of experimental data, the obtained acoustic pressures are frequently compared to the pressures obtained in hydrodynamic solvers. Using BEM and the Farassat 1A formulation, Seol [14] analysed the property of the sound field radiated by a marine propeller without cavitation and with sheet cavitation. For the non-cavitating case, a clear dipole directivity pattern was identified. The sheet cavitation was recognized as the dominant noise source once it occurs. Lloyd et al. conducted several calculations for a single propeller in uniform inflow with a RANS and P-FWH approach. The implementation was validated by good correlations of the obtained results with the corrected measurement data [15]. Then the approach was used for the simulation of the E779A propeller in Reference [16]. The pressure obtained with the hybrid method and RANS solver were compared. It was found that the hybrid method was more sensible to the grid refinement. Numerical disturbances were also observed for the hybrid method, especially at high frequencies and with a far distance from the propeller. They investigated the influence of the end-cap of the permeable FWH surface in another work [17] for the non-cavitating open water case. In this study, the shedding vortices penetrated through the end-cap. It was concluded that closed end-caps leaded to obvious numerical disturbances, while open end-caps resulted in wrong modelling of the monopole contribution, the effect of which is not ignorable. Lidtke et al. [18] coupled the URANS and P-FWH approaches to compute the tonal blade passage noise of the PPTC propeller and coupled LES and P-FWH approaches for the noise generated by a hydrofoil. It was concluded that RANS was unable to accurately account for cavitation dynamics and the associated noise, and LES might help to gain detailed insight into the low-frequency noise generation mechanisms. It was also stated that the incompressibility assumption and the ignorance of fine-scale bubble dynamics would weaken the accuracy. Testa et al. [19] investigated the effectiveness of BEM hydrodynamic data for propeller hydroacoustics by comparing the pressures obtained with the BEM/P-FWH hybrid method, DES/P-FWH hybrid method, and the DES solver for the open water case. Observers were located along the axial direction above the propeller. It was found that the BEM/P-FWH approach worked well for the observers upstream and a little downstream (1D) of the propeller and failed to capture the total sound signals for further downstream positions. The reason was ascribed to the neglect of turbulence induced noise.
Other basic geometries or virtual objects are also studied in the maritime community to gain more understanding of the sound generation mechanism and acoustic analogies. Ianniello [20] developed a concept of rot-pole for the analytical study of rotating sound sources. The properties of the resulting sound field have been analysed, including the dependence of the decay rate on the Mach number and the number of source points, which is a representation of the propeller blade number. It was found that the decay rate out of the rotating disk plane is one order faster than in the disk plane. In the disk plane, the decay rate is firstly r −2 in the near-field and then r −1 in the far-field. Increasing the number of source points leaded to a faster decay rate in the near-field, but did not affect the decay rate in the far-field. Leifvendahl and Bensow [21] investigated the flow-generated sound around different cylinders with both numerical simulation and semi-empirical model. For the numerical simulation, LES and Curle acoustic analogy were coupled together. In-depth analysis of the flow field regarding the sound generation was conducted, for example, the spanwise correlation of the vortex structures. The obtained sound predictions showed reasonable agreement with measurements. Cianferra [22] derived a numerically applicable formula for the volume integration of the quadrupole term in the FWH equation. With the flow field obtained using LES, he compared the flow-generated noise around a square cylinder obtained with the direct FWH and P-FWH approaches. Multiple permeable FWH surfaces have been tested. It was concluded that the direct FWH approach gave the most accurate results and the computational effort was acceptable. It was also found that the persistence of a 2D shaped wake leads to more radiated noise when compared to a 3D wake.

Contributions of Current Work
With the current work, we contribute to two aspects: (1) The appropriate use of the P-FWH approach, especially the justified placement of the permeable surface. In the P-FWH approach a flow domain is enveloped by the permeable surface. When an incompressible solver is adopted for the hydrodynamic calculation, as usually be done, then inside this domain the sound speed is effectively infinite large due to incompressibility. The dimension of this incompressible domain would affect the acoustic results, and cannot be arbitrarily chosen. This effect will be studied by comparing the results obtained with direct FWH approach and P-FWH approach with different permeable surfaces; (2) Understanding the decay rate and main sound source terms of propeller radiated noise in both near-field and acoustic far-field. The distance of observers to the propeller centre would vary from less than one diameter (5 m) to several acoustic wavelengths (550 m). To the authors' knowledge, such a study has not been carried out in the literatures.
Our in-house BEM solver panMARE is coupled with the acoustic analogy for the investigations. A non-cavitating propeller operating in both the uniform (open water case) and non-uniform inflow (behind-hull case) would be analysed. To achieve the research objectives, other relevant issues have also been addressed in this manuscript. These include: (1) The formula to calculate acoustic pressure induced by the potential wake sheet, which models the vorticity flow after the propeller. Such a formula is presented for the first time; (2) The appropriate coupling of the P-FWH approach with BEM for non-uniform inflow conditions. This is found not to be a trivial activity.
Although BEM is unable to capture as many flow details as the finite volume CFD methods-for example, RANS, DES or LES-it has advantages over these methods when the effect of the incompressible domain dimension in the P-FWH approach is investigated. In BEM, there is no 3D mesh in the flow domain, but only panels on the body and wake sheet surfaces. The velocity and pressure in the flow domain are calculated analytically using the solution, that is, velocity potential and its normal derivative, on the panels. Thus, the position of the permeable FWH surface is relatively more flexible. In contrast, in finite volume methods, the permeable FWH surface must be placed in the fine grid structure to prevent the information loss caused by numerical diffusion. Besides, there is also no vorticity distribution in the flow domain in BEM. This ensures that permeable FWH surfaces with different dimensions would have the same sound sources. In contrast, with finite volume methods, a larger permeable surface take into account of more non-linear sound sources than smaller ones. That would make it difficult to interpret the possible differences between the results. Furthermore, in the current work we mainly analyse the acoustic signals in the propeller disk plane at the first BPF (Blade Passage Frequency) for the investigation. It has been shown by other works [13,19] that BEM is capable of providing reliable hydrodynamic data for the hydroacoustic calculation under these circumstances.
The rest of the manuscript is organized as follows-the mathematical formulas, term definitions, and important numerical methods are described in Section 2. In Section 3, numerical results for different cases and the corresponding analysis are presented. Finally, the conclusions are summarized in Section 4.

Methodology
A very brief description of BEM is given first in Section 2.1 to help understand the coupling of BEM and FWH analogy better. Then, the used formulations of the FWH equation are given or derived in Sections 2.2-2.4.
In the current work, the relative positions between the observers and the ship are fixed, that is, the observers move forward together with the propeller in the axial direction. Although this set-up may not be consistent with full-scale experiment, this is a popular treatment in numerical simulations. Part of the formulations in this section are dependent on this assumption.

Boundary Element Method
In BEM, the flow is assumed to be inviscid and irrotational. Thus, the potential theory is used and the main unknown is the disturbed velocity potential φ. Knowing its values and normal derivatives on the boundaries, one can obtain the deterministic distribution of φ in the fluid domain. The potential's jump (equivalent to potential value when φ is assumed to be zero inside the body) and normal derivative are represented by doublet and source strength on the boundaries, and they are the variables to be solved for. The domain boundaries consist of the solid body boundary and the wake sheet. The wake sheet geometry is determined through an iterative process called wake alignment [23], which requires the wake sheet being locally tangent to the fluid velocity. The boundaries are discretized into panels to enable a numerical solution, as shown in Figure 1. There is source and doublet on the solid body panels and only doublet on the wake panels. In the current work, our in-house low order BEM code panMARE is used, which means the source and dipole strength are constant on the panel.
In open water cases, the doublet strengths are constant on each streamwise wake panel strip and are determined according to Kutta condition [25]. In behind-hull cases, the doublet strength changes along the streamwise direction due to the time-varying load on the blades. In current work, the source strengths are pre-evaluated according to impermeable boundary condition. The doublet strengths are obtained by solving a system of linear equations being built with the Dirichlet boundary condition, that is, φ = 0 inside the body [26]. Knowing the source and doublet strengths, the velocity distribution is obtained by adding up the induced velocities and the undisturbed inflow velocity. The pressure in the fluid domain is calculated with the unsteady Bernoulli equation.
More detailed descriptions of BEM can be found in References [23,25,26].

Ffowcs Williams-Hawkings Equation
The current work is based on the FWH equation in the Farassat 1A formulation [9], which is where p (x, t) is the sound pressure at the observer position x and time point t, p T and p L are called the thickness and loading terms, respectively; f = 0 denotes the integration surface for FWH equation (or sound source surface), ρ is the fluid density, ρ 0 is the undisturbed value, c 0 is the sound speed, and p is the absolute pressure; u is the fluid velocity, v is the surface moving velocity, both are defined in a coordinate where the undisturbed fluid (or far-field fluid) is static; M is the vectorial Mach number, that is, M = v/c 0 ; r is the vector from the source point to observer position, r is its length, the subscript ret means that variables inside the square brackets are computed at the corresponding retarded times, that is, at τ = t − r/c 0 ; the subscripted variables denote dot products of a vector and a unit vector implied by the subscript, except for L M = L · M, and n is the unit normal vector on the surface directing into the fluid domain; the dot over a variable means time derivative.
In the Farassat 1A formulation, the quadrupole component-due to the Lighthill stress tensor-is ignored. This ignorance is especially justified in current study, because the Lighthill stress tensor vanishes in an incompressible potential fluid domain. In Section 2.3, we will show that the effect of shedding vorticity can be considered with surface integrals on the wake sheet.
As the observer moves with regard to the static fluid, its position change between the sound generating time τ and the receiving time t must be considered. This is done by calculating the vector r using the formula provided by Bres et al. [27], that is, with and where ∆x, ∆y and ∆z are components of the vector from source point to the observer position at the time of sound generation, and u 0 is the observer's velocity along the x direction. The motion in other directions would not be considered in the current study.
To facilitate the detailed analysis of the components, the thickness and loading terms given in Equations (2) and (3) are further decomposed as where Equations (1)-(8) were directly used for the P-FWH approach. When the formulation is applied to solid surfaces, further simplification can be achieved. The impermeable boundary condition implies u n = v n , which leads to the simplification of the terms U n and L on solid surfaces as

Formulation of FWH Equation on the Potential Wake Sheet
In BEM, the wake sheet is used to represent the shedding vorticity of the propeller blades. The wake sheet has zero thickness and acts as a potential boundary in the solution process, as shown in Figure 1. In this section, we will show that this zero-thickness surface representation of the vorticity field turns the volume integration of the Lighthill stress term into a surface integration on the wake sheet. It should be noted that the following derivation is based on the low-order BEM solver. However, the basic idea would be similar for higher order BEM solvers.
The evolution of the wake sheet can be considered from two perspectives, or two coordinate reference systems, that is, the earth-fixed global reference system and the propeller-fixed local reference system. To make understanding easier, the wake alignment is not considered in the derivation. In the earth-fixed perspective, the propeller rotates and moves forward, and new wake panels are shed out at the trailing edge in every time step. The existing wake panels are static and the doublet strengths on them do not change with time. The wake panels being far enough from the blade are removed to reduce computational cost and mimic the vorticity diffusion. Normally a wake length will be defined to limit the maximum wake panel numbers in the streamwise direction. In the propeller-fixed perspective, both the propeller and the wake panels stay static. The doublet strength is convected panel-wise downstream. Therefore, the doublet strength on the wake panel is time-varying. These two perspectives are depicted in Figure 2. Now we derive the formulation of the FWH equation on the wake sheet in the propeller-fixed perspective. It should be noted that the velocity defined in the earth-fixed reference system will be finally used to calculate the FWH terms. A closed surface is applied to envelop the wake sheet. The surface consists of two subsurfaces parallel to the wake sheet, that is, the upper subsurface and the lower subsurface on both side, and end cap subsurfaces to close the envelope, as shown in Figure 3. Let the distances between the two parallel subsurfaces to the wake sheet approach zero, then they have the same positions but opposite normal directions. The areas of the end cap surfaces become zero and their contribution to the acoustic pressure vanishes. The gap between the upper and lower surface is set to be infinitely close to zero.
We now discrete the upper and lower subsurfaces in the same manner with the wake sheet. Then every FWH term will be checked on the panel. It will be shown that thickness terms cancel out each other, while the loading terms do not totally vanish and form the FWH formulation on the wake sheet. We use the superscript '+' to denote the quantities on the upper subsurface and '-' for those on the lower subsurface. No superscript means the summed effect.

Thickness Terms
Because there is no potential source on the wake sheet, the normal component of the fluid velocity is continuous across the wake sheet. Due to the opposite normal directions on the upper and lower subsurfaces, the following relation holds: Since the upper and lower subsurfaces share the same position with the wake sheet, the surface moving velocities are identical, that is, With these two relations, the following equation can be deduced andU n + Uṅ = dU n /dt ≡ 0.
These imply that all the thickness terms vanish on the wake sheet, that is,

Loading Terms
With the above-mentioned velocity relation, the following equation can be obtained For the first term on the right-hand side, the pressure is calculated with the unsteady Bernoulli equation, that is, where U ∞ is the undisturbed relative inflow velocity, u ind is the induced velocity, φ is the perturbation velocity potential, and p 0 is the far-field pressure. With Equation (23) it can be obtained that The difference (u + ind − u − ind ) is related to the vortex density on the surface. It can be calculated as the opposite of the doublet strength's surface gradient, that is, where µ is the doublet strength on the wake sheet, and ∇ s denotes the surface gradient operation. The velocity U ∞ + (u + ind + u − ind )/2 is the fluid velocity without considering the local dipole gradient, and we denote it as u 0 . Considering the relationship between doublet strength and potential jump on the wake panel, that is, Combing all above together leads to For the second term on the right-hand side of Equation (22), in a real flow the vortices always moves with the fluid velocity, which means (u n − v n ) = 0 on the wake sheet. Thus, the second term would vanish.
Combining the results, we have This can be substituted into Equations (14)- (16) for the calculation of acoustic pressures. In the open water case, u 0 is perpendicular to ∇ s µ andμ ≡ 0, so the wake sheet does not contribute to the acoustic pressure. For the behind-hull case, the term does not vanish and will generate acoustic pressure.
In the earth-fixed perspective, because of the non-time-varying doublet strength on the wake panels, the same analysis will lead to different results. However, there is a dilemma about the newly generated wake panels. It implies a local growth of the integration surface, for which the Farassat 1A formulas are invalid, because the surface divergence term in the original work of Ffowcs Williams and Hawkings ([7], Equation (4.3)) has been ignored or regarded as unity. When the focus is put locally on the newly generated wake panels, the area of the integration surface grows from zero to a finite value, which means an infinite surface divergence. Then the integration on the new panels is undefined. Such a dilemma does not happen in the propeller-fixed perspective.

Coupling of Permeable FWH Approach with BEM
When a permeable surface enveloping the propeller and wake sheet is used as the sound source surface, the Equations (1)-(8) can be applied without any modification. The permeable surface would move forward together with the propeller. All quantities should be defined in the static flow, thus u is the fluid velocity in the earth-fixed reference coordinate system, and v is the constant forward velocity.
There are two specific issues when the P-FWH approach and BEM are coupled for the hydroacoustic calculation for the behind-hull case. The first issue is whether the ship wake velocity should be considered in u. If so, how it should be included. Normally, ship wake velocities are only provided on the propeller disk. This information is not enough to obtain the ship wake velocity on the whole permeable surface. At the same time, because of the nonlinear relation between pressure and local fluid velocity, the addition of even a constant ship wake velocity will lead to differences in all the loading terms. The second thickness term p T,2 is also influenced, while the first thickness term p T,1 is an exception. As the ship wake velocity do not changed with time in BEM solver,U n is not affected. If we set up the permeable surface so that the surface normal vectors do not change with time, then Uṅ will also not be affected. Under the two above mentioned conditions, the term p T,1 keeps the same whether the ship wake velocity is included or ignored on the permeable surface. It will be shown in Section 3.3 that p T,1 is, fortunately, the dominant term when the P-FWH approach is used. Thus it becomes unimportant whether or not to include the ship wake velocity in u.
The second issue is caused by using ship wake velocity on a 2D disk as a 3D distribution in BEM. The ship wake velocity is assumed to be independent of the axial position. This canonical treatment is although proved to work well for the hydrodynamic calculation, it causes problem in the hydroacoustic calculation when the P-FWH approach is used. The velocity field produced in this way is not divergence-free. The source strength on the blade surface, which is calculated according to this velocity field, would therefore not integrate to zero. This will be manifested by a non-negligible volume flux on the permeable surface, which can be calculated as where Q denotes volume flux, and f = 0 denotes the permeable surface. This fictitious volume flux depends on the propeller's azimuthal position and is thus time-varying. It will generate pure numerical acoustic pressures and contaminate the results. A simple but effective solution is adding a virtual source point at the propeller centre to cancel out the fictitious volume flux. The fictitious volume flux is calculated in every time step, and its opposite is used to set the strength of the virtual source point, that is, Then the induced velocity of the virtual source point is superposed on existing velocities on the permeable surface. That is to say, on every panel of the permeable surface, the following calculation is carried out, where r vir is the vector directing from the virtual source position (propeller centre) to the permeable surface panel centre. Such a virtual source correction procedure ensures a significant reduction of the fictitious volume flux. It should be noted that the virtual source is not used for the hydrodynamic calculation, but only used to correct the velocity on the permeable surface. Thus, it has no influence on the hydrodynamic field.

Numerical Results
In this section, the open water case and behind-hull case are analysed with the direct FWH 1 approach in Sections 3.1 and 3.2, respectively. The hybrid method combining BEM and P-FWH approach is investigated in Section 3.3. The comparison of the sound pressures obtained with the direct FWH and P-FWH approaches and the influence of permeable FWH surface dimensions on results are also described in Section 3.3. Finally, in Section 3.4 we present a numerical verification of the formulation to calculate the sound produced by the wake sheet.
The full-scale KP505 propeller is used in the simulations. KP505 is an open access five-blade fixed pitch propeller designed for the MOERI container ship (KCS) [28]. Its main particulars are listed in Table 1. For both the open water and behind-hull cases, the sound pressure decaying behaviour and the dominant FWH terms at different distances are investigated numerically. For such purposes, a set of observer points are placed in the propeller disk plane. They are located starboard on a horizontal line which having the same depth with the propeller centre. The distance to the propeller centre varies from 5 m up to 550 m. The largest distance 550 m is around 70D or 2.3λ 1 , where D is the propeller diameter and λ 1 is the wavelength for BPF (Blade Passage Frequency). In the calculation, both the open water case and behind-hull case are simulated in the unsteady manner. The hydrodynamic simulation is firstly carried out until convergence, and thereafter the hydroacoustic calculation is turned on. The time-step in the simulation corresponds to a rotation of 6 degrees. The numerical result will be mainly presented for the 1st BPF.

Open Water Case with the Direct FWH Approach
The full-scale KP505 propeller is analysed at the design point. The SPL (Sound Pressure Level) at the BPF are plotted in Figure 4 as a function of the distance to the propeller centre. One curve denotes the SPL calculated with the BEM/FWH hybrid method, the other denotes the result calculated in BEM solver using the unsteady Bernoulli equation. Both results correlate very well with each other in the near-field, and the deviation exists in the far-field. The deviation is caused by the incompressibility assumption in the BEM solver, which has a significant effect on the SPL in the far field.
For further verification, the temporal pressure variations at a near-field observer obtained with the FWH analogy and the BEM solver are compared in Figure 5. The good correlation verifies the current implementation of the coupling of the FWH analogy and BEM solver.
To obtain more details, the contributions of different FWH terms are plotted as a function of the distance to the propeller centre in Figure 6. The thickness terms are always dominant over the loading terms. The term p T,2−2 is dominant in the near-field and p T,2−1 is dominant in the far-field. It can also be observed that, the far-field terms p T,2−1 , p L,1 , and p L,3 decay with the decay rate r −Z in the near-field and r −1 in the far-field. The near-field terms p T,2−2 and p L,2 decays with the rate r −(1+Z) in the near-field and r −2 in the far-field. The decay rate is increased significantly in the near-field due to the interference between acoustic pressures from different blades. The summation of all terms result in that the SP (Sound Pressure) obtained with FWH analogy decays very fast, that is, with the rate of r −(1+Z) , in the near-field. Then in the far-field, the decay rate gradually changes to r −1 .
The decay rates obtained here agree very well with the analytical in-plane decay behaviour of the rot-pole, that is, a sound source consisting of one or more rotating force sources, which is proposed and analysed by Ianniello [20]. Figure 24 and 27 in Reference [20] show a decay rate of r −(1+Z) in the near-field and r −1 in the far-field for both the total thickness term and total loading term.

Behind-Hull Case with the Direct FWH Approach
For the behind-hull case, the propeller operates in a ship wake field. The loads on every blade will be different and vary with time. The potential wake sheet acts also as sound sources and the corresponding acoustic pressures are calculated using the formula developed in Section 2.3. In this section, the SPL decay rate, the influence of the wake sheet length on the sound field, and the dominant FWH terms in both blade and wake sheet induced sound signals will be investigated. The ship wake field being used in current work is shown in Figure 7. It is a relatively smooth wake field with a wake fraction of w = 0.16. The propeller advancing speed V s is set to make the effective advance ratio being The SPL obtained with different wake sheet lengths are given in Figure 8, in the form of separate contributions from blade surface, wake sheet and the combined results. There is no obvious range in which the SP decays very fast, that is, with the rate r −(1+Z) as for the open water case. In the far field, both the blade induced and wake sheet induced SP have a decay rate of r −1 . The blade is the dominant sound source. However, such an observation is highly dependent on special cases and wake sheet lengths. The right plot of Figure 8 shows that, the wake sheet produces a 2 dB difference in the total SPL when the hydrophone locates more than 30 m from the propeller centre. wake 4rev 2 2 total 2rev 2 2 2 wake 2rev Figure 8. SPLs in the behind-hull case calculated with different wake sheet lengths. The right plot is an enlarged view for the region 20 m~100 m. In the legend blade means only sound produced by the blade, wake means only sound produced by the wake sheet, total means the combination of both, and 2 rev~8 rev denote the length of wake sheet in revolutions. Figure 9 shows the contributions of different FWH terms in the blade induced SPL. In the behind-hull case, the loading terms p L,1 and p L,2 are evidently dominant. Unlike in the open water case, in the behind-hull case p L,1 and p L,2 achieve the final decay rates of r −1 and r −2 quite early. The reason is that the pressure and also the time derivative of pressure are not equal on each blade. The loading terms on different blades can not cancel out each other so much as in the open water case. SPL re 1µPa p p T,2 p L,1 p L,2 p L,3 Figure 9. Contributions of different FWH terms in the blade induced sound pressure for the behind-hull case. Figure 10 shows the contributions of different FWH terms in the wake sheet induced SPL for the case with 2 wake revolutions. The thickness terms are always zero and are therefore not shown. The first loading term p L,1 is dominant in the far-field and the second loading term p L,2 is dominant in the near-field. This is similar to the situation on the blade surface, which is shown in Figure 9.

Behind-Hull Case with the BEM/P-FWH Hybrid Method
In this section, the acoustic pressure is calculated with the hybrid method, in which the BEM and P-FWH approach are coupled. The difference between the results obtained with direct FWH and P-FWH approaches and the influence of the permeable surface dimension will be investigated. The wake sheet length is set to 2 revolutions for the calculations in this section. Three cylindrical permeable surfaces with different dimensions are used in the investigation. They are denoted with the names basic, large, and long, respectively. Their size and discretization parameters are listed in Table 2 and depicted in Figure 11. The permeable surface with the name basic is used for the general study, and the others are used for the influence of the permeable surface dimension. The adopted panel sizes listed in Table 2 are chosen through mesh independence study, which will be discussed later. radius front length rear length basic permeable surface large permeable surface long permeable surface Figure 11. Different permeable surfaces and the size parameters. The size parameters are depicted on the basic permeable surface. Two wake revolutions are shown here. Table 2. Size and discretization parameters for the permeable surfaces with names basic, long, and large. The meaning of size parameters are explained in Figure 11. In the following, we first show the fictitious volume flux problem and the performance of virtual source technique described in Section 2.4. Thereafter, sensitivity of the result to the permeable surface's panel size is given. Then, the influence of whether considering ship wake field velocity on the permeable surface is investigated, the relative importances of different FWH terms are also given for this investigation. After that, temporal variations of sound pressure obtained with BEM solver, direct FWH and P-FWH approaches are compared for verification purpose. Finally, the SPL obtained with the direct FWH approach and P-FWH approach using different permeable surface dimensions are presented.

Parameters
According to the analysis in Section 2.4, a non-physical fictitious volume flux may exist on the permeable surface, which is caused by the use of ship wake velocity on a 2D plane as a 3D distribution in BEM. Figure 12 shows the temporal variation of the volume flux on the permeable surface, uncorrected or corrected with the virtual source technique, and the integral potential source on propeller surface. The fictitious volume flux is evident, and its oscillation amplitude is more than 0.2 m 3 . This will act as a fictitious sound source. The volume flux has the same oscillation pattern and amplitude with the integral source strength. This explains that, the fictitious volume flux is caused by the non-zero integral source strength on blades. With the virtual source correction described in Section 2.4, the volume flux is reduced to the order of 1 × 10 −4 m 3 . Figure 13 shows the SPL obtained with and without the virtual source correction, using the basic permeable surface. An evident improvement of the results for almost the whole distance range can be observed. The sound pressure amplitude induced by the volume flux (acts like a oscillating sphere) decays as r −1 . Thus, it leads to the cyan straight line with decay rate r −1 in Figure 13.  The sensitivity of the results to the permeable surface's panel size is investigated with the basic permeable surface. The SPLs obtained with different panel sizes are listed in Table 3. The SPL for BPF at the observer being 550 m from the propeller centre is plotted as a function of the panel size in Figure 14. The obtained SPL converges as the panel size decreases. When the panel size is less than 1.0 m, the change in SPL is less than 0.2 dB. The difference in SPLs obtained with 0.5 m and 0.2 m panel size is less than 0.1 dB. Thus, the panel size 0.5 m (D/16) has been chosen for the basic permeable surface for the calculations.
For other permeable surfaces, the same procedure for the panel size convergence has been carried out. The finally chosen panel sizes for all permeable surfaces are given in Table 2. Generally, a smaller permeable surface requires finer discretization on it.
The SPL calculated with and without considering the ship wake velocity on the permeable surface are plotted in the left-hand panel of Figure 15, as a function of the distance to the propeller centre. The two curves coincide with each other, which implies that whether including the ship wake velocity on the permeable surface has little influence, if any, on the acoustic result. To understand the reason, contributions of different FWH terms on the permeable surface are plotted in the right-hand panel of Figure 15, as a function of the distance to the propeller centre. It can be observed that the first thickness term p T,1 is always dominant. As explained in Section 2.4, with the current configuration, the term p T,1 is not influenced by the ship wake velocity. The other terms will be influenced, but their contributions to the total sound pressure are secondary when compared to p T,1 . This explains the ignorable influence of whether considering the ship wake velocity on the permeable surface on the acoustic results. In all other calculations in this section, the ship wake velocity is not considered on the permeable surface.   Table 2. For further verification of the implementation of the direct and P-FWH approaches, the temporal variation of the acoustic pressure obtained by BEM solver, direct FWH and P-FWH approaches are shown in Figure 16, at the observer which is 8m from the propeller centre. The pressure variations correlate with each other very well. This verifies the correct implementation of both FWH approaches. A little phase difference can be observed between the pressures obtained with BEM and the direct FWH approach, which may due to the time delay caused by compressibility considered in the FWH formulation. To study the effect of the permeable surface dimension, the SPLs calculated with different permeable surface dimensions are compared with those obtained with the direct FWH approach. The results are given in Figure 17 as a function of the distance to the propeller centre, in the left plot for the 1st BPF and in the right plot for the 3rd BPF. For the 1st BPF, the results correlate very well with each other except for the observers locating inside of the large permeable surface. At the observer being 100 m from the propeller centre, the difference in the SPL obtained with different permeable surface dimensions is less than 1.5 dB. For the 3rd BPF, the difference becomes significantly larger. The SPL from the P-FWH approach is around 10 dB larger than the SPL from the direct FWH approach. The results obtained with different permeable surface dimensions also fail to correlate with each other well. The difference can be up to 3 dB. The reason is not clear and needs further investigations. The results here help us to conclude that, when different permeable surface dimensions lead to difference in acoustic results, even combined with viscous solvers, one should not easily ascribe it to the quadrupole terms between them.

Verification of the FWH Formulation on the Wake Sheet
The formula to consider the acoustic contribution of the propeller wake sheet in Section 2.3 is proposed for the first time. Here we show a kind of verification for the proposed formula with help of the P-FWH approach. The sound generated by two revolutions of wake sheets are calculated with both the newly proposed direct formula in Section 2.3 and the P-FWH approach. The SPL for the BPF from both methods are plotted in Figure 18 as a function of the distance to the propeller centre. The propeller blades are excluded during this calculation. The permeable surface with the name basic in last section is used. Except for some deviation in the very near field, the results correlate very well with each other for observers locating beyond 16m (2D) from the propeller centre. This good correlation verifies the correctness of the proposed formula, at least for the calculation of the acoustic pressure in the far-field.

Discussion and Conclusions
The radiated noise of a non-cavitating propeller is calculated by combining BEM and FWH acoustic analogy. Both the open water case (uniform inflow) and behind-hull case (non-uniform inflow) are investigated, and both the direct FWH and P-FWH (permeable FWH) approaches are used for the FWH analogy. The main conclusions are described in the following.
The SPL (Sound Pressure Level) obtained with the direct FWH approach and P-FWH approach with different permeable surfaces agrees well with other for BPF. However, for the 3rd BPF, obvious discrepancies are observed both between the results from the direct FWH and P-FWH approaches and between the results from P-FWH approach with different permeable surfaces. This indicates that the incompressibility assumption inside the permeable surface has larger influence at higher frequencies. This is also physically expected, because larger phase error would be produced for high frequencies than for low frequencies. A general conclusion about the modeling error level associated to P-FWH approach needs further detailed investigations.
The dominant FWH terms for the radiated noise and the decay rate of the noise in the propeller disk plane have been investigated. In the open water case, the thickness terms are dominant. The SP (Sound Pressure) has a very fast decay rate, that is, r −(1+Z) with Z being the blade number, in the relative near field, which is up to 100 m for the studied case. In the behind-hull case, the dominant FWH term are the first loading term, which contains the pressure's time derivative. The fast decay ends within two propeller diameters and the final r −1 decay rate is achieved around 40 m (five propeller diameters).
When using P-FWH with BEM for the behind-hull case, the fictitious volume flux on the permeable FWH surface is observed. A virtual source correction method is proposed to eliminate such a numerical error. The proposed method is simple and works very well. Also observed is that, with the P-FWH approach, the FWH term with normal velocity derivative is dominant.
The formula to consider the potential wake sheet as sound sources is proposed for the first time. The good correlation between the results obtained with the direct FWH and P-FWH approaches in Section 3.4 verifies the correctness of the proposed formula. Further comparison with high-quality viscous simulation to show the effectiveness of the proposed formula will be carried out in the near future. The numerical results reveal that the wake sheet induced noise depends on the length of the wake sheet. Normally in BEM, the length of the wake sheet is chosen according to the convergence of hydrodynamic forces on the blades. For the hydroacoustic simulation, a more physical consideration is necessary. As a recommendation, one needs to take consideration of the evolution and instability of propeller shedding vorticity system, as showed in experiments by Maroi Felli [29] and in numerical simulation by Roberto Muscari [30]. This would give a physical end position of the coherent turbulent structure, and the wake sheet should be modelled up to this position. The exact way to consider such a mechanism in BEM needs further investigations.