水翼渦激振動的數值模擬研究

劉胡濤，張懷新，，姚慧嵐

（1.上海交通大學 船舶海洋與建筑工程學院，上海 200240；2.高新船舶與深海開發裝備協同創新中心（船海協創中心），上海 200240）

水翼渦激振動的數值模擬研究

劉胡濤1，張懷新1，2，姚慧嵐1

（1.上海交通大學 船舶海洋與建筑工程學院，上海 200240；2.高新船舶與深海開發裝備協同創新中心（船海協創中心），上海 200240）

對二維水翼結構進行流固耦合運動分析，利用大渦模擬方法計算高雷諾數下水翼繞流場，流場力作用于二自由度剛體上導致周期性的垂蕩和轉動，龍格庫塔法求解剛體水翼的運動方程，位移參數作為下一時間步流場計算的邊界條件，具體通過編譯自定義函數控制剛體運動和流場網格變化。探討了初始攻角、剛心位置以及來流速度對水翼振動的影響，發現在來流速度增大至某一數值時水翼發生了顫振現象，并在時域與頻域上對顫振的發生機理進行深入探討。

水翼繞流；流固耦合；顫振

0　引　言

隨著水面艦艇和潛艇航速的提升，水彈性振動對船體強度以及疲勞壽命產生的影響已引起了造船界的廣泛關注。艦艇附件和海洋結構部件在粘性流體中運動時，引起流動分離，生成旋渦，造成渦激振動，可能導致結構破壞。利用水彈性的研究結果，分析粘性流體通過翼型振蕩物的振動研究，國內外的研究并不多?，F有文獻結果大多是假定圓柱、機翼靜止不動的情況下進行的。多數情況下水翼不是完全固定不動的，而是在彈性軸處受到非線性約束，水翼在流體的作用下會產生渦的脫落等復雜現象，由此流體力作用于結構物本身產生耦合效應渦激振動。如水翼艇中的水翼，潛艇中的舵翼等船舶構件作為動力學系統，設計不當時，在流場中高速運動時將會產生有害的渦激顫振，導致結構破壞。

流固耦合研究早期集中在氣動彈性力學，發展至今已經有較為成熟的理論體系，其大部分理論結果都可直接應用于水彈性力學中。早在 20 世紀 30 年代Theodorsen［1］提出基于非定常氣動力線性化精確解的不可壓縮、無粘流的 Theodorsen 非定常氣動力模型，該模型同樣可以處理二維水翼的非定常力。Huang等［2］對NACA0012 翼型進行實驗分析，并探討了翼根連接處與翼梢對渦激振動的影響。Jung［3］通過實驗對低雷諾數下機翼繞流形成的渦街，發現機翼振動情況下的渦脫頻率相對于機翼靜止不動時要小。王囡囡［4］采用片條理論的簡化非定常氣動力方法研究機翼氣動彈性系統的動態特性，提出了一種基于實測柔度的顫振主動控制方法。另一方面至 20 世紀 70 年代后期，船舶水彈性力學理論和分析方法取得了顯著的進展，并廣泛應用到船舶與海洋工程領域的流固耦合分析中。早期的工作主要集中在試驗研究，Davidson 實驗室對水翼做了大量的試驗研究，主要集中在以 2 個自由度理想翼型為出發點，觀察測量水翼在水洞試驗中的振動、空泡、壓力分布等。Henry［5］進行了一系列的模型試驗，對水翼各參數以及臨界顫振速度的影響作了理論上探討。余志興［6-7］運用 N-S 方程數值模擬了粘性流場中二維機翼的水彈性振動，在此基礎上，劉曉宙［8］利用有均流時二維運動物體的聲輻射方程，研究了大攻角下二維水翼的聲輻射，結果表明：渦脫頻率和機翼的固有振動頻率一致時，聲輻射達到最大。Chae［9-10］對2個自由度下的二維水翼進行了一系列系統的研究，探討了附加水質量、水動力阻尼、非線性流固耦合響應等對水翼穩定性的影響。

本文建立了 2 個自由度下二元水翼的流固耦合模型，分別利用大渦模擬方法和龍格庫塔法求解流體結構方程，界面數據傳遞通過自編函數實現。對NACA0012 翼型在不同工況下的繞流振動進行分析，并分別在時域和頻域上對振動情況進行分析。

1　數學模型

1.1剛體模型

本文將探討懸臂式、二維翼型在不可壓縮流體中的動態響應與穩定邊界，假定水翼僅在翼展方向存在2 個自由度運動，即對應船舶六自由度中的垂蕩與縱搖，如圖1 所示。

圖1　兩自由度水翼模型Fig.1　Two degrees of freedom hydrofoil models

其中，水翼弦長 c=2b，來流速度為 U，初始攻角為 α。EA 為水翼彈性中心（也即是二維翼型剛心），位于翼弦中點后 ab 處，a 為彈性中心到翼弦中點距離的無量綱量，b 為半弦長。翼型在彈性軸處由一個線彈簧和一個扭轉彈簧支撐在彈性軸處，垂直方向位移與旋轉均是在彈性軸上測量，h（t）向上為正，θ（t）逆時針為正，分別稱為水翼的垂蕩運動和縱搖運動。HC 為水動力中心，位于彈性中心前緣 eb 處，即是流體作用力合力的作用點。重心 CG 位于彈性中心后緣 xθb 處，xθ為重心到彈性中心距離的無量綱量。

兩自由度水翼無量綱運動方程為：

上述方程組可簡化為：

式中：Ms，Cs，Ks分別為水翼慣性、阻尼和剛度矩陣；F 為流場力矩陣；X 為兩自由度上廣義位移矩陣。

1.2流體模型

如圖1 所示，流場分為 3 個區域：水翼近場加密網格為動區域；遠場采用非結構網格為變形區域；延長水翼后方網格以觀察尾流。

圖2　水翼計算網格Fig.2　Grid of hydrofoil

基于商業軟件 Fluent 有限體積法對二維水翼進行瞬態繞流分析，利用 UDF 控制剛體與網格的運動。大渦模擬相較于平均雷諾數法可以提供豐富的大渦旋信息，且能計算出壁面處的壓強脈動。本文在對 Fluent中幾種經典湍流模型進行計算比較后選用 LES 作為后續的計算模型。

流場的連續性方程和不可壓縮 N-S 方程可以寫為：

對于 Smagorinsky-Lilly 模型

1.3流固耦合設置

建立 CFD 模型和計算結構動力學模型后，需要建立有效的界面耦合程序來聯系這 2 個單獨的計算模型。耦合界面不匹配網格間信息的準確交換，直接影響耦合計算精度，是獲得真實物理解的關鍵環節之一。

2　計算結果及分析

水翼結構參數取值如下：總弦長 c=1 m，垂蕩、縱搖固有頻率 ωh=ωθ=2s?1，中心對彈性軸的回轉半徑 r=0.8，質量比 μ=m/（πρb2）=14.0，無量綱系數a=?1,xθ=0.1。計算中流體與結構的時間步長均取為 0.01 s，是結構自振周期的 0.318%。

對 1 m/s 來流下渦激振動采用大渦模擬數值分析，剛心位置取在水翼前緣。圖中流場參數為 Re=106。

圖3　.耦合算法流程圖Fig.3　Coupled algorithm flowsheet

從圖4 可看出，10° 攻角下水翼后緣渦街脫落較為明顯，在水翼前緣和尾翼處均有渦的形成，以平行于來流的方向向尾部脫離，并在尾流中形成渦街。而0°攻角下由于水翼振幅很小，僅在尾翼處有渦的形成，且很快消散。目前關于翼型繞流的渦激振動研究，為了獲得翼型繞流周期性的渦發放，大多集中在大攻角工況下。本文以 10°攻角作為主要分析對象，發現大渦模擬能夠較好的捕捉復雜湍流中的細節，而在一般 k-ω等雷諾平均模型中，小攻角水翼繞流無法捕捉到尾部渦這一特征。

保持其他參數不變，改變初始攻角對比其對水翼繞流振動的影響，垂蕩運動和縱搖運動結果如圖6 所示。

顯而易見，攻角越大，其升力與力矩就越大，相應 2 個自由度上的位移就越大。從圖中可看到，不同攻角下的振蕩平衡位置隨攻角增大而增大。在 0° 攻角下，僅在水翼尾部形成較小的渦，且很快耗散掉，圖5 中可以看出，0° 攻角下，水翼振蕩位移的幅值非常小，且頻率較高，沒有體現出一定的周期性。本文嘗試將 0° 攻角來流速度提升至 4 m/s，2 個自由度上的位移仍保持非常小，且周期性較為紊亂。初步得出結論：0° 攻角在沒有初始擾動的情況下很難產生大幅振動的情況，攻角越大，振動平衡位置越偏離初始位置。且垂蕩運動幅值衰減較明顯，縱搖運動幅值衰減較慢。

圖4　NACA0012 水翼壓力、流場以及渦量云圖，Re=106Fig.4　Pressure，velocity and vorticity contours of hydrofoil Re=106

圖5　不同攻角下振動情況對比 U=1 m/sFig.5　Hydrofoil vibration under different attack angle

圖6　剛心位置對振動情況的影響 U=1 m/s，α=10°Fig.6　Hydrofoil vibration under different stiffness center

根據氣動彈性力學理論，重心位于剛心線之前就不會發生顫振，飛行器機翼設計時多采用在翼面前緣防止配重，使重心前移的方式來避免顫振。當重心位置一定時，應使剛心靠后一些（與重心靠近），以此來提高顫振速度。對于水翼而言，剛心位置同樣對顫振臨界條件有著至關重要的作用。通過升力力矩公式可以看出，繞彈性軸力矩的大小與剛心位置有關，剛心位置的改變直接影響到力矩的大小，進而影響到縱搖運動，由于流體-結構的耦合作用，結構運動狀態的改變將導致流場變化從而導致垂蕩運動的變化。圖6顯示 2 個自由度上位移都出現了一定幅度的衰減，且剛心位置對縱搖振動平衡位置的影響較垂蕩運動更加明顯。

從圖7可見，隨著來流速度增大，水翼在 2 個自由度上振動不斷加劇，可以看到 1.67 m/s 工況下縱搖自由度振幅隨著時間不斷增大，而垂蕩自由度的振幅經過小幅度減小后又繼續增大直到相對保持平穩。初步判定，本文所采用水翼的顫振臨界速度在 1.67 m/s附近。其中，1 m/s 時，升力和力矩在數值上保持微小的波動，可近似看作是固定力作用在水翼上，由于運動過程中的阻尼作用，振動幅度隨時間逐漸衰減。而 2 m/s時，升力和剛心處力矩幅度隨時間不斷增大，從而導致 2 個自由度上的位移也不斷增大，且此時水翼縱搖自由度瞬時水動力與彈性位移之間有約 180°的位相差，水翼已經進入顫振狀態，振動幅度不斷擴大從而導致結構的破壞。為了更直觀地觀察水翼顫振情況，做出以上 3 組速度下的相圖和對應，的譜分析圖（見圖8）。

來流速度較小時，流體對水翼做負功，水動力對水翼運動起阻尼作用，振動逐漸衰減最終穩定于某一平衡位置。速度為 1 m/s 時，由相圖分析看出系統未出現極限環運動，2 個自由度上的振幅不斷減小。由頻譜圖可知，垂蕩和縱搖自由度均由 2 個頻率成分疊加而成，其中垂蕩自由度主要頻率成分為 0.35 Hz、縱搖自由度主要頻率成分為 0.18 Hz 和 0.35 Hz。來流速度為 1.67 m/s 時，此時垂蕩自由度含有 2 個頻率分別為0.195 Hz 和 0.34 Hz?？v搖自由度的主要頻率為 0.195 Hz，由圖8可知，此速度下縱搖方向已不穩定，說明頻率較小的成分起到主要作用時系統會發生不穩定現象。速度為 2 m/s 時，2 個自由度上的主要頻率成分均為0.195 Hz，與 1.67 m/s 縱搖不穩定時的頻率相同，次要頻率成分為 0.3 Hz 且幅值遠小于 0.195 Hz，此時系統發生顫振，也同樣印證了上述說法，即較小頻率成分起主要作用時系統不穩定。表明來流速度增大至越過臨界速度后，水動力做正功，水翼發生自激振動，振幅不斷增大而導致振動發散，顫振是指結構物具有 2個自由度以上的以同一頻率耦合起來的振動現象［11］。由頻譜圖可以看出，隨速度增大，2 個自由度上的主要頻率成分趨于一致，且低頻成分逐漸起到主導作用，另一個頻率成分的影響逐漸降低。

圖7　不同來流速度下振動情況對比，α=10°，a=-1Fig.7　Hydrofoil vibration under different velocity

圖8　不同來流速度下縱搖運動相圖及頻譜Fig.8　Motion phase diagram and velocity spectrum in various velocity

圖9　不同流速下水翼的聲學指向性圖Fig.9　Acoustic directivity chart in various velocity

以水翼壁面上的脈動壓力作為偶極子噪聲源，對水翼渦激振動造成的水動力噪聲進行分析，求得圓心為翼弦中點，半徑為 1.5 m 圓上特征點的聲壓值，作出水翼流噪聲的聲學指向性圖。從圖9 可見，由 1 m/s增至 2 m/s 流噪聲增幅較大，而后增速減緩，這是由于顫振導致振幅的不斷增大導致的。且流噪聲在水翼首、尾位置處較小，這是因為來流和尾流對噪聲起著“掩蔽效應”，來流和尾流在靠近壁面處壓力變化梯度大，流場內的不均勻波動會抵消一部分流噪聲的傳播，因此在這一區域流噪聲數值相對較小。圖10 為特征點（2.0，0）處水翼低頻段聲壓級對比圖?？梢钥闯?，流噪聲主要集中在低頻區域，隨著速度增大，高頻成分逐漸占據主要作用，如 2 m/s 時，聲級主要集中在 0～5 Hz 之間，在 3 Hz 處出現峰值，速度增至 4 m/s時，聲級則集中在 0～20 Hz 之間，并出現了多個峰值。聯合頻譜分析可看出，振動頻率與噪聲頻率都隨來流速度增大而增大。

圖10　不同流速下特征點聲壓值對比圖Fig.10　SPL of feature points in various velocity

3　結　論

本文通過數值方法分析二元水翼在紊流中的振動情況，采用弱耦合的方式聯立大渦模擬湍流模型和剛體運動方程求解了水翼在 2 個自由度上的運動，在此基礎上分別對比了初始攻角、剛心位置以及來流速度對振動的影響，得出以下幾點結論：

1）大渦模擬方法能夠較好的捕捉高雷諾數、復雜湍流中的細節，模擬出小攻角水翼繞流后完整的渦街現象。

2）水翼振動狀態對系統的初始值具有依賴性，在沒有擾動的情況下，水翼在 0° 攻角時始終保持微幅振動，攻角越大，振動平衡位置越偏離初始位置，且垂蕩運動幅值衰減較明顯。

3）航行速度對水翼振動影響較為直觀，且隨著速度的不斷增大，水翼出現顫振現象，即振幅不衰減的自激振動。雖然水中顫振產生的條件比較苛刻，但是隨流速增大而產生的振動加劇現象仍不容忽視。

4）流噪聲隨來流速度增大而增大，其在低頻段的數值高于高頻段的數值，它能夠直觀地反應出結構振動對流場的擾動情況。

［1］THEODORSEN T，MUTCHLER W H.General theory of aerodynamic instability and the mechanism of flutter［R］.NACA Report No.496，NACA，1935.

［2］HUANG R F，LIN C L.Vortex shedding and shear-layer instability of wing at low-Reynolds numbers［J］.AIAA Journal，1995，33（8）: 1398-1403.

［3］JUNG Y W，PARK S O.Vortex-shedding characteristics in the wake of an oscillating airfoil at low Reynolds number［J］.Journal of Fluids and Structures，2005，20（3）: 451-464.

［4］王囡囡.二元機翼顫振及其主動控制的研究［D］.徐州: 中國礦業大學，2013.WANG Nan-nan.Research on flutter and active control for two-dimensional airfoil［D］.Xuzhou: China University of Mining and Technology，2013.

［5］HENRY C J.Hydrofoil flutter phenomenon and airfoil flutter theory［M］.New Jersey: Davidson Laboratory，1961.

［6］余志興.粘性流場中的水彈性計算［D］.上海: 上海交通大學，1999.YU Zhi-xing.Hydroelastic dynamics in viscous flow［D］.Shanghai: Shanghai Jiao Tong University，1999.

［7］余志興，劉應中，繆國平.二維機翼彈簧系統的渦激振動［J］.船舶力學，2002，6（5）: 25-32.YU Zhi-xing，LIU Ying-zhong，MIAO Guo-ping.Vortex-induced vibration of two-dimensional wing-spring coupled system［J］.Journal of Ship Mechanics，2002，6（5）: 25-32.

［8］劉曉宙，繆國平，余志新，等.流體通過渦激振動機翼的聲輻射研究［J］.聲學學報，2005，30（1）: 55-62.LIU Xiao-zhou，MIAO Guo-ping，YU Zhi-xin，et al.Research on the sound generation by flow around vortex-induced vibrated aerofoil［J］.Acta Acustica，2005，30（1）: 55-62.

［9］CHAE E J，AKCABAY D T，YOUNG Y L.Dynamic response and stability of a flapping foil in a dense and viscous fluid［J］.Physics of Fluids，2013，25（10）: 104106.

［10］CHAE E J.Dynamic response and stability of flexible hydrofoils in incompressible and viscous flow［D］.Michigan: University of Michigan，2015.

［11］程貫一，王寶壽，張效慈.水彈性力學——基本原理與工程應用［M］上海: 上海交通大學出版社，2013.CHENG Guan-yi，WANG Bao-shou，ZHANG Xiao-ci.Hydroelasticity——the basics with applications［M］.Shanghai: Shanghai Jiao Tong University，2013.

Research on vortex induced vibration of hydrofoil

LIU Hu-tao1，ZHANG Huai-xin1，2，YAO Hui-lan1
（1.School of Naval Architecture，Ocean and Civil Engineering，Shanghai Jiaotong University，Shanghai 200240，China；2.Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration（CISSE），Shanghai 200240，China）

The objective of this research is to investigate the fluid-structure interaction motion of a two-dimensional hydrofoil in viscous flow.Large-eddy simulation method is used to calculate the flow field around hydrofoil，and Runge-Kutta method is applied to solve equations of rigid body motion.Flow results and motion results exchange data at the hydrofoil wall surface，by compiling a user-defined functions to control rigid body motion and flow field grid update.The influence of the hydrofoil parameters，i.e.attack angle，elastic axis position and flow velocity，on vibration are discussed.Once the velocity increase to a certain value，hydrofoil flutter phenomenon occurred.Also，discussion on time and frequency domain is held to investigate the flutter phenomenon.

two-dimensional wing flow；fluid Structure Interaction；flutter

U661.1

A

1672-7619（2016）06-0007-07

10.3404/j.issn.1672-7619.2016.06.002

2015-12-09；

2016-03-17

國家自然科學基金資助項目 （11272213，51479116）

劉胡濤（1991-），男，碩士研究生，研究方向為船舶水動力學。