基于Lighthill聲類比理論的噴水推進泵不穩定流動誘導噪聲的數值研究

張德勝， 張乃舒， 許 彬， 趙睿杰， 高雄發， 李 寧

(1.江蘇大學 國家水泵及系統工程技術研究中心，江蘇 鎮江 212013； 2. 噴水推進技術重點實驗室，上海 200011)

噴水推進泵是一種新型的動力裝置，主要是通過泵噴出水流的反作用力來獲得前進的動力，在潛艇、高速船舶等方面得到了廣泛的應用。噴水推進泵在實際運行時環境變化復雜，在非標準設計工況下運行時，會出現葉輪與導葉、泵體之間的動靜干涉，以及旋轉失速等情況，從而產生強烈的壓力脈動，導致異常的振動與噪聲[1-3]，影響泵的性能以及艦船行駛的穩定性與隱蔽性。為了保證艦船在各種運行環境下都能穩定運行并且降低噪聲輻射強度，噴水推進泵內部流動的特性以及輻射噪聲的控制便值得進一步進行研究。

對于泵的壓力脈動以及流動誘導的噪聲方面，國內外學者展開了一系列研究。Tsai等[4]研究了不斷運動的物體在流動介質內噪聲測量的問題，從而推導出了廣義Lighthill方程。司喬瑞等[5]等采用了聲學有限元法(finite element method，FEM)與邊界元法(boundary element method,BEM)兩種的方法，對離心泵流動誘導的噪聲進行了預測，結果表明BEM法能快速測出葉頻及其整數倍頻的噪聲大小，但無法求解過于復雜的邊界條件。FEM法對于復雜模型的預測與實際更加吻合，適用于寬頻噪聲研究，但計算時間過長。張德勝等[6]采用瞬態微型壓力傳感器的實驗方法得出不同流量工況下，在葉輪進口的動靜干涉的壓力脈動周期性相似規律最強，葉輪出口的波動規律相似性次之。尹江南等[7]采用軟件LabVIEW和MATLAB對試驗數據進行采集分析，結果表明在小流量時，軸頻、葉頻以及倍軸頻、倍葉頻是頻譜的主要頻率。張明宇等[8]采用了DES(detached-eddy simulation)分離渦模型對噴水推進泵進行數值模擬，結果表明壓力脈動頻率取決于葉輪轉動頻率,壓力脈動幅值沿輪轂到輪緣逐漸增大。Dong等[9]通過研究表明改變葉輪與壓水室之間的間隙及增大材料剛度的方法，可以實現對噪聲的降低。GonzLez等[10]通過實驗驗證表明泵內部的壓力脈動是引起泵振動的主要原因。董亮等[11]在試驗與模擬結合的基礎上對泵做透平的葉片安放角進行改變，從而控制其外聲場噪聲的水平。

對于泵內部壓力脈動誘導的振動噪聲研究，目前主要集中于離心泵[12-14]、軸流泵[15-17]以及貫流泵中，并且對于噴水推進泵聲振耦合的研究少有提及。本文以混流式噴水推進泵為模型，采用CFD/CA(computational fluid dynamics/computational acoustics)相結合的方法，先進行數值模擬計算，并將得出的外特性數值與實驗數據進行對比。同時對泵內部的壓力脈動特性及不穩定流動進行研究。將非定常計算中由葉片產生的旋轉偶極子和導葉產生的固定偶極子導出，對比分析邊界元BEM法和有限元FEM法在噴水推進泵內聲場計算中各自的優劣，并進行葉輪內部聲場及外聲場聲振耦合計算，對比不同工況下內聲場聲壓級與外聲場聲功率水平，結合內部流動情況與壓力脈動數值模擬，揭示了噴水推進泵內部流動與其聲場的關系，對于預測泵流致噪聲以及對其噪聲水平的優化，具有重要的指導意義。

1 數值模擬方法

1.1 流動控制方程

在非定常狀態下的流動聲學問題中，N-S方程可作為基本控制方程，來描述噴水推進泵內部流場。

(1)

(2)

(3)

1.2 湍流模型

為考慮葉片吸力面流動分離的情況，本次計算采用的湍流模型為SST(shear stress transfer)k-ω模型。該模型考慮了湍流剪應力的傳播情況，采用混合函數的方法，將標準情況下的模型與k-ε和k-ω模型結合起來，從而能夠精確預測逆壓梯度下如流動分離等復雜的流動情況，適用于精度較高的數值模擬。該模型輸運方程為

(4)

(5)

1.3 聲振耦合理論與Lighthill聲類比模型

固體結構與流體介質接觸時，會受到流體的激勵而發生振動，從而使固體結構周圍流場發生變化。流場的變化又會進一步使固體結構受到的載荷發生變化，這是常見的流固耦合問題。固體結構在水或空氣中受到介質的激勵作用形成振動產生聲音，流場中聲音的振動反過來對估計結構的振動又會產生影響，這便是聲振耦合[18]。噴水推進泵的內聲場及外聲場聲振耦合通過軟件LMS.Virtual Lab進行計算。對于噴水推進泵水下輻射噪聲的常見計算方法為聲學有限元法或聲學邊界元法，但在求解聲場輻射聲壓的問題上，FEM需要將整體的計算域進行離散，而BEM法只需提取結構表面的網格作為計算邊界，具有輸入數據量少、計算時間短的特點。聲音邊界元法分為直接邊界元和間接邊界元兩種算法，直接邊界元計算要求網格封閉，而間接邊界元計算的網格可以封閉也可以不封閉，同時對于計算過程中離散的格式、模型以及計算所需的邊界條件的要求都比直接邊界元要低，對低頻段范圍噪聲的求解具有明顯優勢，因此本文選擇更為簡便的間接邊界元法進行計算。間接邊界元法的本質為萊特希爾(Lighthill)聲類比法[19]，控制方程為

(6)

式(6)右端的3個分量分別表示四極子、偶極子、單極子聲源對泵內噪聲的貢獻量。偶極子聲源，也是本文對于噴水推進泵的聲場模擬中的主要聲源，是本文研究的重點。偶極子聲源主要表現為由于泵內的葉輪轉動，使得流體獲得能量，但由于葉輪是旋轉機械，因此對于流體能量的施加是不均勻的，并且具有周期性的，不同位置的能量或者壓力不同，并隨著葉輪的旋轉實時發生變化，其特性與壓力脈動的原理類似。當流體從葉輪獲得能量后，沖擊泵體，導葉，使固體域產生振動并向外界輻射出噪聲，這便是偶極子聲源在噴水推進泵中的特性。四極子聲源主要存在于研究物體在高速運動下導致湍流應力變化的情況中，在噴水推進泵這類速度較低的流體中的影響不大，在速度較高的流動問題中有較為明顯的貢獻值，例如在氣動聲學中，速度大于0.1 Ma時的問題，因此在本文中不做研究。在噴水推進泵中，偶極子聲源分為固定偶極子與旋轉偶極子兩種聲源，本文主要研究葉輪上的旋轉偶極子及噴水推進泵整體的外聲場聲振耦合。

2 噴水推進泵物理模型和計算域

2.1 泵設計參數及計算域

本文研究的對象為混流式噴水推進泵，主要設計參數為：設計流量Q=191.25 kg/s，葉輪直徑D2=320 mm，設計揚程H=4.75 m，轉速n=900 r/min，葉輪葉片數為Z=6，導葉數為Zd=7。

噴水推進泵模型通過UG 9.0建立，得到的水體計算域以及葉輪導葉的三維圖如圖1所示，該計算域由進水管、葉輪、導葉及出水管組成。

圖1 噴水推進泵計算模型Fig.1 Model of water jet pump

2.2 網格劃分及網格無關性分析

計算域通過ICEM進行網格劃分。該混流式噴水推進泵的葉輪間隙為單邊1 mm，采用包絡葉片的網格畫法，為了保證高質量的模擬計算，在近壁面、葉輪葉尖以及葉輪與泵體間隙層處進行局部網格加密，將葉輪的Y+的最大值控制在24以下，滿足計算所需的條件。五套無關性驗證網格的各部件網格數量,如表1所示，經過網格無關性驗證，最終確定當網格總數為970萬左右時，計算獲得的揚程變化波動小于2%，為獲得更好的網格質量以及更高的計算精度，本次計算的總網格數選用為970萬。噴水推進泵計算域葉輪、導葉的網格劃分，以及網格無關性驗證如圖2所示。計算用泵各部件的網格數量和質量,如表2所示。

圖2 計算域網格劃分及網格無關性驗證Fig.2 Mesh of computational domain and grid independence

表1 五組不同數量網格Tab.1 Mesh number of different parts

表2 各部件網格數及質量Tab.2 Mesh number and quality of components

2.3 邊界條件及前處理設置

為盡可能與試驗工況達到一致，本次計算采用有限體積法，湍流模型為SSTk-ω湍流模型，導葉與進出水管設置為無滑移壁面,葉輪進出口面動靜交界面為Frozen rotor。定常計算收斂精度設置為10-4。進口邊界條件設置為流量進口，出口邊界條件設置為自由出流。在非定常計算時，將定常計算的結果作為初始值，單次計算步長設置為5.5×10-4s，約為葉輪轉3°所需時長，且一次步長中計算次數最小設置為5次，最大迭代20次。為保證得到穩定的計算結果，將總的計算時長設定為1.46 s，約為葉輪轉22圈所需時長。過流部件固體壁面設置為無滑移邊界條件，且壁面粗糙度設置為光滑壁面，采用標準壁面函數解析黏性底層。動靜交界面設置為瞬態動靜轉子，求解器選擇high resolution，對流項使用中心差分方式，時間項選為二階歐拉差分方式[20]。同時設置葉輪每個旋轉周期輸出120個旋轉偶極子文件，導葉輸出120個固定偶極子文件，并將其中3圈的旋轉偶極子、固定偶極子各360個cgns文件作為聲場計算的壓力源文件。

為獲得噴水推進泵內部壓力脈動的情況，在葉輪的進口處、葉輪內部、葉輪與導葉交界面、導葉內部與導葉出口選取6個面設置監測點。從輪轂到輪緣，監測點編號分別為P1i～P6i。從葉輪進口邊到導葉出口邊，監測點編號分別為Pi1～Pi6。如圖3所示。取非定常計算的最后8個周期，作為壓力脈動頻域特性分析數據；選取非定常0.86～1.01 s內葉輪旋轉兩圈的數據作為壓力脈動時域分析數據。

圖3 監測點設置Fig.3 Monitor points

3 結果與討論

3.1 壓力脈動分析

通過定常計算得出了噴水推進泵不同工況下的外特性數據，同時該泵在荷蘭埃因霍溫科技大學完成了外特性實驗,試驗設施與樣機如圖4所示。通過數值模擬的結果擬合出的流量-揚程曲線與實際試驗數據相對比可知，數值模擬結果與試驗結果吻合較好，如圖5所示。根據圖5的流量揚程曲線可知，噴水推進泵出現了明顯的馬鞍區，馬鞍區最低點位于0.55Q處。在該處泵工作性能較差，泵內流動復雜，因此小流量工況下運行的內部流場及聲場是重點研究對象。

圖4 試驗設備Fig.4 The experimental facility

根據在噴水推進泵內設置的監測點得出的壓力脈動數據，可以分析泵在不同位置、不同工況下內部的流動情況。為了更加方便的分析不同位置與工況壓力脈動之間的關系，本文定義在進行分析時，運用無量綱壓力系數來進行描述。表達式為

(7)

式中：pi為任意時刻監測點的壓力數值，Pa；pave為該監測點在所取的2個轉動周期內的壓力數值平均值，Pa。

圖5 模擬與試驗數據對比Fig.5 Comparison of numerical data & experimental data

圖6為葉輪進口邊監測點P11～P16的壓力幅值分析。由圖可知，2個周期內共有12個波峰，單個周期6個波峰，與葉輪的6葉片數相符。壓力幅值曲線呈現明顯的周期性波動，并且從葉輪輪轂P11處到輪緣P16處壓力脈動幅值逐漸增大。

圖6 葉輪進口面壓力系數特性Fig.6 Pressure coefficient at inlet of the impeller

圖7為1.0Q工況(Q=191.253 kg/s)下，由葉輪進口到導葉出口6個不同截面上同取5號監測點的壓力脈動時域圖。由圖7(a)可知，葉輪進口的壓力脈動幅值最大，壓力脈動的幅值由葉輪進口到葉輪出口逐漸減小。由圖7(b)可知，在導葉中壓力脈動也呈現出周期性變化規律，但波動的規律性不如葉輪中的規律性強，并且壓力脈動幅值依然符合從進口到出口依次減小的規律。

圖8(a)為葉輪中監測點P15處的壓力脈動頻域分析。由于2個周期內的數據無法很好的體現低頻范圍內的壓力系數波動，因此取非定常計算最后8個周期作為分析所用的數據。由圖8(a)可以看出，1.00Q與0.75Q工況下壓力脈動系數的幅值都出現在一階葉頻(90 Hz)處，并且在1.00Q工況與0.75Q工況時，二階和三階葉頻處也出現明顯的峰值，在低頻處軸頻附近也出現較低的峰值。在0.55Q工況下，一階葉頻處也出現壓力系數的最大值，但在低頻段內的軸頻附近，出現了壓力系數的又一峰值，接近一階葉頻處的峰值。在高頻范圍內，1.00Q與0.75Q能看到二階、三階葉頻處的明顯峰值，但是0.55Q工況衰減明顯，在二階葉頻以上幾乎無大范圍波動。壓力脈動系數的幅值在1.00Q工況下最小，0.75Q其次，0.55Q最大，約為1.00Q工況的1.7倍。圖8(b)為葉輪中監測點P35處的壓力脈動頻域分析。與葉輪進口P15相比，P35處的壓力系數幅值明顯減小，低頻段范圍變化更加復雜，但是幅值依然是0.55Q工況下壓力系數最大，1.00Q工況最小。1.00Q與0.75Q工況下變化趨勢與進口處相同，葉頻處都會出現壓力系數的峰值，相比之下0.75Q與0.55Q工況的低頻范圍振動比P15處更加明顯，同時0.55Q工況低頻范圍內波動的峰值已經大于葉頻處的值。經過對監測點壓力脈動的綜合分析，1.00Q與0.75Q工況下壓力系數最大值都出現在葉頻處，0.55Q工況壓力脈動大于另外2個工況，且低頻范圍內壓力脈動明顯增大，接近甚至超過一階葉頻處的壓力脈動幅值。

圖7 噴水推進泵內部壓力系數時域特性Fig.7 Time domain characteristics of pressure coefficient in the water jet pump

圖8 噴水推進泵內部壓力脈動系數頻域特性Fig.8 Frequency domain characteristics of pressure coefficient in the water jet pump

通過分析可知，小流量工況下壓力脈動幅值普遍高于正常工況，且在非標準工況下，靠近葉輪導葉交接面的監測點在低頻范圍內壓力脈動幅值十分明顯，遠離葉輪導葉交接面的監測點壓力脈動幅值在葉頻附近較為明顯。這是由于在小流量工況下，出現明顯的旋轉失速，流道中聚集著大量的失速渦團，如圖9(c)所示。通過圖9(d)觀察吸力面不同時刻的速度流線分布可知，從葉片進口邊至出口邊有渦的初生，發展與脫落，并不斷持續這一過程，吸力面持續發生流動分離，失速渦團阻塞流道，阻礙流體正常運動，產生異常的振動。這些流動特征相似的渦團在小流量工況下，會使葉輪中產生明顯高于額定工況時的壓力脈動幅值。由于葉輪出口流動受到導葉葉片的相干作用，在葉輪與導葉間的無葉區產生了明顯的葉片和導葉動靜相干作用誘導的低頻壓力脈動，且這種現象在小流量工況下更為明顯。而在標準工況下，通過圖9(a)葉片吸力面流線可知，內部流動非常穩定，葉片右側進口處有一處較小的渦，是由于葉片進口處產生的延伸到下一個葉片工作面處的泄漏渦，以及頂間隙處產生的泄露渦造成的。1.00Q葉輪流道流動分析如圖9(a)與圖9(b)所示。0.55Q葉輪流道流動分析如圖9(c)與圖9(d)所示，圖中通過其吸力面2D Streamline中的速度流線變化體現其渦的遷移與生成過程。其中葉片圖示右側位置為進口邊。

圖9 不同工況流動分析Fig.9 Flow analysis under different operating condition

通過對0.55Q工況下渦分析可知，流道存在大量失速渦團阻礙正常流動，并且在小流量工況下，通過靠近動靜交接面處監測點頻域變換分析可知其幅值出現在軸頻附近。分析導葉內的流場變化可知，由于導葉具有整流作用，導葉內渦的演化呈現出一持續變化的過程，可分為產生、發展到消失3個過程。分析非定常計算不同時刻的結果文件，整個過程大約以轉軸轉動360°為一個周期，也就是0.066 7 s，約為15 Hz，這一過程在葉輪的轉動周期內不斷持續生成，這是導致交界面附近監測點低頻范圍內高脈動幅值產生的原因，如圖10所示。

圖10 導葉流道中不同時刻流線分布Fig.10 Streamline in passage

3.2 聲場模擬與分析

在聲場計算中，聲源信息是由非定常計算中輸出的CGNS文件得來的。噴水推進泵中的噪聲源主要為偶極子聲源，偶極子聲源分為旋轉偶極子聲源與固定偶極子聲源。旋轉偶極子聲源主要是來源與旋轉的葉輪與流體相互作用形成的，在LMS軟件的葉輪內聲場計算中定義為扇聲源。固定偶極子聲源主要是從葉輪中流出的流體對導葉等壁面形成沖擊或動靜干涉形成的噪聲聲源。采用CFD與邊界元法結合的計算流程如圖11所示[21]

圖11 聲場計算流程圖Fig.11 Flow chart of acoustic calculation

首先將非定常計算生成的CGNS文件導入LMS中，本次聲場計算取最后兩圈，每3°記錄一次的方式，作為本次計算的聲源文件。新建另一個LMS文件，導入預先提取好的邊界元網格與場點網格，分別對兩種不同的網格類型進行定義，并設定流體材料及屬性。然后導入預先儲存好的聲源文件，將其定義為聲源，并將計算范圍設定在0～540 Hz，也就是9倍葉頻處。三種工況下內聲場的計算結果,如圖12所示。三種工況均在葉頻處出現聲功率的最高值，葉頻的倍頻處也出現明顯的峰值，且隨著頻率的升高，聲功率強度逐漸下降。1.00Q工況的聲功率強度在三種工況里最小，0.75Q與0.55Q工況下聲功率強度的峰值非常接近，但其余頻率下0.55Q工況的聲功率強度要強于0.75Q工況。

圖12 內聲場聲功率計算結果Fig.12 Acoustic power

通過對一階葉頻處(90 Hz)葉輪圓周面上的聲壓級云圖分析可知，將三種工況的標尺調節一致后明顯得出，1.00Q工況的聲壓級水平最低，0.55Q工況下聲壓級水平最高。三種工況下的葉輪軸面聲壓級都呈現出一定的規律性，由輪轂到輪緣逐漸增大，與壓力脈動的趨勢相符合。1.00Q標準工況呈現出較為規律的圓形分布，0.75Q出現了略微變形，但由于0.55Q馬鞍區工況下葉輪內流動情況復雜，出現不穩定流動，旋轉失速等情況，因此聲壓分布呈現出明顯的不規律的情況，如圖13所示。

對于外聲場聲振耦合的模擬，首先要進行噴水推進泵固體域的模態分析，得到不同頻率下的模態響應；同時將固體域模型導入LMS中生成網格并粗化，將粗化后的網格導入邊界元計算模塊，定義網格屬性及結構屬性，本次計算將泵體材料設定為不銹鋼。將泵與電機相連的面與泵的進出口法蘭面設為約束面。在聲振耦合中，定義兩種不同的場點網格。分別為面場點網格與指向性場點網格。為監測不同方向上聲壓輻射水平，設置軸向、徑向2個方向的面場點網格，以及以1 m為半徑，360°圓周方向上每隔10°設置一個監測點，共36個監測點，作為指向性場點網格，對計算結果進行監測。兩種場點網格的設置如圖14所示。

圖13 葉輪聲壓級云圖Fig.13 Sound pressure level

圖14 場點網格布置Fig.14 Field point mesh

以1.00Q工況為例分析，混流式噴水推進泵是對稱結構的旋轉機械，因此在徑向截面的場點網格上，聲場輻射云圖呈現出近似于圓的狀態，如圖15(a)所示。而在軸向截面的場點網格上，聲場輻射云圖呈現出明顯的偶極子特性，且極大值延葉輪徑向分布，極小值延軸向分布，如圖15(b)所示。出現這種偶極子分布的主要原因是在流動中有障礙物存在，流體與固體域發生碰撞同時產生不穩定的反作用力。在噴水推進泵中，流體在葉輪的作用下不斷與泵體發生沖擊，同時流體從葉輪進入導葉，導葉具有整流作用，不斷對流體產生反作用力，因此會出現以葉輪和導葉為圓心，并向葉輪徑向方向上輻射出極大值的分布現象。

通過指向性場點網格，可以分析在不同工況下外聲場聲壓輻射水平。如圖16所示，在一階葉頻處，葉輪軸向截面的指向性場點網格顯示出明顯的偶極子特性，聲壓級的極小值位于90°與270°附近，極大值位于0°與180°附近，這種極值分布與平面場點網格的聲壓級云圖所顯示的極值分布規律一致。隨著流量工況減小，聲壓級輻射水平隨之增強，呈現出大流量工況輻射強度小，小流量工況輻射強度大的特點。這個現象也與之前所分析的壓力脈動現象相符合。由于0.55Q工況位于馬鞍區內，運行不穩定，出現明顯的旋轉失速與動靜干涉，產生異常的激振力，因此聲壓輻射水平最強。該現象也能推斷出，噴水推進泵內部壓力脈動越大的工況，對應聲場的輻射聲壓級就越強，這對于研究噴水推進泵內部流動噪聲以及對噪聲水平的預測及優化，具有重要的指導和參考價值。

圖15 面場點網格聲場云圖Fig.15 Cloud image of sound pressure of plane field point

圖16 不同工況聲壓及指向分布Fig.16 Sound pressure level and distribution under different operation condition

4 結 論

本文針對艦船噴水推進泵不同運行工況下內部不穩定漩渦流誘導壓力脈動和噪聲問題，以某一型號混流式噴水推進泵為研究對象，采用CFD數值模擬與CA計算聲學相結合的方法，分析不同工況下壓力脈動時域、頻域變化規律，同時對噴水推進泵內聲場及外聲場聲振耦合進行計算，得出以下結論：

(1)噴水推進混流泵內部壓力脈動幅值從葉輪進口到出口逐漸減小，沿導葉進口到出口逐漸減小，同一軸向截面面上，壓脈動力幅值沿著輪轂到輪緣逐漸增大；尤其在葉輪和導葉交界無葉區，動靜干涉作用較強，該監測面的壓力脈動曲線出現了不規則的波動。

(2)在1.00Q與0.75Q運行工況下，噴水推進泵葉輪內部的壓力脈動頻域時均在一階葉頻處出現峰值；在0.55Q運行工況(馬鞍區)內運行時，一階葉頻也出現了峰值，但同時出現明顯的低頻脈動，且某些監測點低頻范圍內的峰值已經超過了一階葉頻的峰值，主頻不再是一階葉頻。三種工況下頻域變換的峰值最大值均出現在0.55Q運行工況，1.00Q工況下的壓力脈動峰值最小。

(3)在0.55Q與0.75Q和1.00Q三種運行工況下，噴水推進泵的內聲場聲功率強度趨勢相似度很高，葉頻及其倍頻處都出現峰值；且1.00Q工況下的聲功率強度最低，0.55Q與0.75Q工況下聲功率強度最大值接近，但其余頻率下0.55Q工況聲功率強度大于0.75Q工況。

(4)額定工況下噴水推進泵聲場輻射云圖在徑向截面上聲壓云圖趨近于圓，在軸向截面上，呈現出明顯的偶極子分布規律，且極大值位于徑向方向，極小值位于軸向方向上。通過指向性場點網格可知，噴水推進泵外聲場聲壓在不同工況下均呈現出明顯的偶極子特性，且極小值位于90°與270°附近，極大值出現在0°與180°附近。且壓力脈動程度越強的工況，呈現出的外聲場聲壓輻射水平也越強。