基于面元法的船舶螺旋槳附連水質量與阻尼計算方法研究

鄒冬林 ， 張建波 ， 塔 娜 ， 饒柱石 ，3

（1.上海交通大學 振動、沖擊、噪聲研究所，上海 200240；2.上海交通大學 機械系統與振動國家重點實驗室，上海 200240；3.高新船舶與深海開發(fā)裝備協同創(chuàng)新中心，上海 200240）

0 引 言

船舶推進軸系是船舶水中航行的動力系統，也是船舶組件中非常重要的部件。其由螺旋槳、軸及軸承等組成。由于螺旋槳通常工作在不均勻的流場中，且其總是存在質量偏心，所以推進軸系既受到流體的隨機激勵又受到不平衡載荷的周期性激勵。當激勵頻率與軸系某階固有頻率接近時，會引起軸系強烈的共振，導致船舶振動異常并可能引發(fā)安全事故。因此對推進軸系的動力學特性分析一直以來都是國內外研究的熱點[1-5]，而準確預測軸系的固有頻率和模態(tài)振型是軸系動力學分析中的很重要的工作，同時也是難點。這是因為其影響因素多而且復雜，比如支承特性、聯軸器特性及螺旋槳質量等等。目前在大多數推進軸系動力學分析中，常常把螺旋槳簡化為集中質量。文獻[6]指出對推進軸系固有頻率及模態(tài)振型計算不夠準確有部分原因是無法準確估計螺旋槳的“質量”。這里所指的質量除了螺旋槳本身的金屬質量外還包含了螺旋槳與流體間相互作用的“虛擬質量”。這是因為螺旋槳工作在水介質中，水對螺旋槳產生復雜的流體反作用力。研究表明有部分流體反作用力和螺旋槳的振動加速度成比例，工程上稱之為“流體附加質量”或“附連水質量”。有資料顯示，部分螺旋槳縱向流體附加質量接近螺旋槳的金屬質量，因而不能忽略[7]。因此對螺旋槳質量的準確估計關鍵取決于附連水質量的計算，因為對于金屬質量而言，螺旋槳一旦被設計好后，基本就確定下來。而要正確預測螺旋槳的附連水質量，需要掌握流體與螺旋槳之間相互作用方式及壓力生成機理。目前準確計算螺旋槳附連水質量仍是一個難點。類似地，流體作用力中與速度成比例的部分稱之為“流體附加阻尼”。附加阻尼在軸系動力學分析中也是很關鍵的參數，因為其可以抑制軸系的動態(tài)響應，尤其在共振區(qū)時。

綜上所述，對于推進軸系設計工作者而言，要想準確預測軸系的固有頻率和模態(tài)振型，當務之急是如何準確且快速計算螺旋槳附加質量與附加阻尼。這是一項非常有必要和意義的工作，也是一項挑戰(zhàn)。而本文的工作正是從這一問題出發(fā)的。

目前，對于結構體或螺旋槳在流場中附加質量與附加阻尼的處理上，國內外都有不少文獻。概括起來主要分為三類：第一類是經驗估計法[8-10]。此類方法用一個系數乘以螺旋槳金屬質量來估計附加質量，比如文獻[8]中規(guī)定縱向振動附加質量為金屬質量的60%至100%。文獻[9]中規(guī)定回旋振動附加質量為金屬質量的10%到30%。這類方法操作簡單，所以目前仍被廣泛采用。但是由于其沒有考慮螺旋槳幾何形狀的影響，也沒有考慮不同自由度間的質量耦合影響，因而精度很差，而且不能估計附加阻尼的影響。第二類是實驗法。該方法通過測試軸系在空氣中與水中固有頻率的變化，從而反求出螺旋槳附加質量。如Burrill[11]對大量螺旋槳進行扭轉和軸向固有頻率測試，從而反求出了扭轉與縱向方向的附加質量。Wereldsma[12]和Wijk[13]也做了類似實驗。此類方法通常難以考慮不同自由度間的質量耦合效應，同時受結構阻尼影響，因而難以正確估計流體的附加阻尼。第三類是數值計算法，其中又可以分為兩種思路。第一種思路是通過力反求的方法。假設螺旋槳為一個單自由度系統，在水中以某一確定頻率做簡諧運動，通過數值方法計算出作用在螺旋槳上合力隨時間的變化規(guī)律，再經過FFT變換到頻域，找到對應頻率下的力的幅值與相位，由單自由度理論知識反求出附加質量與附加阻尼。而在計算水動力時通常采用基于勢流體的升力線理論、升力面理論、面元法或更為復雜的基于粘流體的CFD方法。如Hylaride[14]采用升力面理論研究了螺旋槳的附加質量與附加阻尼，其中螺旋槳水動載荷通過升力面理論計算得到。類似地，Schwanecke[15]和Vassilopoulos[16]也采用升力面理論計算了螺旋槳的附連水質量與阻尼。Salehyar[17]分別采用葉元體理論和面元法求解風力發(fā)電機葉片的縱向附加質量與附加阻尼，也是通過力來反求。Gaschler[18]采用面元法對空化螺旋槳進行了附加質量與附加阻尼的近似計算，其中水動載荷通過面元法計算得到。Esch[19]求解了螺旋槳縱向與扭轉方向附加質量與阻尼，考慮尾渦影響，研究表明尾渦對附加質量阻尼有影響，其中螺旋槳水動力通過CFD方法求解。此類方法原理上都屬于頻域法，即使規(guī)定螺旋槳按某一確定頻率簡諧振動，載荷中也會存在倍頻分量，而計算時忽略了載荷的高階倍頻分量，因而帶來一定誤差。因此當假設不同的簡諧頻率運動時，得到附加質量與附加阻尼均不同，附加質量與附加阻尼隨頻率變化。Blevins[20]在其專著指出當不考慮流場壓縮性及粘度時（勢流體），附加質量僅與結構的形狀有關。因此該結論與上述方法的結果有出入，這是由于忽略載荷中的高階倍頻分量導致的。數值方法第二種思路是通過計算流體動能的方法，該方法是計算附加質量的經典方法，在很多教科書中均有提及。該方法基于勢流理論，定義附加質量是由于結構運動導致流體動能變化引起的，從而通過計算出流體動能后求出附加質量。如Ghassemi[21]利用面元法，在槳葉表面分布源與渦，計算了其附加質量矩陣。其中忽略了尾渦影響。類似地，Hutchison[22]利用面元法法研究了導管螺旋槳的附加質量，也沒有考慮尾渦的影響。一方面對于螺旋槳葉片這種升力體結構，如果不考慮尾渦的影響，會造成一定的誤差[19]。另一方面這種基于動能的經典方法不能計算水動力引起的附加阻尼。目前，國內研究螺旋槳附連水質量的公開文獻非常少。郭益民[23]利用Parsons提出的近似公式計算了螺旋槳縱向與扭轉方向附加質量，屬于上文分類方法中第一類經驗估計法。李泓運等[24]利用有限元軟件中的聲單元計算了螺旋槳葉片彈性變形下的附連水質量，忽略軸的存在，沒有考慮尾渦影響。屬于數值方法中的基于動能的方法。另外還有不少針對其它結構附連水質量的研究文獻。文獻[25-27]用的是第三類方法第二種思路中基于動能的經典方法，而文獻[28]是第三類方法第一種思路中基于力反求的方法。

綜上所述，目前國內外針對螺旋槳附連水質量的計算分為經驗估計、實驗測量和數值研究三大類。經驗估計誤差太大，實驗測量成本高，且兩種方法均不能估計耦合質量與附加阻尼。數值研究中，部分是從頻域上基于力反求附連水質量和附加阻尼，此種方法求解精度不夠；部分是從計算流體動能出發(fā)，此種基于動能的經典方法不能計算附加阻尼，而且大部分都沒考慮尾渦影響。針對這些不足，本文基于勢流理論，利用面元法構建一種新的求解螺旋槳附連水質量的方法，該方法不僅可以求解附連水質量，還可以求解附加阻尼，同時也考慮了尾渦的影響。在求解附連水質量上，本文證明了該方法與基于動能的經典方法是等價的。文中最后給出幾個算例，證明本文方法的合理性。

1 面元法簡介

目前已有多種成熟方法預報螺旋槳水動力，如升力線理論、升力面理論、面元法（又稱邊界元方法BEM）、CFD法等。升力線理論含有大量假設，因而只適用于輕載螺旋槳。升力面理論雖然可以較為準確地預報總推力和總扭矩，但由于沒有考慮螺旋槳葉片的厚度，因而其預報的槳葉面壓力分布不夠準確。面元法考慮了葉片厚度的影響，因此不僅可以準確計算總推力與總扭矩，而且可以較準確地計算槳葉表面的壓力分布。同時相比于早期的面元法，目前又有了不少改進，應用范圍也越來越廣泛[29-32]。CFD方法盡管也能準確計算槳葉表面壓力分布，但是由于其計算量大，耗費時間長，因而應用并不廣泛。因此本文基于面元法相關理論來構建螺旋槳附連水質量及附加阻尼的計算方法。

面元法分為基于速度和基于速度勢的。本文采用較為簡捷的基于擾動速度勢的面元法。由于篇幅限制，本文對面元法的描述從簡，詳細請見文獻[33-34]。圖1是螺旋槳軸變形前后示意圖。建立兩套坐標系：OXYZ為慣性靜止坐標系，oxyz為隨體坐標系，附在葉片上。假設軸末端位移表示為：相應的速度可以表示為：

假設：（1）流體無粘、無旋且不可壓縮；（2）螺旋槳浸水足夠深，即不考慮自由液面影響，同時流體域延伸到無限遠；（3）不考慮空化影響；（4）尾渦形狀預先假設，即為線性尾渦。取一足夠大的外部控制面將其封閉在內。如圖2所示。

圖2 螺旋槳及周圍流場示意圖Fig.2 Propeller and fluid around it

流域的邊界面由物面SB，尾渦面SW和外邊界面S∞組成。在該流場中可用擾動速度勢Φ來表示螺旋槳的擾動，在oxyz坐標系中，Φ滿足Laplace方程：

同時擾動速度勢在邊界面S的每一部分上，還要滿足以下邊界條件：

（1）當外控制面距螺旋槳表面極遠時，其上的擾動速度趨于零，即

（2）假設尾渦面厚度為零，兩邊沒有速度跳躍和壓力跳躍，即

式中：Q1是尾渦面上的點；“+”和“-”分別代表尾渦面上、下表面。

（3）在物面上滿足法向速度為零的運動邊界條件，即

式中：Vin=V0+ω×r為不考慮軸振動時的進流速度，V0表示來流速度，ω表示螺旋槳旋轉速度；Vb表示由軸系振動引起的葉片速度；nQ是邊界面上的單位法向量。

在oxyz坐標系中，非定常Bernoulli方程可表示為：

式中：P0為參考點處流體壓力；ρ為流體密度；V=Vin+▽Φ，為總擾動速度。

由于擾動勢Φ滿足線性疊加原理，由（4）式可知，擾動勢可以分成兩部分，即

式中：φ由Vin進流速度產生；φ由Vb產生。因此相應地可以分解成兩個問題：

由于▽φ?Vin、▽φ?Vin且Vb?Vin，忽略高階小量，則相應的非定常Bernoulli方程可分解為

式中：V1=Vin+▽φ。

對于（7）式和（9）式，可以理解為求解剛性螺旋槳在非均勻流場中的水動力，屬于常規(guī)問題，很多文獻[33-34]詳細說明了其求解過程。由于本文采用線性假設，因此對進流速度引起的擾動勢φ與軸系振動速度引起的擾動勢φ的求解是兩個相互獨立的過程。換句話說擾動勢φ對附加質量及附加阻尼不會產生直接影響，不需要考慮（本文只針對Vin引起葉片變形很小的情況，在Vin導致葉片變形很大時，此時葉片的幾何形狀發(fā)生改變，從而會導致附加質量與附加阻尼產生相應變化，此種情況超出本文研究范圍。）。因此，接下來本文只討論求解由（8）式和（10）式所構成的問題，因為附加質量與附加阻尼只與這兩式有關。

因為（8）式滿足 Laplace 方程，根據 Green 定理，當場點P（x，y， ）z在物面上時，擾動速度勢可以表示為（結合邊界條件）：

式中：Δφ（Q1）為通過尾渦面的速度勢跳躍，可以記為

對于復雜形狀結構，（11）式很難求解析解，因此采用數值面元法求解。假設將螺旋槳一個葉片及相應輪轂劃分成N個四邊形面元，將該葉片泄出的尾渦劃分成NW個面元。如圖3所示。

（11）式可離散成：

圖3 螺旋槳的面元分布Fig.3 Panel arrangement of propeller and wake

式中：δij為 Kronecker函數；Cij，Wil和Bij為影響系數。定義如下：

式中：Z表示葉片數。（12）式寫成矩陣形式：

使用Morrino Kutta條件，因此Aij值定義如下：

式中：“±”表示隨邊處上下面元。

（14）式進一步可以表示為：

2 附連水質量與阻尼求解

以第i個面元為例推導附連水質量與阻尼。設第i個面元中心坐標為 （xi，yi，zi），則其由軸末端的平動與轉動引起的合速度為：

則有

由（16）式和（18）式可得：

同時有：

式中：Hij為影響系數，為速度影響系數，文獻[33]均給出了其解析公式。

將（20）式、（21）式代入（10）式得：

所有面元產生的合力與合力矩為：

由（24）式可得附連水質量與阻尼矩陣分別為：

由此可見，[m]和[c]是一個6階方正，代表了6個自由度。非對角線上不為零的數表示不同自由度間質量與阻尼具有耦合效應。以質量矩陣為例，文獻[7]給出其一般形式為：

由此可見，附連水質量矩陣是對稱的。其中m32與-m32等中的負號是由于右手坐標系下角度正負差異引起的[7]。

3 等價性證明

前文提及過計算附連水質量矩陣的經典方法（第二類思路）是從動能角度出發(fā)。此類方法無法求解附加阻尼矩陣。本文以文獻[21，35-36]中的算法為例，闡述通過計算流體動能而求解附連水質量矩陣的經典方法。最后證明在求解流體附加質量矩陣上，本文所構建的方法與基于動能的經典方法是等價的。

設由于物體運動導致的流體動能為：

上式變換中應用了Green定理。因為剛體空間運動有6個自由度，因此將速度擾動勢φ線性分解為6個分量，對應這6個自由度。設第i個面元φ可分解為：

由（8）式、（18）式和（28）式可知：

將（28）式和（29）式代入（27）式，經化簡（詳見文獻[19]），得到質量矩陣[m1]為：

（30）式與（25）式第一式對比可知，兩者完全等價。 因為由（20）式和（28）式可知｝恒成立。

因此本文所提方法與基于動能的經典方法在求解附連水質量矩陣時完全等效。而本文所提方法的一個優(yōu)勢是可以計算流體引起的附加阻尼。

4 數值結果及分析

為了驗證本文方法正確性，以幾種結構為分析對象，與理論結果或其它文獻的計算結果作對比，證明本文方法的有效性。

算例1：球體。半徑為0.1 m的球浸潤在無限流場中，由于球的幾何形狀很規(guī)則，因此附連水質量有解析解。將球表面劃分成1 200個四邊形面元，由于其為非升力體，因此可以不考慮尾渦的影響，面元分布如圖4所示。計算結果如（31）式所示。

圖4 球面元分布Fig.4 Panel arrangement of sphere

由于其對稱形狀，只有m11、m22和m33有值且相等，其余均為零，這與理論解結果一致，文獻[35]給出其理論解：

為了驗證網格收斂性，將球表面分別劃分成三種不同數目面元（1 200，2 400，4 800），結果如表1所示。

表1 面元法與理論解結果比較Tab.1 Comparison of the results between the panel method and analytical results

由表1可知，數值解與理論解之間的相對誤差均很小；隨著面元數目的增加，相對誤差越來越小，數值解收斂于理論解，表明了本文方法的正確性。

算例2：橢球體。設x軸半徑為1 m，y軸與z軸半徑均為0.2 m。劃分7 000個面元。如圖5所示。計算結果如（33）式所示。

由于其對稱形狀，有m22=m33，m44=0及m55=m66，其余均為零，這與理論解結果一致，理論解由文獻[35]給出，如表2所示。

表2 面元法與理論解結果對比Tab.2 Comparison of results between the panel method and analytical results

由表2可知，數值解與理論解之間的相對誤差均很小，表明了本文算法的正確性。

圖5 橢球體面元分布 Fig.5 Panel arrangement of ellipsoid

圖6 螺旋槳及尾渦面元分布Fig.6 Panel arrangement of propeller and wake

算例3：螺旋槳。以DTRC4119螺旋槳為例，該螺旋槳無側斜，無縱傾，3葉片，質量約為4.1 kg。螺旋槳主要幾何參數詳見文獻[33]。假設旋轉速度ω=1 rad/s，進給速度J=0.833。將槳葉表面及尾渦共劃分成1 065個面元，如圖6所示。

（34）式和（35）式分別為計算的附加質量矩陣與附加阻尼矩陣（均為國際標準單位）。

由此可以看出，不考慮耦合效應，主對角線上縱向附加質量約占金屬質量的51%，橫向（垂向）附加質量約占金屬質量的22%，這基本都在文獻[8-9]的經驗系數范圍之內。而文獻[8-9]中的經驗系數范圍過于寬廣，實際選取中完全憑經驗，難以操作，而本文的方法準確地給出了數值。

文獻[22]給出了Parsons和Vorus采用升力面理論計算螺旋槳附加質量矩陣的結果，如（36）式所示：

比較（34）式和（36）式，可得兩者相對誤差的絕對值如（37）式所示：

由（37）式可知附加質量矩陣結果相差不大，最大相對誤差為11.1%。主要原因有：（1）對槳葉表面進行四邊形面元劃分時帶來的離散誤差；（2）Parsons和Vorus用的是升力面理論，前文中指出升力面理論由于沒有考慮螺旋槳葉片的厚度，其預報的槳葉面壓力分布不如面元法準確，因而升力面理論本身存在理論誤差。盡管如此，仍然可以看出，本文方法與Parsons和Vorus的方法誤差均在工程接受范圍之內，因此進一步證明了本文方法的有效性。由于對螺旋槳附加阻尼的計算沒有找到相關參考數據，因此本文沒有對其做比較。

5 結 論

面元法是螺旋槳水動力分析中常用的方法。本文利用面元法，構建了螺旋槳隨軸系在水中振動時的附連水質量與阻尼數值計算方法，并且證明了本文所提方法與經典的基于動能的方法是完全等價的?；趧幽艿慕浀浞椒o法求解附加阻尼，而本文的方法可以方便求解附加阻尼。最后以球體、橢球體及螺旋槳為對象給出幾個算例，并與解析解或其它計算結果比較，附加質量結果吻合良好，證明本文方法的有效性。同時本文方法可以方便擴展到船體、水翼、槳舵及導管螺旋槳等結構的附加質量與阻尼的計算，因此該方法在船舶與海洋工程結構物的流固耦合分析中具有較為廣闊的應用前景。