湯世昕， 沈育靜， 陳紀(jì)康， 段文洋

(哈爾濱工程大學(xué) 船舶工程學(xué)院，黑龍江 哈爾濱 150001)

迄今為止已有多種理論或方法計算船舶螺旋槳的水動力性能，如螺旋槳的升力線理論、升力面理論及面元法等。螺旋槳升力線理論借鑒了機翼升力線理論，其基本方法是用一條從葉根到葉梢的渦線近似代替螺旋槳葉模型。但是由于其在理論上過于簡化，不少學(xué)者相繼發(fā)展了升力面理論加以完善。升力面理論把船舶螺旋槳等假定簡化成為薄物體，即忽略了螺旋槳槳葉厚度的影響，因而該理論中的邊界條件與真實的物面條件并不相同。而且由于槳葉薄翼在導(dǎo)邊處的流體速度趨于無窮，使得該理論無法求解出位于導(dǎo)邊附近的壓力情況，限制了進一步對壓力分布的預(yù)報和優(yōu)化。面元法沒有對螺旋槳的幾何形狀做任何假設(shè)，并在實際槳葉面上滿足物面邊界條件，因此用面元法不僅可以準(zhǔn)確地預(yù)報螺旋槳的推力和轉(zhuǎn)矩，而且可以較準(zhǔn)確地預(yù)報螺旋槳表面的壓力分布。

面元法最早由Hess等[1]提出和應(yīng)用，用于求解三維不可壓、無粘和無旋流體中的無升力體水動力問題。Morino[2]提出了速度勢面元法，并將其應(yīng)用至升力體的水動力計算中。隨后Hoshino等[3-5]發(fā)展出了一系列基于速度勢面元法的數(shù)值計算方法，分析普通螺旋槳、大側(cè)傾螺旋槳以及導(dǎo)管螺旋槳的水動力性能。Maniar等[6-8]應(yīng)用B樣條曲線構(gòu)建升力體的幾何表面并使用高階面元法，預(yù)報結(jié)果更為準(zhǔn)確。此外董世湯等[9-12]國內(nèi)學(xué)者為推動面元法在工程中的應(yīng)用進行了研究以及數(shù)值驗證工作。

現(xiàn)有的基于速度勢的面元法稱為常值面元法(constant panel method, CPM)，它假定面元上的速度勢為常值，通常采用Morino方法計算影響系數(shù)，并使用Yanagizawa發(fā)展的插值方法獲得物體表面的速度分布，這種速度計算方式可以有效減少對影響系數(shù)求梯度消耗的計算時間，但也會不可避免地導(dǎo)致數(shù)值誤差。該誤差在槳葉導(dǎo)邊等幾何變化明顯的位置較大，并且最終導(dǎo)致常值面元法的精度下降[5]。

針對上述插值方法求解誘導(dǎo)速度的缺點，本文分別使用零階和一階泰勒展開邊界元方法(taylor expansion boundary element method, TEBEM)求解誘導(dǎo)速度，預(yù)報螺旋槳定常水動力性能，并與插值方法的結(jié)果和實驗值進行比較。隨著計算機性能的提高，可以通過對影響系數(shù)求梯度獲得物面速度，即零階TEBEM的求解誘導(dǎo)速度方法。一階TEBEM[13]對源偶混合分布方法中的偶極強度進行一階泰勒展開，并將誘導(dǎo)速度作為未知數(shù)，在線性方程組中直接求解，可以有效提高物面突變處計算精度，從而改善槳葉導(dǎo)邊、隨邊處壓力分布預(yù)報。

1 螺旋槳定常水動力的定解問題

1.1 坐標(biāo)系建立

對于在速度為VA的來流中以角速度Ω轉(zhuǎn)動螺旋槳，采用如圖1所示的固定于槳葉的坐標(biāo)系[3]o-xyz，原點o固定于槳盤面的中心點，x軸沿螺旋槳的旋轉(zhuǎn)軸指向下游，z軸沿螺旋槳某一葉片的母線，z軸與x軸和y軸組成右手坐標(biāo)系。同時也采用一柱坐標(biāo)系o-xrθ作為參考坐標(biāo)系，其θ角度量由z軸起始逆時針旋轉(zhuǎn)的角度[14]。

圖1 螺旋槳坐標(biāo)系Fig.1 Coordinate system of propeller

1.2 速度勢與邊界條件

流域的邊界面S由物面SB(包括槳葉表面和槳轂表面)、尾渦面SW和外邊界面S∞組成，且邊界的法向量指向域內(nèi)，如圖2所示[15]。

圖2 升力體及其周圍流場Fig.2 Body and flow field around it

(1)

同時滿足以下邊界條件：

1)在物面上滿足法向速度為零的運動邊界條件：

2)假設(shè)尾渦面的厚度為0并通過它沒有法向速度跳躍和壓力跳躍，即：

式中下標(biāo)+、-表示尾渦的上下面。

3)當(dāng)外邊界面極遠時,其上的擾動速度趨于0：

故擾動勢滿足的定解問題為：

(2)

將邊界條件式(2)代入式(1)，流場中場點P滿足的邊界積分方程可簡化為:

(3)

式中V0=VA+Ω×r。

1.3 壓力Kutta條件

式(3)中Δφ需要壓力Kutta條件來決定，即在三維升力體問題中，應(yīng)滿足隨邊處葉面和葉背間的壓力差為零：

假定沿螺旋槳槳葉徑向網(wǎng)格數(shù)為NR，對于定常流動問題，尾渦面上Δφ沿尾流方向上是不變的，且每個槳葉泄出的尾渦片是相同的。只要確定與螺旋槳隨邊相接的尾渦面上的值(Δφ)j，j=1,2,…,NR，即決定了整個尾渦面上的Δφ。如果根據(jù)等壓條件的離散格式：

(4)

同時令：

(ΔpT.E.)j=fj(Δφ1,Δφ2,…,ΔφNR)

j=1,2,…,NR

建立起對(Δφ)j的NR個方程。由于壓力是未知變量φ的非線性函數(shù)，使方程組不能直接進行求解，故只能通過牛頓迭代法求解。

每次迭代得到(Δφ)j的值，可通過定解問題求解速度勢，進而求解出物面上速度分布，再通過Bernoulli方程求出物面上壓力分布及隨邊處壓力差，重復(fù)上述流程直至式(4)成立。

1.4 螺旋槳水動力計算

利用Bernoulli方程：

(5)

槳葉壓力分布確定以后，通過物面壓力積分的方式得到螺旋槳所受推力和轉(zhuǎn)矩[4]：

式中：n=(nx,ny,nz)為單位法向量；r=(x,y,z)為位置矢量。并得到壓力積分離散形式為：

(6)

式中：ni=(nxi,nyi,nzi)為面元i控制點處單位法矢;ΔSi為面元i的面積;(xi,yi,zi)為面元i中心處坐標(biāo)。此外還須考慮粘性力對推力及轉(zhuǎn)矩的作用：

(7)

式中：vti=(vtxi,vtyi,vtzi)為面元i處局部速度;Cfi為粘性阻力系數(shù)，根據(jù)普朗特-許力汀平板摩擦阻力系數(shù)公式計算。因軸向坐標(biāo)向下游為正，從而總推力、扭矩為:

(8)

2 常值面元法原理簡述

將槳葉表面及槳轂表面總稱為物面SB，共劃分了N個面元，在尾渦面SW上共劃分了N1個面元。將式(3)離散化并令場點P沿面元法線方向趨近于原點，即為常值面元法所求解的離散方程，i=1,2,…,N。

(9)

式中Δφ由壓力Kutta條件決定。使用前文的迭代方法得到速度勢，并使用Yanagizawa發(fā)展的插值方法得到槳葉表面速度分布，計算壓力并得到水動力，具體的求解方法可參考文獻[3]。注意式(9)中影響系數(shù)的求法沒有采用Morino方法計算，具體的求法可參考文獻[15-16]，在下文中將這種求解誘導(dǎo)速度的方法簡稱為插值法。

除了技術(shù)人員下基層、到田間面對面講解油菜高產(chǎn)栽培技術(shù)外，農(nóng)業(yè)部門還利用街天、農(nóng)資打假等時機，加大宣傳引導(dǎo)，讓農(nóng)戶早知道，早受益。

3 泰勒展開邊界元法原理

3.1 一階TEBEM原理

(10)

式中i=1,2,…,N。

為了構(gòu)造一階TEBEM求解的邊界積分方程，需要將式(3)離散，并考察垂直于法線平面內(nèi)相互正交的2個方向上的一階導(dǎo)數(shù)，將偶極強度展開式代入后令場點P沿面元法線趨近于原點可得[16]：

(11)

3.2 零階TEBEM原理

零階TEBEM可作為一階TEBEM的特例，即零階TEBEM將偶極看作均勻分布，并通過式(9)直接求解得到速度勢。為了獲得誘導(dǎo)速度，對式(3)中場點P考察垂直于法線平面內(nèi)相互正交的2個方向上的一階導(dǎo)數(shù)，離散化后令場點P沿面元法線趨近于原點可得：

(12)

隨后根據(jù)式(5)計算壓力，通過式(6)～(8)得到水動力結(jié)果。零階和一階TEBEM中影響系數(shù)的求法同樣可參考文獻[15-16]。

4 網(wǎng)格劃分

4.1 槳葉劃分

為獲得離散方程的數(shù)值解，將槳葉表面沿著弦向和徑向劃分成一系列平面四邊形單元，并在導(dǎo)邊和隨邊及葉梢和葉根處進行網(wǎng)格加密。若R、RH分別為螺旋槳半徑和轂徑，弦向和徑向分別采用余弦分割[10]:

式中：si為葉剖面上點至導(dǎo)邊的弦向距離；cj為半徑rj處葉剖面弦長；Nc為弦向劃分次數(shù)；NR為徑向劃分次數(shù)，且NR1、NR2分別為[RH, 0.9R]、[0.9R,R]內(nèi)的徑向劃分次數(shù)，故NR=NR1+NR2。又有βci和βrj滿足：

4.2 尾渦形狀及劃分

在數(shù)值求解方程之前，首先確需定螺旋槳尾渦面的幾何形狀，由于尾渦面實際上事先未知，而且真實幾何形狀相當(dāng)復(fù)雜，一般是建立一個尾渦幾何模型來處理，用經(jīng)驗方法確定它的幾何形狀模型。

本文參考了文獻[9]中的尾渦模型，尾渦在槳葉隨邊處以螺旋槳幾何螺距角作為它的螺距角下泄，不考慮螺距沿軸向的變化，同時忽略徑向的收縮，即尾渦面是一個徑向變螺距的螺旋面。在劃分尾渦網(wǎng)格時，沿著螺旋線進行等距劃分，在靠近隨邊處進行網(wǎng)格加密。

5 數(shù)值結(jié)果與討論

5.1 計算模型

使用3種誘導(dǎo)速度計算方法預(yù)報了DTMB4119、KP458槳的敞水性能，并與實驗結(jié)果進行比較。螺旋槳的主要幾何參數(shù)列于表1中。

表1 螺旋槳的主要幾何參數(shù)Table 1 Geometrical parameters of propellers

5.2 網(wǎng)格收斂性分析

網(wǎng)格收斂性主要包括槳葉網(wǎng)格收斂性、尾渦網(wǎng)格收斂性和尾渦長度收斂性。其中槳葉網(wǎng)格數(shù)為NC×(NR1+NR2)，尾渦網(wǎng)格數(shù)為NR×NW(NW為沿螺旋線劃分次數(shù))。實際數(shù)值計算時，在尾渦面泄出一定距離后將其截斷；尾渦面長度LW指螺旋面沿x軸的長度，需要考慮其對數(shù)值結(jié)果的影響。

將DTMB4119作為計算模型，以一階TEBEM為例驗證數(shù)值結(jié)果的收斂性，表2～3給出了尾渦收斂性情況，表4給出了槳葉網(wǎng)格收斂性情況。首先在給定槳葉網(wǎng)格數(shù)(NC=28、NR1=18、NR2=5)后，通過改變尾渦劃分次數(shù)NW，在不同進速系數(shù)工況下，發(fā)現(xiàn)當(dāng)尾渦長度不小于1.4D(D為螺旋槳直徑)時，尾渦長度對計算結(jié)果幾乎無影響，故選擇尾渦長度LW為1.4D，尾渦網(wǎng)格數(shù)為23×58(NR=23,NW=58)。

表2 尾渦網(wǎng)格收斂性(J=0.7,NC=28,NR1=18,NR2=5)Table 2 Convergence of the wake panels (J=0.7,NC=28,NR1=18,NR2=5)

表3 尾渦網(wǎng)格收斂性(J=0.9,NC=28,NR1=18,NR2=5)Table 3 Convergence of the wake panels (J=0.9,NC=28,NR1=18,NR2=5)

表4 槳葉網(wǎng)格收斂性(J=0.7，尾渦網(wǎng)格數(shù)為23×58，尾渦長度1.4D)Table 4 Convergence of the blade panels (J=0.7， wake panels are LW=23×58，LW=1.4D)

對于槳葉網(wǎng)格劃分，在確定尾渦劃分方式后，發(fā)現(xiàn)當(dāng)槳葉網(wǎng)格數(shù)選擇NC=28、NR1=18、NR2=5時，計算結(jié)果基本收斂，故在隨后的計算中均采用這種槳葉劃分方式，面元分布情況及尾渦如圖3所示。

圖3 離散網(wǎng)格局部坐標(biāo)系Fig.3 Sketch of the local coordinate system

用同樣的方法分析插值法和零階TEBEM的收斂性情況，計算收斂時插值法選擇槳葉網(wǎng)格數(shù)為36×52(NC=36、NR1=42、NR2=10)，尾渦網(wǎng)格數(shù)為52×58(NW=58)，尾渦長度LW為1.4D。零階TEBEM擇槳葉網(wǎng)格數(shù)為33×36(NC=33、NR1=28、NR2=8)，尾渦網(wǎng)格數(shù)為36×58(NW=58)，尾渦長度LW為1.4D。

5.3 數(shù)值結(jié)果分析

本文3種數(shù)值方法對DTMB4119螺旋槳的計算結(jié)果與實驗值[17]及其他研究所的結(jié)果比較如圖4所示。根據(jù)目前的計算結(jié)果，對于推力系數(shù)，3種方法的結(jié)果相近，且精度相當(dāng)；隨著進速系數(shù)的改變，呈現(xiàn)出線性變化的特點；對于轉(zhuǎn)矩系數(shù)，零階和一階TEBEM的結(jié)果大于插值法的結(jié)果，且更接近實驗值，隨著進速系數(shù)的減小，3種方法計算的誤差會稍大一些。但總體來講，零階和一階TEBEM的結(jié)果與實驗值有較高的一致性，且在設(shè)計點附近的精度要稍高于其他的計算結(jié)果。CSSRC[9]與MHI在預(yù)報螺旋槳定常性能時采用了不同的面元形狀及影響系數(shù)計算方式，因此使用插值法計算速度得出的最終結(jié)果與本文用插值法結(jié)果有一定差異。

圖4 DTMB4119的敞水性能Fig.4 Open water performance of DTMB4119

圖5 DTMB4119網(wǎng)格示意Fig.5 Panel arrangement for DTMB4119

DTMB4119槳的槳葉表面壓力分布由Jessup在文獻[18]中給出，通過激光測速儀測得槳葉表面的速度分布，然后由Bernoulli方程換算得到表面壓力分布。該實驗結(jié)果一直被用來作為校驗計算方法的標(biāo)準(zhǔn)依據(jù)。圖6～7給出了4119槳在設(shè)計點J=0.833下槳葉表面壓力分布的實驗值與3種數(shù)值計算方法的結(jié)果比較。

圖6 r/R=0.7處壓力分布Fig.6 Pressure distribution at r/R =0.7

圖7 r/R=0.9處壓力分布Fig.7 Pressure distribution at r/R =0.9

通過0.7R處壓力分布與實驗值的比較，可以看到3種方法得到的結(jié)果與實驗值在大部分范圍內(nèi)吻合較好；使用零階和一階TEBEM得到的壓力系數(shù)在葉背及葉面接近導(dǎo)邊處變化程度較大，即流體速度發(fā)生了明顯的變化；從實驗結(jié)果中也可以觀察到相同的特征，通過插值法得到的壓力變化更加平滑，導(dǎo)致產(chǎn)生誤差而無法體現(xiàn)出上述變化特征。類似地，在0.9R處壓力分布的比較上，零階和一階TEBEM得到的壓力分布更接近實驗結(jié)果，因此在壓力積分的結(jié)果上，零階和一階TEBEM所得轉(zhuǎn)矩系數(shù)更大，且接近于實驗值。

圖8～9為3種方法所得壓力系數(shù)-CP分布的云圖，可以看出差異主要位于導(dǎo)邊和隨邊處，在其余的位置結(jié)果接近。

圖8 葉面壓力系數(shù)分布Fig.8 Pressure coefficient distribution on the upper surface

圖9 葉背壓力系數(shù)分布Fig.9 Pressure coefficient distribution on the lower surface

圖10～11中比較了3種方法在緊靠導(dǎo)邊的面元沿徑向的壓力分布情況；圖12中比較了緊靠隨邊的面元沿徑向的壓力分布情況。尾渦網(wǎng)格劃分時，3種方法均選擇NW=58、LW=1.4D；槳葉劃分時，均采用NC=32，NR1=21，NR2=6。零階和一階TEBEM得到的壓力分布在導(dǎo)邊、隨邊處有相似的變化趨勢，且一階TEBEM的壓力系數(shù)結(jié)果較小；插值法結(jié)果在導(dǎo)邊處的大部分范圍內(nèi)與前兩者都有明顯的差異，即插值法引起與實驗值的誤差；在接近葉梢的位置，3種方法的預(yù)報結(jié)果均向低壓快速變化，可能是未考慮空化的影響。同時參考圖6～7中0.7R、0.9R處導(dǎo)邊和隨邊的壓力分布實驗值，發(fā)現(xiàn)一階TEBEM結(jié)果的誤差更小，因此相比于零階TEBEM，局部壓力預(yù)報的精度更高一些。

圖10 葉面導(dǎo)邊處壓力分布Fig.10 Pressure distribution at the leading edge on the upper surface

圖11 葉背導(dǎo)邊處壓力分布Fig.11 Pressure distribution at the leading edge on the lower surface

圖12 隨邊處壓力分布Fig.12 Pressure distribution at the trailing edge of the blade

為了進一步分析3種誘導(dǎo)速度計算方法的求解精度，選擇了常見的KP458槳進行計算比較。相較于無縱傾、側(cè)傾的三葉槳DTMB4119，KP458槳在各槳葉剖面處具有不同程度的側(cè)傾，在幾何形狀上更為復(fù)雜。在尾渦網(wǎng)格劃分時，3種方法均選擇NW=58、LW=1.4D；槳葉劃分時，一階TEBEM采用NC=28、NR1=18、NR2=5；零階TEBEM采用NC=28、NR1=25、NR2=8；插值法采用NC=36、NR1=42、NR2=10。螺旋槳網(wǎng)格劃分示意圖如圖13所示，計算結(jié)果與實驗值的比較見圖14～15。采用零階和一階TEBEM來計算具有不同幾何特征的螺旋槳敞水性能，推力系數(shù)的精度相近，而轉(zhuǎn)矩系數(shù)的計算精度相比于插值法有明顯的提升，誤差減少近10%。

圖13 KP458網(wǎng)格示意Fig.13 Panel arrangement for KP458

圖14 KP458推力系數(shù)Fig.14 Calculated thrust coefficients for KP458

圖15 KP458轉(zhuǎn)矩系數(shù)Fig.15 Calculated toque coefficients for KP458

6 結(jié)論

1)尾渦長度對計算結(jié)果的影響很小，通過對不同工況的計算驗證，發(fā)現(xiàn)沿x軸向長度取約1.4倍直徑長度即可。

2)使用平面四邊形面元近似螺旋槳幾何形狀可以取得較好的收斂結(jié)果，3種方法中以一階TEBEM的收斂性最好。

3)采用零階和一階TEBEM計算誘導(dǎo)速度可以提高螺旋槳導(dǎo)邊及隨邊處壓力分布預(yù)報的精度，其中一階TEBEM具有更高的局部壓力預(yù)報的精度；兩者計算的轉(zhuǎn)矩系數(shù)與實驗值的誤差較插值法減少近10%。所得敞水性能曲線可進一步用于船舶航速及功率的評估。