片空化狀態(tài)下對(duì)轉(zhuǎn)螺旋槳噪聲特性仿真分析

王順杰，王易川，戴衛(wèi)國，李斌

片空化狀態(tài)下對(duì)轉(zhuǎn)螺旋槳噪聲特性仿真分析

王順杰，王易川，戴衛(wèi)國，李斌

（海軍潛艇學(xué)院，山東青島266044）

為了研究對(duì)轉(zhuǎn)螺旋槳空化噪聲特性，利用Schnerr-sauer空化模型及RNG k-ε湍流模型，并采用動(dòng)網(wǎng)格模型對(duì)對(duì)轉(zhuǎn)螺旋槳進(jìn)行空化數(shù)值分析。首先，對(duì)DTMB4381螺旋槳進(jìn)行空化敞水性能計(jì)算，所得空泡形態(tài)與公開發(fā)表試驗(yàn)及數(shù)值模擬結(jié)果吻合較好。然后對(duì)對(duì)轉(zhuǎn)螺旋槳空化噪聲進(jìn)行計(jì)算及分析，求解了不同監(jiān)測點(diǎn)處的聲壓和聲壓功率譜密度。從計(jì)算結(jié)果可以看出噪聲功率譜中存在著與前、后槳葉頻相關(guān)的線譜成份。

片空化；對(duì)轉(zhuǎn)螺旋槳；數(shù)值分析

1 引言

螺旋槳空化噪聲是船舶輻射噪聲的主要因素，常見的螺旋槳空泡主要有梢渦空泡、片空泡、云空泡、轂渦空泡等，由于空化物理問題的復(fù)雜性，多采用實(shí)驗(yàn)方法對(duì)空化現(xiàn)象進(jìn)行研究[1-2]。隨著CFD技術(shù)的發(fā)展，使用數(shù)值方法研究螺旋槳空化問題逐漸成為國際水動(dòng)力學(xué)界的前沿和熱點(diǎn)問題。而且對(duì)于聲學(xué)問題的直接數(shù)值模擬，可以在任意場點(diǎn)上對(duì)噪聲進(jìn)行監(jiān)測，有效避免了試驗(yàn)中水聽器安裝及其與流場之間的相互干擾問題。因此，數(shù)值計(jì)算將成為解決實(shí)際聲學(xué)問題的一種重要方法。

劉登成等[3]利用RANS求解其采用非結(jié)構(gòu)化網(wǎng)格對(duì)螺旋槳敞水性能和空泡性能進(jìn)行了數(shù)值計(jì)算。楊瓊方等[4]采用改進(jìn)的Sauer空化模型和修正的SST湍流模型，分析研究了非均勻進(jìn)流對(duì)螺旋槳空化水動(dòng)力性能的影響。胡健等[5]用數(shù)值方法計(jì)算了螺旋槳的擾動(dòng)速度，研究了螺旋槳尾流場的分布規(guī)律。Lindau等[6]用Kunz模型，和用UNCLE-M對(duì)螺旋槳P4831的空泡特性進(jìn)行了分析，在較大范圍的進(jìn)速系數(shù)下，均預(yù)報(bào)出使推力系數(shù)和扭矩系數(shù)急劇減小的臨界空泡數(shù)，并與試驗(yàn)結(jié)果吻合較好。

水下高速航行體普遍采用對(duì)轉(zhuǎn)螺旋槳的推進(jìn)方式。當(dāng)其航速很高或未在大深度航行時(shí)，螺旋槳的局部空化常常不可避免。因此，研究對(duì)轉(zhuǎn)螺旋槳空化噪聲特性對(duì)于目標(biāo)判別具有重要意義。目前，采用CFD技術(shù)對(duì)螺旋槳空化進(jìn)行的研究主要集中在單槳方面，未見公開發(fā)表的利用CFD技術(shù)對(duì)對(duì)轉(zhuǎn)螺旋槳空化噪聲特性進(jìn)行分析的相關(guān)文獻(xiàn)。本文采用Schneer-sauer空泡模型及RNG k-ε湍流模型，利用動(dòng)網(wǎng)格模型技術(shù)對(duì)對(duì)轉(zhuǎn)螺旋槳進(jìn)行空化噪聲特性進(jìn)行了數(shù)值分析。

2 數(shù)學(xué)模型

2.1 控制方程

假定氣、液兩相間無速度滑移，則氣、液兩相的連續(xù)方程和動(dòng)量方程分別為：

式中：m表示混合物，ρm為混合物密度，p為壓力，ui為直角坐標(biāo)系中絕對(duì)速度矢量的i分量，μm為混合物粘度，μt為湍流粘度。

2.2 Schnerr-sauer空化模型

采用Schnerr-sauer空化模型處理螺旋槳空化過程，并考慮相間滑移和水中未溶解氣體對(duì)空化的影響，空化流動(dòng)的連續(xù)性方程為：

水蒸氣空泡相：

式中：ρv為水蒸氣密度，a為水蒸氣體積分?jǐn)?shù)，Re和Rc分別為考慮蒸發(fā)和凝結(jié)的源項(xiàng)。根據(jù)Rayleigh-Plesset方程描述空泡的增長和潰滅過程，不同空化模型源項(xiàng)有不同的表達(dá)形式。

Schnerr-sauer空化模型為：

式中：ρl為液體密度，RB為空泡半徑，pv為飽和蒸汽壓。

2.3 RNG k-ε湍流模型

RNG k-ε湍流模型是基于重整化群的理論提出的，湍動(dòng)能k和湍動(dòng)耗散率ε的運(yùn)輸方程如下：

式中：湍流粘度μt=ρCuk2/ε，模型常數(shù)C2=1.68，Cu=0.085，σk=σε=0.717 9。β=0.012，η=Sk/ε，

與標(biāo)準(zhǔn)k-ε模型相比，RNG k-ε湍流模型通過修正湍流粘度，考慮了平均流動(dòng)中的旋轉(zhuǎn)及旋流流動(dòng)情況，適用于強(qiáng)流線彎曲、漩渦和旋轉(zhuǎn)流。

2.4 噪聲模型

FLUENT中采用Ffowcs Williams和Hawkings提出的FW-H方程模擬聲音的產(chǎn)生與傳播，此方程中采用了Lighthill的聲學(xué)近似模型。噪聲模塊的噪聲處理主要采用FW-H式，該式由連續(xù)性方程與Navier-Stokes方程導(dǎo)出，具體形式如下：

式中：ui為xi方向流體速度分量，un為垂直于面f=0的流體速度分量，vi為xi方向面運(yùn)動(dòng)速度分量，vn為垂直于面f=0的面運(yùn)動(dòng)速度分量，δ（f）為delta函數(shù)，H（f）為Heaviside函數(shù)，P′為遠(yuǎn)場聲壓（P′=P-P0），面f=0為源面，可為物體表面或具有滲透性的非物體表面，ni為指向外域（f＞0）的垂直單位矢量，c0為遠(yuǎn)場聲速，Tij為Lighthill張力，定義為：

（9）式的推導(dǎo)利用了格林函數(shù)δ（g）/4πr，方程右邊三項(xiàng)分別表示單極子、偶極子和四極子。該方程的解包括兩個(gè)面積分和一個(gè)體積分，其中面積分反映了單極子、偶極子和部分四極子的貢獻(xiàn)，而體積分反映了控制面外四極子的貢獻(xiàn)。在水中聲速比螺旋槳旋轉(zhuǎn)速度大很多，因此四極子項(xiàng)可以被忽略，F(xiàn)LUENT中既如此，因此壓力場求解式為：

3 數(shù)值方法驗(yàn)證

近年來，由泰勒水池（DTMB）提供的DTMB4381螺旋槳模型被大量研究者引用，作為檢驗(yàn)其螺旋槳敞水性能及空化數(shù)值仿真方法的計(jì)算對(duì)象。據(jù)此，本文同樣采用該螺旋槳模型作為檢驗(yàn)數(shù)值仿真方法可靠性及計(jì)算精度的標(biāo)準(zhǔn)，同時(shí)，可以很方便地參考已公開發(fā)表的DTMB4381螺旋槳相關(guān)的文獻(xiàn)制定計(jì)算策略。表1為DTMB4381螺旋槳模型的主要參數(shù)。

表1 DTMB4381槳模主要參數(shù)Tab.1 Main parameters of DTMB4381 propeller model

3.1 計(jì)算域及網(wǎng)格劃分

計(jì)算區(qū)域設(shè)置來流段取3D（D為螺旋槳直徑），去流段取7D，徑向直徑為4.2D。分為螺旋槳所在的旋轉(zhuǎn)區(qū)域及其以外的固定區(qū)域兩部分。使用GAMBIT軟件對(duì)計(jì)算域進(jìn)行網(wǎng)格劃分，其中螺旋槳所在旋轉(zhuǎn)區(qū)域采用非結(jié)構(gòu)網(wǎng)格，固定區(qū)域則使用結(jié)構(gòu)網(wǎng)格。將GAMBIT輸出的網(wǎng)格文件導(dǎo)入到TGID中，對(duì)螺旋槳表面附近劃分邊界層網(wǎng)格，邊界層網(wǎng)格類型為三棱柱。

圖1 槳葉表面網(wǎng)格布置Fig.1 Mesh division of blade surface

3.2 邊界條件

進(jìn)口邊界設(shè)為速度進(jìn)口邊界條件，出口邊界設(shè)為壓力邊界條件。旋轉(zhuǎn)區(qū)域與固定區(qū)域之間通過設(shè)置交界面完成兩個(gè)區(qū)域之間的數(shù)據(jù)傳遞。采用FLUENT中提供的動(dòng)網(wǎng)格模型控制旋轉(zhuǎn)區(qū)域的網(wǎng)格按照給定螺旋槳轉(zhuǎn)速整體運(yùn)動(dòng)，這種旋轉(zhuǎn)區(qū)域整體運(yùn)動(dòng)的處理方式可以避免因網(wǎng)格重構(gòu)造成的運(yùn)動(dòng)壁面附近較大的計(jì)算誤差，同時(shí)也無需考慮因網(wǎng)格重構(gòu)容易產(chǎn)生負(fù)體積對(duì)計(jì)算時(shí)間步長的限制。槳轂和槳葉設(shè)為壁面。

3.3 計(jì)算結(jié)果

依據(jù)國際船模水池會(huì)議的CFD不確定度分析臨時(shí)規(guī)程[7]中收斂性研究，為進(jìn)行網(wǎng)格收斂研究至少需要三套網(wǎng)格，并且各網(wǎng)格之間須滿足一定的細(xì)化率要求。本文三種網(wǎng)格數(shù)量均是按照槳葉表面網(wǎng)格以0.7左右的比例進(jìn)行細(xì)化，這組網(wǎng)格分別為Mesh1(1 488 654個(gè)網(wǎng)格)、Mesh2(2 236 952個(gè)網(wǎng)格)、Mesh3(3 240 491個(gè)網(wǎng)格)。

表2為取進(jìn)速V=3.048 m/s，頻率n=14.28 rps時(shí)對(duì)非定常條件下螺旋槳水動(dòng)力性能進(jìn)行模擬的推力系數(shù)KT、轉(zhuǎn)矩系數(shù)KQ的計(jì)算結(jié)果與試驗(yàn)值[8]比較結(jié)果。定義收斂因子：

其中：εT32為Mesh3與Mesh2推力系數(shù)解之差，εT21為Mesh2與Mesh1推力系數(shù)解之差，εQ32為Mesh3與Mesh2轉(zhuǎn)矩系數(shù)解之差，εQ21為Mesh2與Mesh1轉(zhuǎn)矩系數(shù)解之差。

表2 推力系數(shù)KT和轉(zhuǎn)矩系數(shù)KQ結(jié)果Tab.2 Calculation result of push coefficient KTand torque coeffiecient KQ

計(jì)算可得推力系數(shù)及轉(zhuǎn)矩系數(shù)的收斂因子RKT=0.33、RKQ=0.72，由于兩收斂因子均小于1，因此三套網(wǎng)格是單調(diào)收斂的，網(wǎng)格收斂性滿足，即網(wǎng)格數(shù)量增加到一定時(shí)，RKT與RKQ數(shù)值解開始收斂，網(wǎng)格數(shù)再增加對(duì)計(jì)算結(jié)果影響較小。

圖2為三種網(wǎng)格在進(jìn)速系數(shù)J=0.7，空泡數(shù)σ=3.5條件下的螺旋槳葉面片空化大小的數(shù)值計(jì)算結(jié)果。圖3（a）、（b）分別是Boswell[8]利用實(shí)驗(yàn)觀察及Kim[9]數(shù)值計(jì)算得到的DTMB4381螺旋槳在該條件下槳葉表面空化情況，圖3（a）中黑色區(qū)域代表發(fā)生空化的位置，圖3（b）中槳葉白色區(qū)域代表發(fā)生空化的位置。通過比較可以看出，本文采用的三種網(wǎng)格對(duì)片空化位置的預(yù)報(bào)基本一致，但對(duì)片空化程度的預(yù)報(bào)不一致；網(wǎng)格Mesh1對(duì)片空化程度的預(yù)報(bào)與實(shí)驗(yàn)結(jié)果相差最大，Mesh2與Mesh3對(duì)片空化程度的預(yù)報(bào)與文獻(xiàn)中的結(jié)果較為接近。

圖2 螺旋槳空化數(shù)值計(jì)算結(jié)果Fig.2 Result of propeller cavitation numerical calculation

圖3 Boswell試驗(yàn)及Kim計(jì)算得到的DTMB4381螺旋槳空泡圖Fig.3 Pictures of DTMB4381 propeller cavity from Boswell’s experiment and Kim’s calculation

上述計(jì)算結(jié)果表明，采用RNG k-ε湍流模型配合Schnerr-sauer空化模型在合適的網(wǎng)格尺度下能準(zhǔn)確地預(yù)報(bào)螺旋槳水動(dòng)力性能及螺旋槳空化。雖然三種網(wǎng)格尺度均能較好地預(yù)報(bào)螺旋槳的水動(dòng)力性能，但考慮到Mesh3預(yù)報(bào)的的水動(dòng)力性能誤差最小且所用的網(wǎng)格尺度能很好地模擬空化區(qū)域，因而后文中對(duì)對(duì)轉(zhuǎn)螺旋槳片空化狀態(tài)下的噪聲特性仿真分析以Mesh3的網(wǎng)格尺度為標(biāo)準(zhǔn)劃分網(wǎng)格。

4 對(duì)轉(zhuǎn)螺旋槳片空化狀態(tài)下噪聲特性分析

本文所采用對(duì)轉(zhuǎn)螺旋槳的主要參數(shù)如表3所示，其模型如圖4所示。考慮對(duì)轉(zhuǎn)槳的主要用途，在對(duì)轉(zhuǎn)槳前方布置一航行體。設(shè)置來流段為2倍航行體長度，去流段為5倍航行體長度，圓柱半徑為12倍航行體半徑。對(duì)轉(zhuǎn)槳區(qū)域的網(wǎng)格劃分按第2節(jié)中Mesh3的網(wǎng)格大小進(jìn)行，整個(gè)計(jì)算域網(wǎng)格總量為5 143 322。

前、后螺旋槳的轉(zhuǎn)速均為n=33 rps，來流速度V=43 kns，通過調(diào)節(jié)出口壓力控制空化數(shù)，使對(duì)轉(zhuǎn)槳發(fā)生片空化。圖6、圖7即為前、后槳發(fā)生空化的區(qū)域。

本文利用FLUENT14.0中的聲學(xué)計(jì)算模塊Acoustics-Ffowcs-Williams&Hawkings進(jìn)行聲場分析，F(xiàn)W-H方程由連續(xù)方程與Navier-Stokes方程導(dǎo)出。首先存儲(chǔ)足夠多的時(shí)間步長上的測量點(diǎn)的壓力脈動(dòng)值，作為聲場計(jì)算的聲源項(xiàng)輸入，可得到該點(diǎn)的聲壓值隨時(shí)間的變化。對(duì)此聲壓時(shí)間序列做FFT變換，即可得到功率譜密度在頻率上的變化情況。同時(shí)，為便于分析對(duì)轉(zhuǎn)螺旋槳空化噪聲的聲場分布等特性，在噪聲模塊中設(shè)置了6個(gè)聲壓監(jiān)測點(diǎn)，監(jiān)測點(diǎn)的位置如圖8所示。

圖4 對(duì)轉(zhuǎn)螺旋槳模型Fig.4 Counter-rotation propeller model

圖5 槳葉表面網(wǎng)格布置Fig.5 Mesh division of blade surface

圖6 前槳空化區(qū)域示意圖Fig.6 Sketch map of front propeller cavitation area

圖7 后槳空化區(qū)域示意圖Fig.7 Sketch map of rear propeller cavitation area

4.1 噪聲聲壓脈動(dòng)特性

圖9至圖12為1、2、4、5號(hào)監(jiān)測點(diǎn)的聲壓脈動(dòng)圖。從圖中可以看出，監(jiān)測點(diǎn)1、2具有相同的聲壓變化規(guī)律，其幅值隨著距離的增加不斷減小；監(jiān)測點(diǎn)4、5具有相同的聲壓變化規(guī)律，其幅值隨著距離的增加不斷減小。

4.2 噪聲聲壓功率譜特性分析

前、后槳轉(zhuǎn)速均為33 rps，因此前槳葉頻F1=198 Hz，后槳葉頻F2=231 Hz。分別觀察6個(gè)點(diǎn)在0-1 200 Hz頻率范圍內(nèi)的功率譜特性。其中監(jiān)測點(diǎn)1-3均具有圖13（a）、（b）所示的1號(hào)監(jiān)測點(diǎn)的功率譜密度分布。從圖13（a）中可以看出396 Hz=2F1和462 Hz=2F2頻率對(duì)應(yīng)的功率譜密度最大，另外從圖13（b）中可以看出在198 Hz=F1，231 Hz=F2，660 Hz=F1+2F2，825 Hz=3F1+F2，792 Hz=4F2，1 056 Hz=3F1+2F2，1 089 Hz=2F1+3F2處存在7條較為明顯的與前、后槳葉頻相關(guān)線譜。監(jiān)測點(diǎn)4-6均具有圖14（a）、（b）所示的4號(hào)監(jiān)測點(diǎn)的功率譜密度分布。從圖中可以看出在198 Hz=F1、231 Hz=F2、429 Hz=F1+F2、594 Hz=3F1、627 Hz=2F1+F2、858 Hz=2F1+2F2、891 Hz= F1+3F2、1 122 Hz=F1+4F2處存在8條較為明顯的與前、后槳葉頻相關(guān)線譜。因此，在對(duì)轉(zhuǎn)螺旋槳空化噪聲功率譜中存在著前槳葉頻線譜成份、后槳葉頻線譜成份以及前、后螺旋槳共同作用下的線譜成份。

圖8 監(jiān)測點(diǎn)布置圖Fig.8 Picture of monitor layout

圖9 監(jiān)測點(diǎn)1的聲壓脈動(dòng)圖Fig.9 Acoustic pressure pulse at monitor 1

圖10 監(jiān)測點(diǎn)2的聲壓脈動(dòng)圖Fig.10 Acoustic pressure pulse at monitor 2

圖11 監(jiān)測點(diǎn)4的聲壓脈動(dòng)圖Fig.11 Acoustic pressure pulse at monitor 4

圖12 監(jiān)測點(diǎn)5的聲壓脈動(dòng)圖Fig.12 Acoustic pressure pulse at monitor 5

圖13 1號(hào)監(jiān)測點(diǎn)功率譜密度Fig.13 Power spectral density at monitor 1

圖14 4號(hào)監(jiān)測點(diǎn)功率譜密度Fig.14 Power spectral density at monitor 4

5 結(jié)論

本文利用計(jì)算流體力學(xué)方法對(duì)對(duì)轉(zhuǎn)螺旋槳空化噪聲進(jìn)行了數(shù)值計(jì)算。計(jì)算時(shí)，采用的Schneersauer空泡模型、RNG k-ε湍流模型及動(dòng)網(wǎng)格技術(shù)模擬螺旋槳空化，通過本文研究，可得以下結(jié)論：

（1）Schnerr-sauer空化模型及RNG k-ε湍流模型的組合方式可以較好地模擬螺旋槳槳葉表面空化問題。

（2）在對(duì)轉(zhuǎn)螺旋槳空化噪聲功率譜中存在著前槳葉頻線譜成份、后槳葉頻線譜成份以及前、后螺旋槳共同作用下的線譜成份。該線譜成分可用公式fk，m=kB1f1+mB2f2，其中B1、B2分別為前、后槳葉片數(shù)、f1、f2分別為前、后槳軸頻，為今后利用該特點(diǎn)對(duì)水下高速目標(biāo)進(jìn)行判別打下了基礎(chǔ)。

[1]Kubota A,Kato H,Yamaguchi H.Unsteady structure measurement of cloud cavitation on a foil section using conditional sampling technique[J].Journal of Fluids Engineering,1989,111(2):204-210.

[2]王波,張永坤,高霄鵬.基于模型試驗(yàn)的螺旋槳空泡噪聲計(jì)算[J].海軍工程大學(xué)學(xué)報(bào),2010,22(6):87-90. Wang Bo,Zhang Yongkun,Gao Xiaopeng.Calculation of ship propeller cavitation noise based on model experiment[J]. Journal of Naval University of Engineering,2010,22(6):87-90.

[3]劉登成,洪方文,張志榮等.螺旋槳片狀空泡的CFD分析[J].艦船科學(xué)技術(shù),2009,31(1):43-46. Liu Dengcheng,Hong Fangwen,Zhang Zhirong,et al.The CFD analysis of propeller sheet cavitation[J].Ship Science and Technology,2009,31(1):43-46.

[4]楊瓊方,王永生,張志宏等.非均勻進(jìn)流對(duì)螺旋槳空化水動(dòng)力性能的影響[J].水動(dòng)力學(xué)研究與進(jìn)展,2011,26(5):538-550. Yang Qiongfang,Wang Yongsheng,Zhang Zhihong,et al.Calculation of highly-skewed propeller’s load noise using BEM numerical acoustics method in frequency domain[J].Acta ArmamentarⅡ,2011,32(9):1118-1125.

[5]胡健,馬騁,黃勝等.螺旋槳尾流場的數(shù)值計(jì)算[J].哈爾濱工程大學(xué)學(xué)報(bào),2008,29(12):1255-1260. Hu Jian,Ma Cheng,Huang Sheng,et al.Numerical analysis of flow fields behind propellers[J].Journal of Harbin Engineering University,2008,29(12):1255-1260.

[6]Lindau J W,Boger D A,et al.Propeller cavitation breakdown analysis[J].Journal of Fluids Engineering.2005,127:995-1002.

[7]朱德祥,張志榮,吳乘勝等.船舶CFD不確定度分析及ITTC臨時(shí)規(guī)程的初步應(yīng)用[J].水動(dòng)力學(xué)研究與進(jìn)展,2007, 22(3):363-370. Zhu Dexiang,Zhang Zhirong,Wu Chengsheng,et al.Uncertainty analysis in ship CFD and the primary application of ITTC procedures[J].Journal of Hydrodynamics,2007,22(3):363-370.

[8]Boswell R J.Design,cavitation performance and openwater performance of a series of research skewed propellers[R].Naval Ship Research and Development Center,Wasington,DC,Report 3339,March 1971.

[9]Kim S E.Multiphase CFD simulation of turbulent cavitating flows in and around marine propuslors[R].David Taylor Model Basin,U.S.A.,June 1982.

Numerical analysis for sheet cavitation noise characteristic of countra-rotating propeller

WANG Shun-jie,WANG Yi-chuan,DAI Wei-guo,LI Bin
(Navy Submarine Academy,Qingdao 266044,China)

In order to study cavitation noise characteristic of contra-rotating propeller,the numerical analysis of counter-rotation propeller cavitation was accomplished by using Schnerr-sauer cavitation model, RNG k-ε turbulence model and dynamic mesh model.The open-water performance was done with DTMB4381 propeller cavitation,the predicted cavity shapes were in agreement with both the experimental observation and numerical simulation result on open literature.And,by calculating and analyzing the contra-rotating propeller cavitation,the sound pressure and power spectrum density of different monitors was solved.The result shows that line-spectrum which is related to the blade frequences of front and rear propeller existing in the power spectrum and the non-uniformity of rear propeller flow is more obviously than front propeller.

sheet cavitation;contra-rotating propeller;numerical analysis

U661

A

10.3969/j.issn.1007-7294.2014.07.006

1007-7294（2014）07-0778-08

2013-12-10

國家自然科學(xué)基金資助項(xiàng)目（60901054）

王順杰（1984-），男，海軍潛艇學(xué)院博士研究生，E-mail：xiaojie000@163.com；王易川（1976-），男，海軍潛艇學(xué)院副教授。