Sound field properties of non-cavitating marine propeller

The sound field properties for non-cavitating marine propellers are studied using FW-H 1 (Ffowcs William-Hawkings) acoustic analogy and BEM (Boundary Element Method) approach. For 2 the first time, the FW-H formula to calculate acoustic pressure generated by the propeller wake sheet 3 is proposed. The corresponding sound signal can be dominant up to several kilometres. Secondly, 4 non-physical fictitious volume flux is observed when the permeable FW-H approach is used together 5 with BEM. The reason is explained and a virtual source correction method is proposed to solve this 6 problem. Furthermore, analytical analyse show that the second thickness term and the first two 7 loading terms in Farassat 1A formula are important, and the others are negligible. The numerical 8 studies show that the permeable FW-H approach produces underestimation when compared to the 9 direct FW-H method, and the underestimation is severer for larger permeable surface and higher 10 frequency. It is also found that the first loading term in Farassat 1A formula is the dominant source 11 term in the direct FW-H method, while the first thickness term is dominant when the permeable 12 FW-H approach is used. 13

The current work is based on the FW-H equation in the Farassat 1A form [7], 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 142 the thickness and loading terms, respectively; f = 0 denotes the integration surface for FW-H equation 143 (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 145 coordinate where the undisturbed fluid (or far-field fluid) is static; M is the vectorial Mach number, i.e. 146 M = v/c 0 ; r is the vector from the source point to observer position, r is its length, the subscript ret 147 means that variables inside the square brackets are computed at the corresponding retarded times, i.e.

148
at τ = t − r/c 0 ; the subscripted variables denote dot products of a vector and a unit vector implied by 149 the subscript, except for L M = L · M, and n is the unit normal vector on the surface directing into the 150 fluid domain; the dot over a variable means time derivative.

151
In equation (1) the quadrupole component due to the Lighthill stress tensor is ignored, because it vanishes in an incompressible potential fluid domain. As the observer is moving relative to the static fluid, the position change between the retarded time τ and the observer time t must be considered. This is done by calculating the vector r using the formula provided by Bres et al. [21], i.e. with and where ∆x, ∆y and ∆z are components of the vector from source point to the observer position at the 152 retarded time, and u 0 is the observer's velocity along the x direction. The motion in other directions 153 would not be considered in the current study. Normally, the global coordinate can be chosen to meet 154 this motion direction requirement.

155
To facilitate the detailed analysis of the components, the thickness and loading terms given in equations (2) and (3) are further decomposed as p T = p T,1 + p T,2 = p T,1 + (p T,2−1 + p T,2−2 ) (9) where Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 September 2020 doi:10.20944/preprints202009.0117.v1 On the solid surface, the impermeable condition implies u n = v n . Thus the term U n and L on the solid surface can be simplified to U n = v n (17) In BEM, the wake sheet is used to represent the shedding vorticity of the blades. The wake sheet 157 has zero thickness and acts as a potential boundary in the solution process, as shown in Figure 1. In 158 this section, we will show that the zero-thickness surface representation of the vorticity field turns the 159 volume integration of the Lighthill stress term into a surface integration on the wake sheet.

160
The evolution of the wake sheet can be considered in two perspectives, or two coordinate reference 161 systems, i.e. the earth-fixed global reference system and the propeller-fixed local reference system. To

178
We now discrete the upper and lower subsurfaces in the same manner with the wake sheet. Then 179 every FW-H term will be checked on the panel. It will be shown that thickness terms cancel out each 180 other, while the loading terms do not totally vanish and form the FW-H formula on the wake sheet. 181 We use the superscript '+' to denote the variable values on the upper subsurface and '-' for those on 182 the lower subsurface, no superscript means the summed effect. Because there is no potential source on the wake sheet, so the normal component of the fluid velocity is equal on the upper and lower subsurfaces. Due to the opposite normal directions, the following relationship holds: Because the upper and lower subsurfaces share the same position with the wake sheet, the surface moving velocities are the same, i.e. v + = v − .
With these two relationships, the following results can be deduced These imply further that on the wake sheet p T = p T,1 = p T,2 = 0 (21)

185
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, i.e.
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 the unsteady Bernoulli equation, it can be obtained that Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 September 2020 doi:10.20944/preprints202009.0117.v1 The difference (u + ind − u − ind ) is related to the vortex density on the surface. It is zero on the wake panels in a lower order BEM code. Considering the relationship between doublet strength and potential jump on the wake panel, i.e. −µ = (φ + − φ − ), we have For the second term on the right-hand side of equation (22), the difference (u + − u − ) is related to 186 the vortex density on the surface, which is zero on the wake panels in a lower order BEM code. Thus, 187 the second term would vanish.

188
Combination of these two terms leads to This can be substituted into equations (14)  case, the term does not vanish and will generate acoustic pressure.

191
The numerical results will show that the term p L,1 is dominant among the wake sheet induced acoustic pressure. To enable a further detailed analysis, the term p L,1 is separated further to p L,1−1 and p L,1−2 , i.e.
In the earth-fixed perspective, because of non-time-varying doublet strength on the wake panels,  The second thickness term p T,2 is also influenced. The first thickness term is an exception. If we choose 211 the permeable surface so that the normal vectors on the surface do not change with time, then Uṅ will 212 not be influenced. If the ship wake velocity is timely constant, which is a canonical assumption in 213 the potential solver, thenU n is also not influenced. Thus, under the two above mentioned conditions, 214 the term p T,1 keeps the same when the ship wake velocity is included or ignored on the permeable permeable formula is used. Thus it becomes unimportant whether to include the ship wake velocity in 217 u. 218 The second problem is caused by the use of 2D ship wake velocity information as a 3D velocity 219 field. This results in a non-divergence-free velocity field. The potential source on the blade surface, 220 which is calculated according to this velocity field, would therefore not integrate to zero. This will  In this section, we conduct a rough mathematical analysis for the orders of magnitude of the FW-H terms on the blade surfaces. Without loss of generality, the propeller's rotation and forward speed are assumed to be constant, and the propeller is regarded as a rigid body. Under these assumptions, the normal surface motion velocity v n do not change with time, i.e. dv n dt ≡ 0 (29) Before the formal analysis, two concepts are introduced, i.e. the far-field term and the near-field 231 term, which have the decay rates of r −1 and r −2 , respectively. It is obvious that with enough distance, 232 i.e. in the far-field, the far-field terms are dominant over the near-field terms.  For the first thickness term p T,1 , we havė Furthermore, the following proportional relations exist where D is the propeller diameter, v is the magnitude of v, and ω is the rotation speed. Together with the condition M r << 1, and considering the integrand of p T,2−1 in equation (12), i.e.
it can be obtained that With the same procedure, the proportional relation for p T,2−2 can be expressed as Then, for the loading terms, the following relations are used where θ is the azimuthal angle. Using these relations we obtain With these proportional relations in equations (31) to (35), some insights into the relative 239 significance of the FW-H terms can be achieved. Firstly, considering the situation where ωD/c 0 ∝ 240 M << 1, it can be concluded that p L,3 is always secondary when compared to p L,1 + p L,2 for the 241 marine propeller problem. Then, only four terms remain to be investigated, which are p T,2−1 , p T,2−2 , 242 p L,1 , and p L,2 .

243
It can be observed that p T,2−1 and p L,1 are far-field terms and p T,2−2 and p L,2 are near-field terms. Besides, p T,2−1 and p L,1 are proportional to the same kernels, and the terms p T,2−2 and p L,2 also. If p T,2−1 and p L,1 are compared using their original integrands in equations (12) and (14), it turns out to be the comparison ρ 0 v 2 vs. p + dp/dθ.
For the near-field terms p T,2−2 and p L,2 it would be the comparison Taking the KP505 propeller as an example, the surface motion velocity v is around 30m/s near the 244 tip, and ρ 0 v 2 approaches the order of magnitude 10 6 . For a ship propeller, p is in the same order of 245 magnitude with the atmospheric pressure, i.e. 1 × 10 5~2 × 10 5 . For vehicles working in deep water, 246 p would have a different order of magnitude. The term dp/dθ is quite small for open water case 247 and smaller than p when ship wake field is smooth, but can also achieve the order of 10 5 when the 248 propeller works in a severe ship wake field. Thus, the thickness terms p T,2−1 and p T,2−2 seem to be 249 more significant than the loading terms p L,1 and p L,2 , with integrands having one order higher of magnitude. The numerical results in section 4.1 will show that it is really the case for the open water 251 case.

252
However, in unsteady cases, the motion velocity has the same value on different blades, while 253 p and dp/dθ are highly θ-dependent. The interference then leads to a much faster decay rate of the 254 thickness terms than that of the loading terms. This result in the dominance of the loading terms over 255 the thickness terms in unsteady cases.

256
When comparison are conducted between the near-field and far-field terms, i.e. between p T,2−1 in 257 (31) and p T,2−2 in (32) or between p L,1 in (33) and p L,2 in (34), it can be found that the distance where 258 the dominant term changes is related to c 0 /ω = λ 1 /(2πZ), where λ 1 is the acoustic wavelength for 259 the 1st BPF. However, this formula can not directly be used to calculate the critical distance. Other 260 factors, e.g. the interference, must be considered.

262
To study the decaying behaviour and determine the dominant FW-H term at different distances, a 263 set of observer points are placed in the propeller disk plane. They are located on the right side of the 264 propeller and the distances to the propeller centre are from 5m to 550m, which is around 70D or 2.3λ 1 .

265
The numerical result will be presented mainly for the 1st BPF.

279
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 280 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 281 rate r −(1+Z) in the near-field and r −2 in the far-field. That means the decay rate is increased in the and r −1 in the far-field for both the total thickness term and total loading term. This is the combined 288 effect of the behaviours of all the terms shown in Figure 5 here. The sound pressure level obtained with different propeller wake sheet lengths are given in Figure   294 7, in the form of separate contributions from blade surface, wake sheet and combined results. There the wake sheet has a close relation to the finite sound speed, which is not captured in the Bernoulli 300 pressure signal.

301
The wake sheet induced SPL increases with the wake sheet length. For short wake sheet lengths, i.e.

302
1 and 2 revolutions, the blade induced SPL is dominant. For long wakes, i.e. no less than 3 revolutions, 303 the wake sheet induced SPL is dominant in a large distance rage (< 400m for 3 revolutions). The wake 304 sheet induced SP has a decay rate of r −2 , which is faster than that of the blade induced SP. This implies 305 that the blade would always be the dominant role in the very far field, but the distances might need 306 to be larger than several kilometres. Figure 8 shows the variation of the SPL up to 4km for the cases 307 with 2 and 3 wake revolutions. With 3 wake revolutions, the wake sheet induced SPL crosses the blade 308 induced SPL around 700m and the wake sheet effect can be ignored firstly at the distance around 3km.   Unlike in the open water case, in the unsteady case p L,1 and p L,2 achieve the final decay rates of r −1 314 and r −2 quite early, i.e. at the distance around 16m (2D). The reason is that the pressure and also the 315 time derivative of pressure are not equal on each blade, which is caused by the inhomogeneous ship 316 wake field. Therefore the interference effect becomes quite weak.  shown. Beyond 40m (5D), the dominant term is p L,1 , more exactly p L,1−2 , which is defined in equation 320 (28). Beyond 5D the decay rate of p L,1−2 is firstly r −2 and then changes to r −1 at a position, which is 321 dependent on the wake sheet length. The term p L,1−1 decays always with the rate r −1 beyond 5D. In 322 the both cases, p L,1−2 achieves the decay rate of r −1 before it becomes smaller than p L,1−1 . Thus the 323 dominant term beyond 5D is always p L,1−2 .  In this section, the acoustic pressure is calculated with the FW-H permeable surface. The wake 337 sheet length is set to 3 revolutions for the study in this section.

338
Without using the virtual source correction, we have firstly studied the sensitivity of the acoustic 339 pressure to the permeable surface size. Four different cylindrical permeable surfaces have been studied.

340
They are denoted with the names basic, small, large, and long, respectively. Their size and discretization 341 parameters are listed in table 2 and depicted in Figure 11. The adopted panel sizes are chosen through 342 mesh independence study for every permeable surface so that the mesh convergence is achieved.

343
Generally, when the permeable surface is small, it is located near to the potential boundaries (blades and wake sheets), so the mesh on it should be fine to capture the flow field variation correctly. When 345 the permeable surface is large, the mesh can be relatively rough due to the smooth flow field.  Then we come back to the premature arrival of the r −1 decay rate shown in Figure 12. According to the analysis in section 2.4, due to the use of 2D ship wake velocity information as a 3D velocity field in BEM, the potential source on the propeller surface do not integrate to zero. This leads to a non-physical fictitious volume flux on the permeable surface. The sound pressure amplitude induced by the volume flux (like a monopole) decays as r −1 . When it dominates the acoustic pressure field, the total SP would have a decay rate of r −1 . Figure 14 shows the timely variation of the integral potential source on propeller surface and the volume flux on the permeable surface, i.e. large permeable surface they are 2e-3m 3 and 2e-4m 3 , respectively. For the basic one they are 0.015m 3 369 and 0.02m 3 , respectively. Generally, larger permeable surface and finer discretization lead to smaller 370 volume flux. Figure 15 shows the SPL obtained with and without the virtual source correction, using 371 the basic permeable surface. An improvement for almost the whole distance range can be observed.  With the virtual source correction, the influence of the permeable surface size is restudied for both 373 the 1st and 3rd BPF. The results are given in Figure 16. When compared to the direct FW-H results, a 374 general underestimation can be observed. Larger permeable surfaces have more underestimation. The 375 underestimation is more obvious for the higher frequency. This is because the wavelength is shorter, 376 and then the phase difference due to the incompressible assumption inside the permeable surface

382
The property of the sound field calculated with FW-H equation for non-cavitating propeller as 383 well as the property of different FW-H terms are studied with the potential BEM (Boundary Element 384 Method) approach, taking its advantage of zero numerical diffusion and the lack of spatial vorticity 385 distribution, which could be hard to be handled correctly.  In the unsteady case, the dominant terms on the blade surface and the wake sheet are both the 398 first loading term in Farassat 1A formula, which contains the pressure's time derivative. They are both 399 related to the azimuthal derivative of the ship wake velocity. Thus, a precise description of the ship 400 wake field regarding the azimuthal continuity and variation is a prerequisite for the correct acoustic 401 pressure prediction.

402
The way to use ship wake field in BEM, i.e. using the velocities on a 2D plane as a 3D distribution, 403 produces a non-divergence-free velocity field. This leads to a non-vanishing integral source strength on the blade surface, which produces a fictitious volume flux when the permeable FW-H approach is 405 used. A virtual source correction method is proposed to cancel out such a numerical volume flux. The

406
proposed method is quite simple and proved to work very well.

407
With the permeable FW-H approach, the thickness term is dominant. The term with the time 408 derivative of normal velocity is especially important. It is also found that the permeable FW-H always 409 produces an underestimation, which is around 3dB for the BPF (Blade Pass Frequency) and becomes 410 severer for higher frequencies. A larger permeable surface leads to further underestimation. This 411 might be due to the ignoring of sound travel time inside the permeable surface.

412
It needs to be emphasized that the formula to consider the acoustic contribution of the propeller team are also acknowledged for providing the easy-to-use BEM code, which make the current study 420 possible.