基于Theodorsen方法的螺旋槳空泡數值計算研究

趙 鵬 周玉龍 趙鵬飛 楊 超

（江蘇科技大學 船舶與海洋工程 鎮江212003）

0 引 言

避免或推遲螺旋槳空泡的產生對螺旋槳設計意義重大。如今，應對螺旋槳空泡問題的主流方法是采用空泡斗校核［1-2］，該方法雖然簡單易行，但對于槳葉剖面是否會出現空泡，只能進行粗略判斷，無法細化到空泡產生的具體位置；而面元法［3-5］雖然能較好地預報槳葉上的空泡形狀與位置，但計算量相對較大，運算周期長［6］。

目前已廣泛應用升力線、升力面理論方法來設計螺旋槳，可初步得到一個槳葉剖面形狀（見下頁圖1）及各項有關性能參數，然而無法得到槳葉剖面上的壓力分布，因此不能直接對螺旋槳空泡性能進行預報。針對這一缺陷，本文使用Theodorsen法［7］及邊界層修正［8］，來計算槳葉剖面有關各點的壓力，將所得槳葉剖面各點壓力系數與該點空化數進行比較，從而對空泡進行初步預報，判斷槳葉剖面上空泡出現的位置，為螺旋槳設計者調整各項參數以改善空泡性能提供依據。主要內容如下：

圖1 初步提供的槳葉剖面

圖2 進行坐標調整后的槳葉剖面

（1）整理出由Theodorsen法和邊界層修正來預報空泡的整體思路；

（2）對有攻角的槳葉剖面進行坐標轉換，使其滿足Theodorsen方法的要求；

（3）對近似圓中坐標點所需進行的線性擬合的目標函數定義為三角函數；

（4）邊界層修正時，為便于計算，將邊界層厚度的修正轉換到極坐標系中進行；

（5）以HSP槳作為算例，計算槳葉r/R=0.7處，剖面在旋轉至 0°、90°、180°、270°時的非定常壓力系數，對該剖面在各角度處是否出現空泡或空泡位置作出推斷；計算了槳葉r/R=0.7處，剖面在0.25弦長、0.4弦長、0.6弦長、0.8弦長位置旋轉一周的壓力變化。

1 計算方法

1.1 Theodorsen法

Theodorsen法的基本思想是先把給定的槳葉剖面周線通過儒可夫斯基變換式轉換為近似圓，再利用特定的級數形式的變換函數將近似圓變換為準確圓，求得槳葉剖面繞流速度分布與槳葉剖面坐標的關系。由伯努利方程可知，在已知速度分布的情況下，可求得槳葉剖面各點的壓力分布。

在ζ平面上給定一個槳葉剖面，其周線為C，弦長l=4a。以弦線中點為原點建立坐標系，令前緣點A點坐標為（-2a，0），T 點坐標為（2a，0），見圖 2。

通常由升力線、升力面初步計算得出的槳葉剖面形狀均有一定的攻角α，并不適合Theodorsen法坐標系的建立。因此需對該槳葉剖面初始坐標（s，t）點進行坐標旋轉，轉換如下：

式中：

然后，通過下列儒可夫斯基變換式將槳葉剖面周線C變換為z′平面上的近似圓。

令z′平面上復數為：

式中：Ψ和ω是z′平面上復坐標變量。將式（4）代入式（3），建立變換的對應關系：

所以有

從式（6）可解得ω和ψ為

因此，根據槳葉剖面坐標（x，y），便可通過式（7）、式（8）確定近似圓上對應點處極坐標變量ω和ψ，并對離散點Ψ（ω）進行擬合、求導，最終得到。

通過對NACA66 a=0.8（Mod）的某槳葉剖面的計算，得到數據如圖3所示。

圖3 槳葉剖面各點在z′平面內的坐標分布

通過圖3對離散點分布的觀察，可將葉背、葉面曲線擬合的目標函數分別定義為：

對葉背曲線的擬合如圖4所示。

圖4 離散點曲線與擬合的函數間的關系

本文采用 Levenberg-Marquardt法［9］對目標函數各變量進行計算，通常其相關系數均可達到0.995以上。

下一步，利用級數形式函數式（10）與式（11），將z′平面上近似圓變換到z平面上為一準確圓。

對式（10）利用求三角級數系數的方法，可得常數

從而可得z平面上圓半徑

根據共軛三角函數關系知γ與ψ是共軛的（因χ0是常數）。 如果已知 ψ（θ），則 γ（θ）就可確定。 由泊松公式可知：

式中：γ0= （ω-θ）0表示 θ-θ0（特定值）時的 γ 值。 然而，我們知道的是ψ（ω），尚不知道ψ與θ之間的關系。但因為近似圓C1與圓S相差并不太大，可假定θ與 ω 之差不大，而令 ψ=ψ（ω）=ψ（θ）作為第一級近似，因此，將它代入式（14），可得到 γ=ω-θ的第一級近似的 γ0（1）。

求出 γ0（1）后，于是得到 ω=θ+γ0（1）。 將 ω=θ+γ0（1）代入 ψ=ψ（ω）=ψ（θ+γ（1））后，可得 γ 的第二級近似：

以此類推。但是，實際上第一級近似已有滿意的結果，然后求出 γ（2）=γ（2）（ω）的數值關系，以及的分布。

第四步，根據式（12）、（13）求得 χ0、r0。

第五步，根據Kutta條件，對應于z平面上圓柱繞流，應規定這一點（θ=-γ處）為后駐點，以保證槳葉剖面繞流在尾緣處的速度為有限值。

第六步，將上述求得的結果代入式（18）

由式（18）便可求得槳葉剖面上速度與ω之間的關系。式中α為來流攻角，可在合理范圍內任意給定。因式（7）、式（8）已確定 ψ、ω 與槳葉剖面上（x，y）的關系，所以沿槳葉剖面上速度分布也隨之求出。

1.2 邊界層修正

勢流理論是基于理想流體對速度分布、壓力分布進行求解，該解必須經過粘性修正才能與實驗結果基本相符。本文采用邊界層修正來近似替代粘性對槳葉剖面的影響。

假定使用Theodorsen法進行計算得到的速度勢流解作為邊界層外緣處的速度分布，引入邊界層位移厚度 δ1（x）和邊界層動量厚度 δ2（x）：

并定義動量厚度雷諾數

引入形狀因子

為便于使用，將邊界層動量積分方程寫為：

式中：Cf為摩擦系數

單獨由動量積分方程式（23）無法求解，因此必須補充方程。在計算層流邊界層時，本文使用Thwaites法來求解動量積分厚度和位移厚度；對于紊流邊界層的計算，則采用Head法所提出的補充方程進行求解，用修正后的米歇爾（Michel）經驗公式進行轉折點的判斷：

式中：RexT

為Michel基于實驗提出的轉折點位置XT處的雷諾數。

下面分別對層流和湍流邊界層進行計算。

1.2.1 使用Thwaites法對層流邊界層進行計算

最終便可求得位移厚度 δ1= δ2·H12。

1.2.2 使用Head法對紊流邊界層進行計算

Head法針對動量方程提出兩個補充方程：

（1）卷吸積分關系式：

由實驗數據確定其中變量

（2）摩擦系數關系式：

摩擦系數Cf與動量厚度雷諾數Re2有關，利用Ludwieg和Tillmann提出的關系式：

由式（23）、（29）、（30）、（31）可組成封閉方程，從而可求出 δ2、H12、δ1。

在求得ζ平面上的位移厚度δ1（x）后，將其轉化為 Z′平面內的位移厚度 δ1z′（ω），然后對原槳葉剖面進行粘性修正，其葉背、葉面處的修正分別為：

將物面外推δ1z′（ω），對包括排擠厚度的槳葉剖面再作勢流解，重復Theodorsen法各個步驟，便是粘性修正后的初次近似結果。將該近似結果繼續進行邊界層修正，反復迭代，直至迭代結果滿足精度要求。式（32）中n為迭代次數。得到滿意的速度分布后，可由伯努利方程推導出槳葉剖面壓力系數分布：

1.3 空泡校核

空化數計算公式：

由此，可比較槳葉剖面各點處的壓力系數與空泡數大小，從而判斷該點處是否會產生空泡。

2 結果分析

對HSP槳槳葉0.7半徑處剖面進行壓力分布計算。表1所示為HSP螺旋槳的主要參數；表2所示為該槳的試驗工況。

表1 HSP螺旋槳的主要參數

表2 HSP槳的試驗工況

表3所示為該槳的槳葉幾何參數，圖5是該槳的伴流分布。

表3 HSP槳葉型值

圖5 HSP的伴流分布

圖6為HSP槳0.7半徑剖面的非定場壓力系數（0°）。由圖6可推斷，槳葉轉到此角度時，槳葉剖面自接近導邊至x/C=0.45處壓力系數均大于空泡數，因此可推知此段槳葉會產生空泡。

圖6 HSP槳0.7半徑剖面的非定場壓力系數（0°）

圖7為HSP槳0.7半徑剖面的非定場壓力系數（90°）。由圖7可推斷，槳葉轉至此角度時，槳葉剖面各處壓力系數均小于空泡數，因此可推知不會產生空泡。

圖7 HSP槳0.7半徑剖面的非定場壓力系數（90°）

圖8為HSP槳0.7半徑剖面的非定場壓力系數（180°）。由圖8可推斷，槳葉轉至此角度時，槳葉剖面在x/C=0.1～0.35處葉背開始產生空泡。

圖8 HSP槳0.7半徑剖面的非定場壓力系數（180°）

圖9為HSP槳在0.7半徑剖面處的非定場壓力系數（270°）。由圖9可推斷，槳葉轉至此角度時，槳葉剖面不會產生空泡。

圖9 HSP槳在0.7半徑剖面的非定場壓力系數（270°）

HSP 槳在 0.7 半徑剖面、（0.25、0.4、0.6、0.8）弦長位置處旋轉一周的壓力變化如圖10～圖13所示。

圖10 HSP槳在0.7半徑剖面、0.25弦長位置旋轉一周的壓力變化（r/R=0.7，x/C=0.25）

圖11 HSP槳在0.7半徑剖面、0.4弦長位置旋轉一周的壓力變化 （r/R=0.7，x/C=0.4）

圖12 HSP槳在0.7半徑剖面、0.6弦長位置旋轉一周的壓力變化（r/R=0.7，x/C=0.6）

圖13 HSP槳在0.7半徑剖面、0.8弦長位置旋轉一周的壓力變化（r/R=0.7，x/C=0.8）

3 結 論

本文以第22屆ITTC推進器技術委員會在1998年發布的HSP螺旋槳為算例，計算槳葉r/R=0.7處剖面在旋轉至 0°、90°、180°、270°時的非定常壓力系數，以及 r/R=0.7 處剖面在（0.25、0.4、0.6、0.8）弦長位置旋轉一周的壓力變化，并判斷在各角度下是否會產生空泡以及空泡發生的位置。結果顯示，計算值與實驗值符合性較好，為螺旋槳設計者提供了一種可與螺旋槳升力線、升力面設計法相結合，初步判斷槳葉剖面的空泡位置和空泡形態的方法。

［1］TERRYB.Mininum pressure envelopes for modified NACA-66 section with NACA a=0.8 camber and busfips typeⅡsection［J］.Science.1965，148：12-64.

［2］黃汝道、陳才源.考慮來流速度場不均勻性影響的螺旋槳設計計算方法［J］.艦船性能研究，1982（4）：33-54.

［3］胡健.螺旋槳空泡性能及低噪聲螺旋槳設計研究［D］.哈爾濱：哈爾濱工程大學，2006.

［4］方會東.船舶螺旋槳空泡性能理論預報［D］.哈爾濱：哈爾濱工程大學，2008.

［5］王超.螺旋槳水動力性能、空泡及噪聲性能的數值預報研究［D］.哈爾濱：哈爾濱工程大學，2010.

［6］王國強，董世湯.船舶螺旋槳理論與應用［M］.哈爾濱：哈爾濱工程大學出版社，2005：235-313.

［7］王獻孚.船用翼理論［M］.北京：國防工業出版社，1998：50-57.

［8］劉岳元，馮鐵城，劉應中.水動力學基礎［M］.上海：上海交通大學出版社，1990：230-247.

［9］JEFFERY J L著.張威，劉志軍，李艷紅，等譯.數值分析與科學計算［M］.北京：清華大學出版社，2008：441-443.