基于Theodorsen方法的螺旋槳空泡數(shù)值計(jì)算研究

趙 鵬 周玉龍 趙鵬飛 楊 超

（江蘇科技大學(xué) 船舶與海洋工程 鎮(zhèn)江212003）

0 引 言

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

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

圖1 初步提供的槳葉剖面

圖2 進(jìn)行坐標(biāo)調(diào)整后的槳葉剖面

（1）整理出由Theodorsen法和邊界層修正來預(yù)報(bào)空泡的整體思路；

（2）對(duì)有攻角的槳葉剖面進(jìn)行坐標(biāo)轉(zhuǎn)換，使其滿足Theodorsen方法的要求；

（3）對(duì)近似圓中坐標(biāo)點(diǎn)所需進(jìn)行的線性擬合的目標(biāo)函數(shù)定義為三角函數(shù)；

（4）邊界層修正時(shí)，為便于計(jì)算，將邊界層厚度的修正轉(zhuǎn)換到極坐標(biāo)系中進(jìn)行；

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

1 計(jì)算方法

1.1 Theodorsen法

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

在ζ平面上給定一個(gè)槳葉剖面，其周線為C，弦長(zhǎng)l=4a。以弦線中點(diǎn)為原點(diǎn)建立坐標(biāo)系，令前緣點(diǎn)A點(diǎn)坐標(biāo)為（-2a，0），T 點(diǎn)坐標(biāo)為（2a，0），見圖 2。

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

式中：

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

令z′平面上復(fù)數(shù)為：

式中：Ψ和ω是z′平面上復(fù)坐標(biāo)變量。將式（4）代入式（3），建立變換的對(duì)應(yīng)關(guān)系：

所以有

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

因此，根據(jù)槳葉剖面坐標(biāo)（x，y），便可通過式（7）、式（8）確定近似圓上對(duì)應(yīng)點(diǎn)處極坐標(biāo)變量ω和ψ，并對(duì)離散點(diǎn)Ψ（ω）進(jìn)行擬合、求導(dǎo)，最終得到。

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

圖3 槳葉剖面各點(diǎn)在z′平面內(nèi)的坐標(biāo)分布

通過圖3對(duì)離散點(diǎn)分布的觀察，可將葉背、葉面曲線擬合的目標(biāo)函數(shù)分別定義為：

對(duì)葉背曲線的擬合如圖4所示。

圖4 離散點(diǎn)曲線與擬合的函數(shù)間的關(guān)系

本文采用 Levenberg-Marquardt法［9］對(duì)目標(biāo)函數(shù)各變量進(jìn)行計(jì)算，通常其相關(guān)系數(shù)均可達(dá)到0.995以上。

下一步，利用級(jí)數(shù)形式函數(shù)式（10）與式（11），將z′平面上近似圓變換到z平面上為一準(zhǔn)確圓。

對(duì)式（10）利用求三角級(jí)數(shù)系數(shù)的方法，可得常數(shù)

從而可得z平面上圓半徑

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

式中：γ0= （ω-θ）0表示 θ-θ0（特定值）時(shí)的 γ 值。 然而，我們知道的是ψ（ω），尚不知道ψ與θ之間的關(guān)系。但因?yàn)榻茍AC1與圓S相差并不太大，可假定θ與 ω 之差不大，而令 ψ=ψ（ω）=ψ（θ）作為第一級(jí)近似，因此，將它代入式（14），可得到 γ=ω-θ的第一級(jí)近似的 γ0（1）。

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

以此類推。但是，實(shí)際上第一級(jí)近似已有滿意的結(jié)果，然后求出 γ（2）=γ（2）（ω）的數(shù)值關(guān)系，以及的分布。

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

第五步，根據(jù)Kutta條件，對(duì)應(yīng)于z平面上圓柱繞流，應(yīng)規(guī)定這一點(diǎn)（θ=-γ處）為后駐點(diǎn)，以保證槳葉剖面繞流在尾緣處的速度為有限值。

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

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

1.2 邊界層修正

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

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

并定義動(dòng)量厚度雷諾數(shù)

引入形狀因子

為便于使用，將邊界層動(dòng)量積分方程寫為：

式中：Cf為摩擦系數(shù)

單獨(dú)由動(dòng)量積分方程式（23）無法求解，因此必須補(bǔ)充方程。在計(jì)算層流邊界層時(shí)，本文使用Thwaites法來求解動(dòng)量積分厚度和位移厚度；對(duì)于紊流邊界層的計(jì)算，則采用Head法所提出的補(bǔ)充方程進(jìn)行求解，用修正后的米歇爾（Michel）經(jīng)驗(yàn)公式進(jìn)行轉(zhuǎn)折點(diǎn)的判斷：

式中：RexT

為Michel基于實(shí)驗(yàn)提出的轉(zhuǎn)折點(diǎn)位置XT處的雷諾數(shù)。

下面分別對(duì)層流和湍流邊界層進(jìn)行計(jì)算。

1.2.1 使用Thwaites法對(duì)層流邊界層進(jìn)行計(jì)算

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

1.2.2 使用Head法對(duì)紊流邊界層進(jìn)行計(jì)算

Head法針對(duì)動(dòng)量方程提出兩個(gè)補(bǔ)充方程：

（1）卷吸積分關(guān)系式：

由實(shí)驗(yàn)數(shù)據(jù)確定其中變量

（2）摩擦系數(shù)關(guān)系式：

摩擦系數(shù)Cf與動(dòng)量厚度雷諾數(shù)Re2有關(guān)，利用Ludwieg和Tillmann提出的關(guān)系式：

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

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

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

1.3 空泡校核

空化數(shù)計(jì)算公式：

由此，可比較槳葉剖面各點(diǎn)處的壓力系數(shù)與空泡數(shù)大小，從而判斷該點(diǎn)處是否會(huì)產(chǎn)生空泡。

2 結(jié)果分析

對(duì)HSP槳槳葉0.7半徑處剖面進(jìn)行壓力分布計(jì)算。表1所示為HSP螺旋槳的主要參數(shù)；表2所示為該槳的試驗(yàn)工況。

表1 HSP螺旋槳的主要參數(shù)

表2 HSP槳的試驗(yàn)工況

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

表3 HSP槳葉型值

圖5 HSP的伴流分布

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

圖6 HSP槳0.7半徑剖面的非定場(chǎng)壓力系數(shù)（0°）

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

圖7 HSP槳0.7半徑剖面的非定場(chǎng)壓力系數(shù)（90°）

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

圖8 HSP槳0.7半徑剖面的非定場(chǎng)壓力系數(shù)（180°）

圖9為HSP槳在0.7半徑剖面處的非定場(chǎng)壓力系數(shù)（270°）。由圖9可推斷，槳葉轉(zhuǎn)至此角度時(shí)，槳葉剖面不會(huì)產(chǎn)生空泡。

圖9 HSP槳在0.7半徑剖面的非定場(chǎng)壓力系數(shù)（270°）

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

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

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

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

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

3 結(jié) 論

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

［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］黃汝道、陳才源.考慮來流速度場(chǎng)不均勻性影響的螺旋槳設(shè)計(jì)計(jì)算方法［J］.艦船性能研究，1982（4）：33-54.

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

［4］方會(huì)東.船舶螺旋槳空泡性能理論預(yù)報(bào)［D］.哈爾濱：哈爾濱工程大學(xué)，2008.

［5］王超.螺旋槳水動(dòng)力性能、空泡及噪聲性能的數(shù)值預(yù)報(bào)研究［D］.哈爾濱：哈爾濱工程大學(xué)，2010.

［6］王國(guó)強(qiáng)，董世湯.船舶螺旋槳理論與應(yīng)用［M］.哈爾濱：哈爾濱工程大學(xué)出版社，2005：235-313.

［7］王獻(xiàn)孚.船用翼理論［M］.北京：國(guó)防工業(yè)出版社，1998：50-57.

［8］劉岳元，馮鐵城，劉應(yīng)中.水動(dòng)力學(xué)基礎(chǔ)［M］.上海：上海交通大學(xué)出版社，1990：230-247.

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