周　衛(wèi)，張吉萍，許頌捷

（1.浙江海洋學院船舶與海洋工程學院，浙江舟山　316022；2.浙江省海洋漁業(yè)船舶交易服務(wù)中心，浙江杭州　310007；3.浙江省近海海洋工程技術(shù)重點實驗室，浙江舟山　316022；4.上海羽翼船舶設(shè)備有限公司，上?！?01600）

零航速減搖鰭升力特性數(shù)值研究

周衛(wèi)1，2，張吉萍1，3，許頌捷4

（1.浙江海洋學院船舶與海洋工程學院，浙江舟山316022；2.浙江省海洋漁業(yè)船舶交易服務(wù)中心，浙江杭州310007；3.浙江省近海海洋工程技術(shù)重點實驗室，浙江舟山316022；4.上海羽翼船舶設(shè)備有限公司，上海201600）

為了探究零航速減搖鰭周期擺動時產(chǎn)生的升力與轉(zhuǎn)動頻率、轉(zhuǎn)角幅度、轉(zhuǎn)軸位置以及來流速度幾方面的關(guān)系，本文運用滑移網(wǎng)格技術(shù)對二維NACA0012翼型的周期擺動進行數(shù)值預(yù)報。研究表明：Transition SST湍流模型較其他常見湍流模型更適合分析大攻角周期擺動問題；對于轉(zhuǎn)動方式滿足余弦函數(shù)的減搖鰭，產(chǎn)生的升力隨轉(zhuǎn)角幅度、轉(zhuǎn)動頻率增加而增大，來流速度對升力的影響大于轉(zhuǎn)軸位置。論文研究對零航速減搖鰭選型設(shè)計具有一定指導(dǎo)意義。

零航速；減搖鰭；升力特性；Transition SST；NACA0012

隨著人們對低航速下的舒適性要求的提高，零航速減搖概念自提出便受到國內(nèi)外機構(gòu)及學者的關(guān)注，其工作原理是快速擺動鰭體，產(chǎn)生抵抗船舶橫搖力矩，實現(xiàn)低航速或停泊狀態(tài)的減搖。目前關(guān)于零航速減搖鰭的研究主要以試驗研究為主，美國及荷蘭相關(guān)機構(gòu)在1998年開展了探索性試驗研究［1-2］，昆騰公司利用試驗方法開展全航速減搖鰭的初步設(shè)計［3］。在國內(nèi)，金鴻章等［4］根據(jù)實驗數(shù)據(jù)及勢流理論構(gòu)建零航速減搖鰭的力學模型，表明升力受旋轉(zhuǎn)角速度、角加速度等因素的影響。

零航速減搖鰭的周期性擺動實際是以旋渦產(chǎn)生、耗散為主導(dǎo)的湍流運動，對于這種復(fù)雜的湍流運動，試驗研究受試驗器材和設(shè)備的限制，往往僅能提供有限的積分變量，而CFD數(shù)值方法能提供更完善的流場信息。同時伴隨湍流模型的完善，及計算精度與穩(wěn)定性的提高，運用數(shù)值方法進行零航速減搖鰭的研究具有較明顯的優(yōu)勢。

本文以二維NACA0012翼型作為研究對象，運用Fluent軟件自帶的UDF函數(shù)編譯各種運動方式，并結(jié)合滑移網(wǎng)格技術(shù)進行數(shù)值預(yù)報。針對試驗工況，對各類湍流模型、網(wǎng)格大小等相關(guān)設(shè)置的影響進行分析，表明數(shù)值方法能較好的預(yù)報大攻角動態(tài)失速問題。其次針對減搖鰭運動特征，分析了轉(zhuǎn)角幅度、轉(zhuǎn)動周期、轉(zhuǎn)軸位置、來流速度幾方面對于升力的影響。

1　數(shù)值方法

1.1湍流模型

直接計算法和雷諾平均法（RANS）是進行數(shù)值分析的兩種主要方法，RANS由于存在表示湍流脈動對平均流動影響的雷諾應(yīng)力項τij，因此需要添加湍流模型進行封閉。常用的k-ω和k-ε湍流模型，由于存在較大耗散會制約旋渦的產(chǎn)生，無法準確的預(yù)報旋渦的產(chǎn)生的位置、旋渦擴散的過程［5］。因此，本文主要采用Transition SST模型進行分析減搖鰭的周期擺動問題。

Transition SST湍流模型是基于k-ω SST模型，同時考慮湍流轉(zhuǎn)捩問題，因此運輸方程還需構(gòu)造關(guān)于動能厚度雷諾數(shù)Reθ和間歇因子γ的運輸方程，并通過控制方程與SST k-ω模型結(jié)合起來。

式（1）、（2）為SST k-ω湍流模型，式（3）、（4）為間歇因子γ及轉(zhuǎn)捩雷諾數(shù)的運輸方程；式（1）、（2）中：ρ為流體密度，μj對應(yīng)速度矢量，xj對應(yīng)位置矢量，Gk和Gω為關(guān)于k和ω的產(chǎn)生項，Гω、Гk為有效擴散項；Yω、Yk為各自的湍流耗散項，而Sω、Sk則表示為可自定義的源項。式（3）中Pγ1、Eγ1為過渡源項，Pγ2、Eγ2為破壞源項。

轉(zhuǎn)捩模型與SST湍流模型結(jié)合的控制方程為：

G%k與Y*k為SST模型中的初始的生成項與耗散項。

1.2研究對象

目前，關(guān)于零航速減搖鰭的試驗數(shù)據(jù)有限，因此文本選用LEE等［6］針對大攻角失速研究采用的NACA0012翼型作為研究對象。選用該試驗數(shù)據(jù)進行驗證分析主要有兩方面原因：首先是該實驗涉及大攻角翼型的雷諾數(shù)范圍與減搖鰭擺動過程的雷諾數(shù)范圍相近，同時空氣與水的流動特性相近。其次，該試驗得到的各項系數(shù)已剔除三維流動對是實驗結(jié)果的影響，因此試驗數(shù)據(jù)適用于二維的分析。NACA0012翼型的擺動原理具體設(shè)置可參見文獻［6］。表1為數(shù)值計算涉及工況的主要參數(shù)。

表1　工況參數(shù)Tab.1　Parameters of Different Loadcases

1.3網(wǎng)格及計算設(shè)置

文中所計算的兩種工況采用相同的網(wǎng)格進行計算，網(wǎng)格利用ICEM軟件的O型網(wǎng)格形式進行劃分，翼型表面初始網(wǎng)格高度為0.02mm，滿足y+值小于1要求，網(wǎng)格沿法線方向的高度以1.2比率增長，確保能捕獲近壁面附近的小尺度脈動同時降低網(wǎng)格數(shù)量，最終確定的初始網(wǎng)格數(shù)量為12萬。在此基礎(chǔ)上，弦長方向上的網(wǎng)格加密1.5倍，分別生成18萬和27萬網(wǎng)格。

圖1　NACA翼型表面網(wǎng)格Fig.1　Surface Mesh of NACA Airfoils

圖2　計算域及邊界條件Fig.2　Calculation Domain and boundary conditions

數(shù)值計算采用的邊界條件與試驗條件相同，具體設(shè)置如下：入口位置：翼前18倍弦長，設(shè)置為速度入口；遠場條件：距翼型18倍弦長，設(shè)置為速度入口；出口位置：翼后24倍弦長位置，設(shè)置為壓力出口；Interface邊界：距翼型10倍弦長。

本文涉及的計算均采用雙精度基于壓力法求解器進行求解，對流項采用二階迎風差分格式，壓力速度差分方程組的偶合求解采用SIMPLE格式，采用中心差分格式進行離散。對于DDES湍流模型，采用有界中心差分格式。運用滑移網(wǎng)格技術(shù)以及UDF函數(shù)實現(xiàn)不同運動方式的控制。

2　計算結(jié)果與分析

2.1驗證工況分析

圖3為時間步長是0.000 2 s下，不同湍流模型計算得到升力系數(shù)。對于所有湍流模型對于從-5°到25°的上揚階段準確性均優(yōu)于25°到-5°的下擺階段，無論是升力變化趨勢還是與試驗結(jié)果的差值。

上擺階段：角度滿足-5°≤α≤25°范圍時：RNG k-ε、DDES、Transition SST三種湍流模型均做出較好的預(yù)報，在20°≤α≤25°范圍，三種湍流模型表現(xiàn)較大差異，DDES模型預(yù)報的失速角度為21°，且升力系數(shù)峰值較實驗值偏小15.6%；RNG k-ε模型則升力系數(shù)隨轉(zhuǎn)角增大而保持震蕩增加，并未出現(xiàn)明顯的失速現(xiàn)象，這是由于RNG k-ε模型無法模擬前緣渦分離現(xiàn)象。Transition SST模型預(yù)測的失速角度為24.5°，25°時的升力系數(shù)下降約0.9，表明Transition SST模型過度預(yù)報分離現(xiàn)象，將直接影響下擺過程的準確性。

下擺階段：采用加強避免函數(shù)處理的RNG k-ε得到的升力系數(shù)變化較為平穩(wěn)，與試驗結(jié)果吻合較好；Transition SST模型和DDES模型在15°≤α≤25°范圍，升力系數(shù)波蕩較大，這是由于過度預(yù)報分離現(xiàn)象，造成的二次分離預(yù)報不準確。在15°≤α≤25°，Transition SST模型預(yù)報的升力系數(shù)的震蕩現(xiàn)象減小明顯，趨勢與實驗結(jié)果吻合。

圖4為12萬、18萬、27萬網(wǎng)格采用Transition SST模型、時間步長采用0.0002 s時計算得到的升力系數(shù)。上揚階段：三種網(wǎng)格得到的升力系數(shù)幾乎貼合，在20°至25°過程，12萬網(wǎng)格與其他網(wǎng)格得到升力系數(shù)數(shù)值差異較大；下擺階段：三種網(wǎng)格得到升力系數(shù)變化趨勢整體保持一致，27萬網(wǎng)格與試驗結(jié)果的數(shù)值差異最小。整體而言，18萬網(wǎng)格策略為較優(yōu)的網(wǎng)格策略，在后續(xù)的計算分析中采用18萬網(wǎng)格進行研究分析。

圖3　不同湍流模型的升力系數(shù)Fig.3　Lift Parameters of Different Turbulence Models

圖4　不同網(wǎng)格劃分的升力系數(shù)Fig.4　Lift Parameters of Different Mesh

2.2零航速數(shù)值分析

零航速減搖鰭的運動方式主要有兩種，一是采用加速-勻速-減速的運動方式［7］，二是滿足余弦運動規(guī)律，文中采用第二種運動方式進行分析，并對轉(zhuǎn)角幅度、轉(zhuǎn)動周期對升力影響進行分析。

圖5　不同轉(zhuǎn)角幅度的升力變化曲線Fig.5　Lift Curve of different Angle Ranges

圖6　不同轉(zhuǎn)角周期的升力變化曲線Fig.6　Lift Curve of different Angle Periods

圖5為零航速減搖鰭的轉(zhuǎn)動周期為1 s，轉(zhuǎn)角幅度分別為35°、40°、45°、50°時的升力變化曲線，由圖5可以看到，四種轉(zhuǎn)角幅度的升力曲線在上揚或下擺階段升力峰值均經(jīng)歷先增大后減小的過程。下擺階段除轉(zhuǎn)角幅度為50°時升力峰值位于11°附近，其他轉(zhuǎn)角幅度的升力峰值出現(xiàn)在0°（最大角速度位置）附近；上揚階段所有轉(zhuǎn)角幅度的升力峰值均出現(xiàn)在12°附近，即超前于最大角速度所在位置，這是由于附連水質(zhì)量力的影響。此外，轉(zhuǎn)角幅度增加5°，下擺階段升力峰值平均增幅約35%，上揚階段升力峰值平均增幅約22%。

圖6對應(yīng)的是轉(zhuǎn)角幅度為40°，轉(zhuǎn)動周期分別為0.5 s、1 s、2 s時的升力變化曲線。通過比較可以發(fā)現(xiàn)，三種周期下的升力峰值出現(xiàn)的位置基本一致，升力峰值與轉(zhuǎn)動周期的平方成反比，即升力峰值與轉(zhuǎn)動角速度平方成正比。

圖7　不同轉(zhuǎn)軸位置的升力變化曲線Fig.7　Lift Curve of different Rotor Positions

圖8　不同流速下的升力變化曲線Fig.8　Lift Curve of different Velocitys

圖7為四種轉(zhuǎn)軸位置的升力變化曲線，鰭軸位置向前移動，升力變化趨勢基本保持不變，上擺過程的升力峰值增大，下擺過程的升力峰值減小。由圖8可以看到，流速達到0.2m/s時的升力變化趨勢和升力峰值均與航速低于0.1m/s時發(fā)生明顯變化，上擺和下擺過程升力變化趨勢相似、升力峰值相近，表明零、低航速時減搖鰭的升力變化受航速影響較大。

3　結(jié)論

本文利用數(shù)值方法對零航速減搖鰭上下擺動過程中涉及的大攻角動態(tài)失速問題進行研究，通過上述的分析可以得到以下結(jié)論：

（1）對于大攻角動態(tài)失速問題，Transition SST湍流模型的預(yù)報結(jié)果優(yōu)于常用的RNG k-ε湍流模型，DDES湍流模型存在過度預(yù)報旋渦分離現(xiàn)象。

（2）運動規(guī)律滿足α=α1cos（2πf·t）的零航速減搖鰭，在轉(zhuǎn)動過程中減搖鰭產(chǎn)生的升力與轉(zhuǎn)角幅度成正比、同時與轉(zhuǎn)動周期的平方成反比。改變鰭軸位置對于升力的變化影響不大，來流速度對于減搖鰭的性能影響明顯。

（3）對于零航速減搖鰭設(shè)計初期，可利用數(shù)值方法輔助開展翼型評估、轉(zhuǎn)軸位置分析、轉(zhuǎn)動形式優(yōu)化等方面工作。

［1］DALLINGA R P.Roll stabilization of motor yacht：Use of fin stabilizers in anchored conditions［Z］.Project 99，Amsterdam，1999.

［2］DALLINGA R P.Roll stabilization of motor yacht：Hydrodynamic aspects of the comparison of anti-roll tanks and fins［Z］. Project 2002，Amsterdam，2002.

［3］OOMS J.The use of roll stabilizer fins at zero speed［Z］.Project 2002，Amsterdam，2002.

［4］金鴻章，張曉飛.零航速減搖鰭升力模型研究［J］.海洋工程，2007，25（3）：83-87.

［5］ANSYS.ANSYS FLUENT 12.0 Theory Guide［Z］.Ansys Inc，USA，2009.

［6］LEE T，GERONTAKOS P.Investigation of flow over an oscillating airfoil［J］.Journal of Fluid Mechanics，2004，512：313-341.

［7］宋吉廣.基于升力反饋的全航速減搖鰭研究［D］.哈爾濱：哈爾濱工程大學，2010.

Research on Lift Characteristics of Zero-speed Fin Stabilizer

ZHOU Wei1，2，ZHANG Ji-ping1，3，XU Song-jie4
（1.School of Naval Architecture and Ocean Engineering，Zhejiang Ocean University，Zhoushan316022；2.Zhejiang Province Marine Fishing Vessel Trading Service Center，Hangzhou310007；3.Key Laboratory of Offshore Engineering Technology of Zhejiang Provine，Zhoushan316022；4.Yoyeah Marine，Shanghai201600，China）

To explore the influence from rotational frequency，angle range，rotor position and flow velocity to the cyclic swing lift of Zero-speed Fin Stabilizer，the numerical forecast to the cyclic swing of two-dimensional NACA0012 airfoils is carried out with the sliding mesh technique.The following conclusions can be obtained：Transition SST turbulence model is more appropriate to analysis the cyclic swing at high angles of attack that other turbulence models；to the fin stabilzer with the rotating mode meeting the cosine function，the lift increased with the increasing of angle range and rotational frequency，and the influence from the flow velocity to the lift is more than the influence from the rotor position.The research is of guiding significance to the type selection and design of zero-speed fin stabilizer.

zero-speed；fin Stabilizer；lift characteristics；Transition SST；NACA0012

U664.72

A

1008-830X（2016）02-0150-05

2015-11-23

浙江省自然科學基金（LY14E090002）；國家自然科學基金（51409232）

周衛(wèi)（1987-），男，碩士研究生，研究方向：農(nóng)業(yè)機械化.