WEI Ying-san, WANG Yong-sheng, CHANG Shu-ping, FU Jian

Department of Mechanical Engineering, Naval University of Engineering, Wuhan 430033, China,

E-mail: weiyingsan@163.com

NUMERICAL PREDICTION OF PROPELLER EXCITED ACOUSTIC RESPONSE OF SUBMARINE STRUCTURE BASED ON CFD, FEM AND BEM*

WEI Ying-san, WANG Yong-sheng, CHANG Shu-ping, FU Jian

Department of Mechanical Engineering, Naval University of Engineering, Wuhan 430033, China,

E-mail: weiyingsan@163.com

A mesh-less Refined Integral Algorithm (RIA) of Boundary Element Method (BEM) is proposed to accurately solve the Helmholtz Integral Equation (HIE). The convergence behavior and the practicability of the method are validated. Computational Fluid Dynamics (CFD), Finite Element Method (FEM) and RIA are used to predict the propeller excited underwater noise of the submarine hull structure. Firstly the propeller and submarine's flows are independently validated, then the self propulsion of the“submarine+propeller” system is simulated via CFD and the balanced point of the system is determined as well as the self propulsion factors. Secondly, the transient response of the “submarine + propeller” system is analyzed at the balanced point, and the propeller thrust and torque excitations are calculated. Thirdly the thrust and the torque excitations of the propeller are loaded on the submarine, respectively, to calculate the acoustic response, and the sound power and the main peak frequencies are obtained. Results show that: (1) the thrust mainly excites the submarine axial mode and the high frequency area appears at the two conical-type ends, while the torque mainly excites the circumferential mode and the high frequency area appears at the broadside of the cylindrical section, but with rather smaller sound power and radiation efficiency than the former, (2) the main sound source appears at BPF and 2BPF and comes from the harmonic propeller excitations. So, the main attention should be paid on the thrust excitation control for the sound reduction of the propeller excited submarine structure.

submarine, propeller excitation, self propulsion, underwater noise, Computational Fluid Dynamics (CFD), Finite Element Method (FEM) and Refined Integral Algorithm (RIA)

Introduction

When the Boundary Element Method (BEM) is applied to calculate the near field sound pressure, with its singularity problem[1-3], the integral coefficients of Helmboltz Integral Equation (HIE) is not accurate, especially for a coarse Boundary Element (BE) mesh. Visser[4]proposed a local adaptive BEM quadrature scheme to deal with this problem, but the algorithm is complex and not very convenient to perform the mesh refinement. In this article, a global self-adaptive BEM quadrature named the Refined Integral Algorithm (RIA) is proposed to deal with the above problem.

Recently, the propeller and the flow excited submarine vibration and the noise radiation attracted much attention[5,6]. Erik[7]and Brouwer[8]analyzed the propeller induced submarine acoustic responses, but their focus was mainly on the submarine structure and the acoustic responses under a unit force of the propeller, not so much on the actual propeller excitation. To capture the fluctuated pressure of the flow on the underwater vehicle boundary, the panel method was often adopted. With the development of the Computational Fluid Dynamics (CFD), the CFD technique is gradually adopted to simulate the performance of the propeller and the submarine. Shen and Su[9]used the CFD method to simulate different unsteady hydrodynamic performances related with propellers, and the numerical results agree well with the experiment data. Wei and Wang[10]and Zhang et al.[11]calculated the fluctuated pressure of the submarine hull via CFD. In this article, the actual propeller excitation in the submarine viscous non-uniform wake is simulated via CFD, and the acoustic response of the submarine under the propeller fluctuated thrust and torque excitation is analyzed.

1. Refined integral algorithm for near-field acoustic response

1.1 Boundary element formulation

In a homogeneous medium, for the 3-D linear time-harmonic problem of the external acoustics with Neumann Boundary Condition (BC), the Helmholtz equation solution is where x is a general field point, y is the source point, u is the acoustic pressure, n(y) is the unit normal aty∈?D and directed into D, and D stands for the domain of the propagation, ?D is the boundary of D, α(x)is a geometry related coefficient, G(r )=-eikr/4π is the Green function in the free space, with. In order to obtain a unique solution of Eq.(1), N CHIEF points pointing in the normal direction of the interior boundary are introduced as shown in Fig.1, N is the number of BE.

Fig.1 Configuration of CHIEF points. “·” Source point, “○”CHIEF point

1.2 Numerical discretization

Divide ?D into N BEsγ1,γ2,…,γN, so Eq.(1) is discretized as

where vn(y) is the normal velocity. Eq.(2) is a set of 2N linear over-determined equations with N un-

knowns, and

By a variational operation, the over-determined system of Eq.(2) can be expressed as follows

For Eq.(4), the general minimal residual (GMRES) algorithm can be used to solve the over-determined equations, and once the sound pressure on the boundary is known, the pressure at any point in the exterior sound field can be determined by the following equation

In Eq.(2), when the source point coincides with a field point, aiiand biiwould become singular. As aiican be obtained with an analytical integration[12], so some advanced method like the adaptive Gauss quadrature is adopted to calculate bii.

1.3 Self-adaptive quadrature of quadratic quadrilateral boundary element

Visser[4]proposed a local mesh refinement sche me to solve the Helmholtz singular and nearsingular integrals. This scheme is effective, but complex and not very convenient in performing the mesh refinement. In this article, a global mesh refinement scheme is proposed to calculate bii. Adopting the quadrilateral (QUAD8) boundary element, by a coordinate transformation, the domain of the integral can be transformed to a standard square with the help of the shape function on each node. So the variables and their derivatives on QUAD8 can be interpolated via the shape functionNj, as

where (xj,yj)are the node coordinates,Nξj, Nηjthe partial derivatives of the shape function withrespect to ξ,η, respectively.

Generally, aij, bijand other coefficients can be written in a unified form as

Fig.2 Flowchart of the RIA

Here, an adaptive mesh refinement scheme is proposed to solve Eq.(6). Figure 2 shows the flowchart of this scheme, in which the initial BE (father element) is refined into four children elements, which can further be refined in the same manner. So, the integral on a coarse BE is turned to the integrals on these refined elements. At each stage of the mesh refinement, the sum of the integrals on all children elements of the same level serves as the approximate solution of that level. Generally speaking, the multirefinement of the BE is needed to get a convergent solution at a given allowable tolerance, and simultaneously, the temporary children elements in each level are just temporarily used, so the adaptive quadrature scheme is a kind of mesh-less techniques. With the help of Gauss quadrature, Eq.(5) can be further expressed as

where (ξi,ηj)are the Gaussian points, Hijis the weight,Ng=4Nris the number of Gaussian points, Nr is the refined level. Next, we are going to get the Gaussian points and the weight of the children elements at each level.

1.4 Multi-level boun dary element refinement algorithm

In the first level of the element refinement as shown in Fig.2, the initial father element P1,P2,P3,P4is divided into four children elements with five new nodes being born. The Gaussian points of each element can be determined by the nodes and calculated via the node transform matrices T1,T2,T3,T4, as

where X=[x1,x2,x3,x4]Tis the initial nodes. Once the children element nodes are determined, the Gaussian points (one point per element is chosen) can be determined via the Gaussian point transform matrices, as

where Xg0is the initial element Gaussian point, Xg1are thefirst level element Gaussian points. In order to determine the Gaussian points of the Nthlevel children elements, four operators are defined:“?” represents the matrix left-multiplication with thecell matrix, “?”represents the matrix right-multiplication with the cell matrix, “@”represents thecell matrix multiplication with the cell matrix, “^” represents the exponential multiplication of the matrix.

For example:

So, the node coordinates and the element Gaussian points of each level can be expressed as:

where Tr=[T1,T2,T3,T4], Xni,Xgiare the ithlevel element nodes and Gaussian point matrix,respectively. The weight of theNrthlevel children elementsw=S/4Nr, with Sbeing the initial square Nf f element area. So, once Xn0isgiven, XgNcan be determined and Eq.(7) canbe easily evaluated at each level.

Similarly, when the field point and the source point are very close to each other but not coincided, the RIA can also be applied to get the near-singular quadrature.

Fig.3 Integration on refined singular element

Fig.4Adaptive quadrature convergence behaviour whenkernel function: f(r)=1/r

1.5 Co nvergence analysis

Take the singular part of G(r) as an example, that: f(r)=1/r. Assume that the field point locates at the center of the initial square element as shown in Fig.3, and the integral on this BE is singular. As the element refined number N increases, the refined quadrature of Eq.(6) is converged at a constant value, see Fig.4. At each level, an approximate solution is compared to the previous level. The convergence behaviour of the integral solution in each refinement level is defined by ERRO, as

When the relative error ERRO is convergent at a given error ERROmax, usually, 10-4, no further refinement is necessary, also the maximum refinement number Nmaxis prescrib ed to prevent excessive temporary elements. It is shown that Nmaxshould be between 4and 8[4]. Figure 5 shows theconvergence behaviour and CPUrun time of RIA. So to ensure the accuracy and at the same time to keep the time-cost at somelevel, the optimum refinement number can be determined as Nmax=6.

Fig.5 Convergence behaviour and CPU run time of RIA to perform singular integration

Fig.6 Sound radiation pattern of the source

Fig.7 Sound dir ectivity

2. Numeri cal validation

In thissection, the performance of the RIA is investigated. The pulsating sphere is widely used as a validation case for the boundary element method, since the analytic solution is available. Firstly the BEM model of the sphere is made, then the RIA is used to calculate the HIE integral. Figure 6 shows the sound pressure in thenear-field and the far-field from 10 Hz to 1K Hz, Fig.7 shows the sound directiions, and a good agreement of the numerical solutions is obtained.

3. Application to the propeller excited acoustic responses of submarine underwater structure

It is shown that the propeller excitation canno-uniform wake of the submarine, the thrust and the torque of the propeller are fluctuated in the time do-the vibratory force transmitted through fluid is onlyforcedue to the fluid force is between 10-50% of that due to the shaft force. Therefore, in this article, only the propeller excitation transmitted through shaft is considered, as much stronger acoustic response is excited by it than by fluid force. In most of the studies of the propeller excited submarine noise, the propeller excitation is assumed to be constant or is substituted with a dipole[15,16]. However, when evaluating the acoustic performance of the propeller excited submarine noise, the actual excitation of the propeller should be known. Merz noted that the combination of CFD, FEM and BEM is the trend in the studies of the propeller induced submarine hull vibration and the underwater noise radiation in future[16].

In this section, CFD, FEM, RIA are respectively used to predict the submarine underwater vibration noise induced by the propeller excitation of the“SUBOFF submarine +4 381 propeller” system. Firstly the open water characteristics of the 4 381 propeller are calculated via CFD, and then the numerical results are compared with the experimental results to validate the credibility of the propeller results. Secondly the steady flow characteristics of the SUBOFF coefficients on the submarine hull as well as the sail and rudder section lines, and the wake velocity field of the submarine are compared with the experimental results to validate the credibility of the submarine results. Thirdly the self propulsion of the “submarine + propeller” system is simulated via CFD to obtain the self propulsion factors such as the advance ratio of the propeller, the thrust deduction factor, and the wake fraction, so in this way the balanced point of the system can be identified. Fourthly the transient flow of the “submarine + propeller” system is simulated viaCFD analysis to obtain the propeller fluctuated thrust and torque. Lastly the vibration and underwater noise radiation of the submarine excited by the propeller is calculated via FEM and RIA tools, respectively.

Fig.8 Mesh of the propeller single passage

Fig.9 Calculation domain of open water propeller

3.1 Validation of independent propeller’s flow

In this article, a five blade propeller model with di a meter D =0.25m is chosen to thrust the SUBOFF submarine. Firstly the single passage of the propeller is meshed with a structure mesh, with the local areas such as thetip and the root refined, as shown inFig.8. The total number of the single passage meshes is 380 000, they-plus value on the boundary of the blade and the root is about 84. Thecalculation domain of the open water propeller is shown in Fig.9, and the angular speed of the propeller is 650 rpm on thissection. Here, the advance ratio of the propeller changes from 0.1 to 1,so in this way, the thrust coefficient Kt, the torque coefficient Kq, and the open water efficiency η can be calculated via CFD.

Fig.10Open water characteristics of propeller

The numerical model is based on the RANS equations, with the body forces due to the rotation of the blade being treated based on the quasi-steady Multiple-Frame-of-Reference method. The turbulent flow within the blade is formulated in a rotating reference frame, the Shear Stress Transport (SST) turbulent model is adopted in this article. Figure 10 shows the open water characteristics of the propeller, and the Ktvalue agrees well with the experimental result, while the Kqvalue differs but a little with the experimental value.

Fig.11 Calculation domain of “submarine + propeller” system

Fig.12Mesh of the “submarine + propeller” system

3.2 Validation of independent submarine’s flow

In this section, the hydrodynamics of thesubmarine without propeller are simulated. The total length of the submarine is L=4.356m, the inflow velocity Vs=3.036m/s, corresponding to the Reynolds Number ReL=1.31×107. The fluid domain and the mesh of the single submarine is obtained from the“submarine +propeller” system as described in Fig.11 a nd Fig.12, where Dmaxis the max imum diameter of the submarine. Here, the SST turbulent model isadopted, and to capture the flow detail on the boundary of the submarine, the meshes near the boundary, the sail and the rudders of the submarine are refined. By a mesh independent analysis, when the total mesh number reaches2 060 000, the resistance of the submarine converges.

Table 1 Comparison of measured and calculated resistance values

Table 1 shows a comparison of the submarinevalidate the accuracy of the submarine flow, the local pressurecoeffi cients (Cp) are compared with experimental results[17]as well as the results in Ref.[17], as shown in Fig.13, and the wake velocity fieldofthe submarine in Fig.14 is also compared with the experimental results, which shows a good agreement.

Fig.13 Cpdistributions along the submarine local area

Fig.14 Wake velocity field of the submarine

3.3 Self propulsion of “submarine + propeller”

As the flowcharacteristics of the independent propeller and submarine are validated, so in this section, the propeller mesh and the submarine mesh are combined together to determine thebalanced point of the “submarine + propeller” system. The balanced point is calculated via a steady flow analysis by CFD, and the calculation domain and the mesh are shown in Fig.11 and Fig.12. The global number of nodes and elements is 4 129 092 and 3 970 770, and the y-plus value is about 80, and the SST turbulent model is adopted, with the advection term and the momentum equation discretized by the second order upwind scheme. The inflow velocity is set to 3.036 m/s. Changing the propeller rotation speed, the curve of the propeller thrust and the submarine resistance versus the rotation speed can be determined as shown in Fig.15. Then the balanced point is determined at the intersection of the two curves at the propeller rotation speed of 665 rpm, and in this case the thrust T, the resistance R, and the self propulsion factors such as the propeller advance ration J, the wake fraction w and the thrust deduction t can be determined as shown in Table 2.

Fig.15 Determination of the balanced point

Table 2 Self propulsion factor of the system

Fig.16 Fluctuation in time domain

3.4 Propeller excitation calculation

In this section, the fluctuated thrust and torqueof the propeller and the resistance of the submarine aresi mulated in a transient way at the balanced point as determined in the previous section. Firstly, the steady flow of the “submarine + propeller” system is simulated, and the result is used as the initial value of the transient analysis. In the transient analysis, the time step is set to 0.001 s, and thetotal simulation time is 10 s. The SST turbulent model is adopted in this case, the advection term and the momentum equation are discretized by the second order upwind scheme, the transient formulation is solved by the second order implicit scheme.

Fig.17 Propeller excitations in frequency domain

Fig.18 Submarine FE/BE model

Figure 16 shows the fluctuations of the thrust,the torque and the resistance in the time domain in a part of simulation time, which indicates a distinct periodic feature of the fluctuations. Figure 17 shows the propellerfluctuated thrust and torque in the frequencydomain. It can be seen that the main peak frequenciesof theexcitation spectrum are near the area where the propeller blade passes the frequencies of BPF and 2BPF, due to the interaction between the propeller and the non-uniform wake of the submarine, and the maximumexcitation in amplitude appears at BPF.

3.5 Acoustic response analysis

In this section, the underwater acoustic response of the submarine induced by propeller excitations are calculated by FEM and RIA, respectively. Firstly, the finite element model of the submarine is built as shown in Figs.18(a) and 18(b). Here, the submarine hull is divided into five cabins with four bulkheads and stiffened with ring ribs as well as longitude ribs. To ensure the local integrity of the hull, the sail and the stern are also stiffened. The hull and the bulkheads are represented by 26 739 SHELL63 elements, and the ribsby 5 517 BEAM188 elements. The thickness of the hull and the bulkheads is 0.006 m, the rib section on the hull is of inverted T shape, (3×20)/(4×12), the rib section on the appendage is of H shape, (2×20)/[2*(2×6)], and the submarine’s material is steel.

Then the propeller thrust and torque are, respectively, loaded on the submarine end bulkhead centre and the hull vibration responses are calculated. With the norm velocity of the hull elements as the boundary condition of HIE, and the finite elements of the hull are directly taken as the boundary elements as shown in Fig.18(c). In this way, the nodes and elements of the boundary element model are coincided with those of the finite elementmodel, and so no error of the data projection will be introduced.

Fig.19 Propeller excited submarine underwater sound power and equivalent sound power

Here, the RIA is adopted to calculate the submarinenear-field and far-field sound pressure. Toquantitatively evaluate the submarine underwater radiated noise, both the total sound power P and the Equivalent Radiated Power(ERP)are considered,that is[1]

where Re{} denotes the real part of thesummation,*denotes the conjugate complex, p is the boundary sound pressure, vnis the norm velocity. ERP refers to thevibration energy of the hull, while P re fers to the sound energy that the hull radiates, and the ratio ofPto ERP reflects the radiation efficiency of the submarine. Figure 19 shows the submarine sound power excited by thepropeller thrust and torque, Fig.20 shows the color map of the submarine boundary sound pressure obtained by using different excitation models. It is indicated that:

(1) Both P and ERP have peak frequencies at Blade Passing Frequency (BPF) and 2BPF, and their maximum values at BPF, and as the frequency increases from about 100 Hz, the sound power decreases.

This canbe due to the high thrust and torque excitation of the propeller at BPF and 2BPF, as shown in Fig.17, and at the high frequency, only the submarine’s local modes are excited, and the local modes contribute but little to the sound radiation, soP and ERP both decrease with the freque ncy at the high frequency range.

(2) The spectrum shape of P and ERP coincides with the thrust and torque spectrum shape.

This implies that the sound radiation pattern has a close relationship with the excitation pattern. So, it is advisable to control the propeller excitations to reduce the submarine’s underwater noise.

(3) Compared with the thrust, the torque excites the submarine’s radiated sound power level and the radiation efficiency in a much less extent.

Fig.20 Submarine boundary sound pressure at BPF in dB

This differences can be due to different hull vibration modes excited by the thrust and the torque. Asmainly vibrates in the axial direction and the normal response at the submarine’s two ends dominates the total hull’s vibration, and the main vibration and sound energy comes from the two ends of the submarine, as depicted in Fig.20(a). While the propellervertical direction, as shown in Fig.20(b). So, both the sound pressure level and the acoustic radiation efficiency excited by the propeller thrust are higher than those excited by the torque, and from that point of view, the torque excitation of the propeller can be neglected in the consideration of main excitations.

4. Conclusion

An accurate prediction of the propeller excitation is a key in the prediction of the propeller induced submarine underwater noise, and in this process, the self propulsion is the most important factor that affects the results. The submarine different vibration modes are excited by the propeller thrust and torque, as well as the sound modes, and the main sound source appears at BPF and 2BPF. The control of the propeller thrust excitation is advisable for submarine noise reductions.

[1]MARBURG S. Computational acoustics of noise propagation in fluids-finite and boundary element methods[M]. Berlin: Heidelberg Publication, 2008, 411-432.

[2]KOLM P., ROKHLIN V. Numerical quadratures for singular and hyper-singular integrals[J]. Computers and Mathematics with Applications, 2001, 41(1): 327-352.

[3]CARLEY M. Numerical quadrature for singular and hyper-singular integrals in boundary element methods[J]. SIAM Journal on Scientific Computing, 2007, 29(3): 1207-1216.

[4]VISSER R. A boundary element approach to acoustic radiation and source identification[D]. Ph. D. Thesis, Enschede, The Netherlands: University of Twente, 2004.

[5]CARESTA M., KESSISSOGLOU N. Low frequency structural and acoustic responses of a submarine hull[J]. Acoustics Australia, 2008, 36(2): 47-52.

[6]LV Shi-jin, YU Meng-sa and LI Dong sheng. Prediction of hydrodynamic radiation noise of underwater vehicle[J]. Journal of Hydrodynamics, Ser. A, 2007, 22(4): 475-482(in Chinese).

[7]ERIK V. Recent developments in predicting propeller induced hull pressure pulses[C]. Proceedings of the 1st International Ship Noise and Vibration Conference. London, UK, 2005, 1-8.

[8]BROUWER J. Ship propeller induced noise and vibrationsprediction and analysis[D]. Ph. D. Thesis, Enschede, The Netherlands: University of Twente, 2005.

[9]SHEN Hai-long, SU Yu-min. Study of propeller unsteady performance prediction in ship viscous non-uniform wake[J]. Chinese Journal of Hydrodynamics, 2009, 24(2): 232-241(in Chinese).

[10]WEI Ying-san, WANG Yong-sheng. Prediction of underwater vehicle vibration-noise based on singularitydecomposition and self-adaptive bem quadrature[J]. Chinese Journal of Computational Mechanics, 2011, 28(6): 858-863,925(in Chinese).

[11]ZHANG Nan, SHEN Hong-cui and YAO Hui-zhi. Large eddy simulation of wall pressure fluctuations of underwater vehicle[J]. Chinese Journal of Hydro- dynamics, 2010, 25(1): 106-112(in Chinese).

[12]WEI Ying-san, WANG Yong-sheng. A meshless selfadaptive boundary element quadrature scheme for underwater noise prediction[C]. Proceedings of 20th International Congress on Acoustics-ICA 2010. Sydney, Australia, 2010, 1-7.

[13]CARESTA M. Structural and acoustic response of a submerged vessel[D]. Ph. D Thesis, Sydney, Australia: University of New South Wales, 2009.

[14]CARESTA M., KESSISSOGLOU N. J. Acoustic signature of a submarine hull under harmonic excitation[J]. Applied Acoustics, 2010, 71(1): 17-31.

[15]KINNS R., THOMPSON I. Hull vibratory forces transmitted via the fluid and the shaft from a submarine propeller[J]. Ships and Offshore Structures, 2007, 2(2): 183-189.

[16]MERZ S., KESSISSOGLOU N. and KINNS R. excitation of a submarine hull by propeller forces[C]. 14th International Congress on Sound Vibration. Cairns, Australia, 2007, 1-8.

[17]ALIN N., FUREBY C. and SVENNBERG S. 3D unsteady computations for submarine-like bodies[C]. 43rd AIAA Aerospace Sciences Meeting and Exhibit. Reno, Nevada, 2005.

July 4, 2011, Revised November 18, 2011)

* Biography: WEI Ying-san (1984-), Male, Ph. D. Candidate Corresponding author: WANG Yong-sheng,

E-mail: yongshengwang666@126.com