NACA 0012擺動水翼水動力特性的二維數值模擬

呂元博 ，田新亮 ，李欣 ，宋春輝

1上海交通大學海洋工程國家重點實驗室，上海200240

2高新船舶與深海開發裝備協同創新中心，上海200240

0 引言

自Hine等［1］率先研制出波浪滑翔機以來，作為一種直接將波浪能轉化為前進動能的海洋裝置，波浪滑翔機被廣泛應用于海上長時間觀測作業。由于波浪滑翔機具有續航能力強、噪聲小等優點［2-3］，其在檢測海洋環境要素、海洋災害預報、海洋科學研究等方面呈現出廣闊的應用前景。

波浪滑翔機的水翼是影響其航行性能的重要因素。針對水翼的水動力性能，國內外學者進行了大量研究。Kraus［4］通過建立模型，對波浪滑翔機各部位進行了六自由度模擬，確定水翼最大擺角為20°時最優，模擬得出滑翔機在不同海況下的各項水動力學參數。賈立娟［5］利用FLUENT軟件研究了翼型、擺角等對水翼水動力學特性的影響，發現擺角為18°時水翼的水動力性能較好。Sun等［6］針對果蠅翅膀的主動擺動進行了數值模擬，發現動態平均升力系數可以達到準靜態升力系數的2倍。Andersen等［7］對某一對稱翼型進行了主動擺動和主動垂蕩運動的數值與實驗研究，得出2種運動模型尾流場相似，但對于垂蕩運動在低頻高幅值時，在1個振蕩周期內將產生2對對稱的尾窩。胡合文［8］通過勢流與粘流結合的方法對滑翔機整體及水翼部分進行了計算，得到水翼的最佳擺角在 15°附近。張曉慶等［9］和劉煥興等［10］以魚類尾鰭擺動推進為背景，利用計算流體動力學（Computational Fluid Dynamics，CFD）軟件，模擬了二維剛性水翼擺動的非定常水動力性能。

本文選用NACA 0012作為波浪滑翔機的翼型。針對該翼型，Ohtake等［11］和 Yonemoto等［12］均進行過大量研究，獲得了定常流場下機翼的性能參數。Read等［13］和 Triantafyllou等［14］通過一套水翼主動擺動裝置進行實驗，分析了NACA 0012型水翼在主動擺動和升沉運動下的水動力性能，得出有效攻角的不平滑變化會使水翼性能降低。Schouveiler等［15］在麻省理工學院（MIT）拖曳水池通過實驗給定NACA 0012翼型正弦規律的升沉運動和主動擺動，得出最佳推進效率的斯特哈爾數St值為0.25～0.35，St較大時水動力性能會降低。由于波浪滑翔機是通過垂蕩方向的運動引起機翼的被動旋轉，因此，通過水翼的主動運動模擬機翼的實際運動可能存在較大誤差。被動旋轉則需要考慮水翼的擺動運動與流場運動的瞬態耦合，且要考慮水翼在流場中運動時所受力矩以及限位角的約束條件，存在一定的仿真難度，但對水翼擺動過程有更好的模擬，而目前針對非定常流場中的水翼被動擺動的研究少有涉及。

本文將基于STAR-CCM+軟件中的DFBI模塊，以NACA 0012翼型為研究對象，對波浪滑翔機水翼做垂蕩運動時引起的繞自身中心軸的被動旋轉運動進行模擬，對水翼主動旋轉與被動旋轉2種情況下的推力系數進行對比，并分析限位角、波浪參數等因素對水翼性能的影響，以期為波浪滑翔機的設計與研究提供參考。

1 理論概述

1.1 波浪滑翔機運動原理

波浪滑翔機的運動原理如圖 1（a）［16］所示。波浪滑翔機由水面母船和水下滑翔機2部分組成。當水面母船隨波浪上升時，纜繩會拉緊，帶動水下滑翔機向上運動，由于水的作用力，會使水下滑翔機的水翼逆時針翻轉。當水面母船隨波浪下降時，纜繩松弛，水下滑翔機自由下降，水翼受到水的反作用力而順時針擺動。水翼在不斷的上、下擺動過程中會產生向前的推力，拖拽水面母船不斷前進。

圖1（b）所示為波浪滑翔機水翼在波浪起伏過程中的受力分析結果。水翼在隨波浪上升的過程中，受垂直于翼面向下的水動力作用，會產生水平方向的分力Fx和垂直向下的分力Fy。水翼在隨波浪下降的過程中，受垂直于翼面向上的水動力作用，同樣會分解為水平方向的分力Fx和垂直向上的分力Fy。兩個過程中的水平分力Fx成為使波浪滑翔機向前運動的推進力。在此，將整個研究對象簡化為圖1（b）中的水翼模型。

1.2 水翼被動運動方程

波浪滑翔機水翼的運動分為2部分：一部分是在波浪驅動下的y方向的平移運動；另一部分是在平移過程中受水動力作用的擺角運動。

水翼的垂向運動h(t)為

式中：t為運動時間；f為運動頻率；h0為垂向運動幅值。

水翼擺角的被動運動控制方程為

式中：J為繞重心的轉動慣量；?為水翼被動旋轉的角加速度；M為作用在水翼上的水動力力矩，可以根據水翼表面壓強積分得到。對式（2）進行積分，可以得到角速度?和擺角θ。

在設計波浪滑翔機的水翼時，當水翼擺角θ達到限位角θ0時，擺角達到最大值，直到力矩反向，水翼開始反向擺動。即

1.3 水翼主動運動方程

對于水翼擺動問題，如果在垂直方向和擺動方向均采用主動運動模型，則其擺動方向的運動規律為

式中，ψ為平移運動和擺動運動的相位差。

對于主動旋轉工況，有效攻角α(t)為水平來流速度和平移速度的合速度與水翼的擺角之間形成的實際攻角，其表達式為

式中：h'(t)為平移速度；U為來流速度，即水翼前進速度。

在此，引入流體力學中的相似準則斯特哈爾數St：

式中，A為尾渦區域在平移方向的寬度，近似表示為平移幅度的2倍，即A=2h0。代入式（6），得St=2fh0/U。與式（5）聯立，可得

式中，αmax為最大有效攻角。

1.4 推力系數

水翼的瞬時推力系數為

式中：Fx(t)為水翼的瞬時推力；ρ為水的密度；b為水翼弦長；s為水翼展長。

式中，T為水翼旋轉周期。

平均推力系數為

2 數值模型及驗證

2.1 模型參數

水翼采用NACA 0012翼型剖面，最大厚度位于弦長1/3處，如圖2所示。本文取水翼特征弦長b=0.1 m，擺動軸位于最大厚度處。

2.2 數值模型

如圖3所示，二維擺動水翼的計算區域為50b×30b。翼表面為圓形計算域，其直徑為5b，用于翼的旋轉運動。對于邊界條件的設定，左邊界和上、下邊界均為速度入口，右邊界為壓力出口。

計算采用SST k-ω湍流模型，利用STAR-CCM+軟件的切割體網格對流域進行分割。如圖4所示，在水翼表面，采用了網格無相對變形且數值交換較好的重疊網格，并在計算域及尾流區域進行局部加密。邊界層取為30層，控制y+嚴格小于1。在水翼運動方面，采用STAR-CCM+軟件中的流體與固體相互作用（Dynamic Fluid Body Interaction，DFBI）模塊，通過輸入上、下表面邊界的正弦邊界條件，模擬波浪起伏過程，然后通過DFBI模塊計算出水翼在升沉運動下的被動旋轉運動情況，并通過單自由度旋轉模塊控制水翼可旋轉的最大角度，即限位角θ0。

2.3 數值模型驗證

為驗證所建立的數值模型是否滿足計算要求，對網格尺寸、時間步長、計算結果可靠性進行驗證。

2.3.1 網格尺寸驗證

考慮到計算的準確性及高效性，計算分析了4種網格尺寸固定角度下的升力系數CL，并進行比較。表1為在固定角度5°時4種網格尺寸A1～A4對應的CL值及相對誤差（指計算值與實驗值的相對誤差），從中可觀察到從A2開始計算值趨于穩定。根據對翼型的經驗分析，CL實驗值一般略小于CFD計算值，且可觀察到計算值與理論值十分接近。

表1 網格尺寸驗證Table 1 The mesh size validation

圖5為4種網格尺寸的水翼上、下表面壓力系數Cp分布曲線。從中可觀察到：A2，A3和A4均重合較好；A1在上表面約1/3處出現了波動，結合該處的流場進行分析，發現此處的邊界層出現了明顯分離；A1邊界層網格尺寸較大，無法精確捕捉。綜合考慮計算量及計算準確性，選取A2作為網格模型。

2.3.2 時間步長驗證

在A2網格模型下，通過STAR-CCM+軟件中的Motion模塊，對水翼施加主動升沉運動疊加繞旋轉中心的旋轉運動，時間步長Δt分別取T0/100，T0/200，T0/300，和T0/400共4種工況，計算升力系數峰值CLmax（表2）。由表2可知，除T0/100時間步長以外，其他工況的CLmax均較接近。

表2 網格尺寸A2時間步長驗證Table 2 Time step validation of A2

圖6為不同時間步長下，升力系數CL與時間T的關系曲線。與其他時間步長相比，T0/100在波峰波動處出現了較大的相位偏差，且波峰數量減少。綜合考慮計算量及計算準確性，最終選用Δt/T0=1/200為本計算模型的時間步長。

2.3.3 可靠性驗證

為驗證模型的準確性，在固定擺角、相同雷諾數下，將計算值與文獻［12］中的2個實驗值進行對比，攻角α分別在-5°～5°之間每隔1°選取，Re=105。圖7為固定擺角下水翼的升力及阻力系數。由圖7（a）可見計算值與實驗值吻合較好，且CL比Yonemoto的計算值更精確、數據點范圍更廣，圖中的直線為斜率a=2π的理論值。而從圖7（b）可知，阻力系數CD和Yonemoto計算值均出現了略低于實驗值的情況，且計算結果基本一致。因此，本文數值模型可準確用于計算水翼的水動力特性問題。

3 計算結果與分析

針對水翼主動與被動旋轉2種工況，在STAR-CCM+軟件中分別對水動力性能進行仿真，并對比仿真結果，著重討論限位角、波高、波浪頻率等對水翼推力系數的影響。主動與被動旋轉工況分別如表3和表4所示。

表3 水翼主動擺動計算工況Table 3 The calculation conditions for active oscillating wing

表4 水翼被動擺動計算工況Table 4 The calculation conditions for passive oscillating wing

3.1 水翼主動旋轉與被動旋轉對推進性能的影響

為了比較水翼主動旋轉與被動旋轉對推進性能的影響，選取進流速度U=0.4 m/s，最大有效攻角αmax=15°，h0=0.8 m，對不同St數下的主動與被動工況進行分析，結果如圖8所示。

從圖8可以發現，不論是主動擺動運動還是被動擺動運動，在一個周期內水翼的平均推力系數均隨St數的增大而增大，且主動擺動的平均推力系數大于被動擺動，采取主動擺動方法時其平均推力系數將比被動擺動高出約30%。分析原因，主要在于2種方式水翼擺動的旋轉速率不同。采取被動擺動方式模擬時，水翼在水中受到水反作用力的力矩，在升沉時瞬間完成翻轉，旋轉速率在瞬間達到最大，達到最大限位角后變為0，之后以此限位角做平移運動。而傳統的主動擺動方法給定水翼的正弦擺動速率，水翼在靠近平衡位置處旋轉速率最快，從平衡位置到波峰的過程中水翼并非瞬間完成翻轉，而是以正弦速率周期性擺動，所以此正弦速率的擺動過程會對水產生向后的額外推力，從而使水翼獲取額外的向前推進力，使平均推力系數增大。從能量的角度分析，被動擺動水翼的輸入能量僅僅為垂蕩運動所需能量；而主動擺動的輸入能量則為主動垂蕩運動與主動旋轉所需的能量之和，相較于被動擺動，增加了主動旋轉部分的能量，故也將產生更大的推力。

由圖8可以發現，在St數較小時，兩者的平均推力系數差值較小，隨著St數的逐漸增大，兩者的誤差逐漸增加并趨于穩定，二者平均推力系數的差別可達30%。可見，主動旋轉與被動旋轉對推進性能的影響較大，對于波浪滑翔機水翼的水動力性能研究，采用傳統的主動擺動運動模型會增加計算誤差，不能準確模擬波浪滑翔機水翼的實際運動狀態。

3.2 水翼擺動限位角θ對推進性能的影響

為了研究水翼擺動限位角對推進性能的影響，選取5種限位角進行分析，其他參數設置如下：進流速度U=0.4 m/s，波浪頻率f=0.2 Hz，波高h=0.8 m，推力系數時歷如圖9所示。

從圖9（a）可以分析得到，當水翼被動擺動的最大限位角θ0在 10°，15°和 20°時，推力系數隨時間做周期性的類正弦變化，且隨著最大限位角增加 C（t）也逐漸增加，在限位角θ0=20°時，推力系數峰值達到最大。

圖9（b）為水翼被動擺動的最大限位角θ0在20°，25°和30°工況下對應的C（t）-T曲線。由圖可知，在前半個周期，推力系數的峰值仍隨限位角的增加而增大，但在后半周期，在限位角θ0=25°和θ0=30°時 C（t）峰值出現了突變和尖點，沒有呈現很好的周期性。分析突變產生的原因，當擺角大于25°時，水翼開始逐漸發生失速，會受遠處來流與水翼上表面流動分離的相互影響，對推力性能產生較大影響。因此當擺角超過某一臨界值時，水翼推力系數將不再繼續升高，會出現明顯的波動。

圖10為在不同限位角下水翼在T/4時刻被動擺動運動時的渦量圖。由渦量圖可知，在限位角θ0=10°～20°時，水翼上表面分離區較小，呈現出流線型的繞流。但隨著擺角逐漸增大，從25°開始，水翼的尾緣產生明顯的渦旋，翼型上表面的分離區面積逐漸增加，后緣附近開始有翼尖渦，一定程度上降低了水翼的水動力性能。如圖10所示，水翼的脫落渦均呈逗號形式，尾部產生的射流增加了水翼的推力，隨著限位角的增大，脫落渦逐漸減小。

綜合以上分析，選取θ0=20°作為水翼被動擺動的最大限位角。

3.3 波高h和波浪頻率f對水翼推進性能的影響

確定限位角θ=20°后，在來流速度U=0.4 m/s，波浪頻率f=0.2 Hz的條件下，選取波高h=0.1，0.4和0.8 m共3種工況，得到了對應的推力系數C（t）關于周期T的曲線（圖11）。

由圖11可見，C（t）隨波高h的增大而增大，且在波高h=0.1和0.4 m時，推力系數曲線未呈現出明顯的周期性，C（t）的大部分小于0，說明流場并沒有趨于穩定。當h=0.8 m時，曲線開始出現明顯的周期性變化。此外，在波浪頻率f=0.2 Hz時，隨著波高h的增加，推力系數逐漸增大。且在一定來流速度和波浪頻率下，波高需達到一定的高度才可產生穩定向前的推力，流場和推力系數也可趨于穩定。

上述對于波高h的討論限定在單一頻率，即f=0.2 Hz條件下，規律不具有普遍性。為找到波高h對水翼推進性能的影響規律，分別選取多組工況進行了計算對比。

由于波高h=0.8 m時推力系數周期化變化，故取波高h在0.8 m附近變化，波浪頻率分別取f=0.1，0.2和0.4 Hz共3種情況，研究不同波高h對C（t）的影響規律（圖12）。由圖12可見，在不同波浪頻率下，推力系數C（t）均隨波高h的升高而增大。

由上述計算結果可知，在其他條件一定時，從波高h=0.8 m開始，推力系數C（t）趨于穩定。故在波高h=0.8 m，來流速度U=0.4 m/s時，選取不同波浪頻率進行計算。

圖13是h=0.8 m時不同波浪頻率下被動擺動水翼的C（t）-T曲線。由圖可知，當波浪頻率f=0.2 Hz時，推力系數在C（t）=0附近波動，波浪頻率f每增加一倍，C（t）增幅遠高于一倍。當f=0.8 Hz時，C（t）峰值達到20。可見，波浪頻率f對推力系數C（t）的影響更為明顯。在波高h=0.8 m時，推力系數C（t）隨波浪頻率f的增大而增加。

同樣，上述對于波浪頻率的討論也限定在單一波高，即h=0.8 m條件下，不具有普遍性，為找到波浪頻率f對水翼推進性能的影響規律，分別取多組工況進行了計算對比。

圖14所示為不同波浪頻率f對C（t）的影響規律。由圖可知，在不同波高下，推力系數C（t）隨f增加而增大，當波浪頻率從f=0.2 Hz擴大一倍至f=0.4 Hz時，C（t）曲線的峰值增加了 5～10倍。在波浪頻率較低時，即f=0.1 Hz和f=0.2 Hz時，推力系數C（t）峰值相差不大。由此可知，在波浪頻率較小時，波浪頻率f的改變對推力系數的影響較小。

4 結 論

本文以波浪滑翔機的水下滑翔部分水翼為研究對象，對NACA 0012型水翼在波浪中的被動擺動情況進行了水動力學分析及計算，得到如下主要結論：

1）通過自主建立的STAR-CCM+軟件的重疊網格與DFBI模塊相結合的模型，實現了對水翼垂直方向主動平移狀態下的被動旋轉的模擬。

2）相同條件下，主動擺動和被動擺動運動對水翼的推進性能有較大影響，二者的平均推力系數可相差30%。

3）通過分析計算結果與渦量圖，得到滑翔機水翼的最大限位角在20°附近時滑翔機的水動力性能達到最優。

4）在其他條件一定時，滑翔機水翼的推進性能受波浪參數影響較大，且隨波高h和波浪頻率f的增加而增加，在波浪頻率較小時，波浪頻率f的改變對推力系數影響不大。

參考文獻：

［1］HINE R，WILLCOX S，HINE G，et al.The wave glid?er：a wave-powered autonomous marine vehicle［C］//OCEANS 2009，MTS/IEEE Biloxi-Marine Technology for Our Future：Global and Local Challenges.Missis?sippi：IEEE，2009：1-6.

［2］徐春鶯，陳家旺，鄭炳煥.波浪驅動的水面波力滑翔機研究現狀及應用［J］.海洋技術學報，2014，33（2）：111-117.XU C Y，CHEN J W，ZHENG B H.Research status and applications of wave gliders［J］.Journal of Ocean Technology，2014，33（2）：111-117（in Chinese）.

［3］楊海，劉雁集，張凱.實驗尺度下無人水下滑翔機設計與試驗［J］.中國艦船研究，2016，11（1）：102-107，120.YANG H，LIU Y J，ZHANG K.Design and experi?ment for laboratory-scale autonomous underwater glid?ers［J］.Chinese Journal of Ship Research，2016，11（1）：102-107，120（in Chinese）.

［4］KRAUS N D.Wave glider dynamic modeling parame?ter identification and simulation［D］.Hawaii：Universi?ty of Hawaii，2012.

［5］賈立娟.波浪動力滑翔機雙體結構工作機理與動力學行為研究［D］.天津：國家海洋技術中心，2014.JIA L J.Study of operation principal of two-part archi?tecture and dynamic behavior of wave glider［D］.Tian?jin：National Ocean Technology Center，2014（in Chi?nese）.

［6］SUN M，TANG J.Unsteady aerodynamic force genera?tion by a model fruit fly wing in flapping motion［J］.Journal of Experimental Biology，2002，205（1）：55-70.

［7］ANDERSEN A，BOHR T，SCHNIPPER T，et al.Wake structure and thrust generation of a flapping foil in two-dimensional flow［J］.Journal of Fluid Mechan?ics，2017，812：R4-1-R4-12.

［8］胡合文.波浪滑翔機的水動力分析［D］.哈爾濱：哈爾濱工程大學，2015.HU H W.Hydrodynamic analysis of wave glider［D］.Harbin：Harbin Engineering University，2015（in Chi?nese）.

［9］張曉慶，王志東，張振山.二維擺動水翼仿生推進水動力性能研究［J］.水動力學研究與進展（A輯），2006，21（5）：632-639.ZHANG X Q，WANG Z D，ZHANG Z S.Hydrodynam?ic study of bionic propulsion for 2D flapping foil［J］.Journal of Hydrodynamics（Ser.A），2006，21（5）：632-639（in Chinese）.

［10］劉煥興，蘇玉民，龐永杰.非正弦擺動對水翼水動力性能的影響［J］.華中科技大學學報（自然科學版），2016（2）：15-20.LIU H X，SU Y M，PANG Y J.Effects of the nonsinu?soidal oscillating on the hydrodynamic performance of the hydrofoil［J］.Journal of Huazhong University of Science and Technology（Natural Science Edition），2016（2）：15-20（in Chinese）.

［11］OHTAKE T，NAKAE Y，MOTOHASHI T.Nonlinear?ity of the aerodynamic characteristics of NACA0012 aerofoil at low Reynolds numbers［J］.Journal of the Japan Society for Aeronautical and Space Sciences，2007，55（644）：439-445.

［12］YONEMOTO K，TAKATO K，OCHI H，et al.Kutta condition violation in two-dimensional NACA0012 airfoil at low Reynolds number［C］//26th AIAA Ap?plied AerodynamicsConference.Hawaii： AIAA，2008：6399.

［13］READ D A，HOVER F S，TRIANTAFYLLOU M S.Forces on oscillating foils for propulsion and maneu?vering［J］.Journal of Fluids and Structures，2003，17（1）：163-183.

［14］TRIANTAFYLLOU G S，TRIANTAFYLLOU M S，GROSENBAUGH M A.Optimal thrust development in oscillating foils with application to fish propulsion［J］.Journal of Fluids and Structures，1993，7（2）：205-224.

［15］SCHOUVEILER L，HOVER F S，TRIANTAFYL?LOU M S.Performance of flapping foil propulsion［J］.Journal of Fluids and Structures，2005，20（7）：949-959.

［16］Liquid Robotics.Converting wave motion into propul?sion［EB/OL］.（2014-09）［2017-06］.https：//www.liq?uid-robotics.com/platform/how-it-works/.