低質量比舵系統流致振動特征研究

王人鳳， 尤云祥， 陳 科， 胡曉峰

(上海交通大學 海洋工程國家重點實驗室 高新船舶與深海開發裝備協同創新中心，上海 200240)

顫振是一種由彈性力、慣性力、阻尼力和流體自激力共同作用引起的彈性不穩定現象。由于這類振動有可能演化為發散振動，故該現象一旦發生，結構就存在損壞的風險。經典顫振理論認為質量比(μ=m/πρb2l，m為系統質量，ρ為流體密度，b為舵的半弦長，l為舵的展長)低于臨界質量比時顫振便不會發生，但一些試驗表明，當水下升力體的質量比較小時也有可能發生顫振[1]。顫振理論最初應用于航空領域，后又被推廣到船舶領域。機翼的臨界顫振速度隨著機翼質量、質心到剛心的距離以及結構阻尼系數的增加而增大。當結構具有立方非線性剛度時，系統會在一定的速度范圍內做穩定的極限環運動，并且極限環的幅值會隨著流速的增加而增大。但是，當流速增大到一定程度時，系統進行發散振動。在求解振動方程時，不同的初始積分條件也會對系統極限環振動的幅值產生影響[2]。

張琪昌等[3]針對含立方非線性剛度二元機翼系統顫振時的極限環振動開展了相關研究，對分岔點類別和平衡點的性質進行了定性的分析，進而研究了系統參數尤其是機翼外形變化對極限環顫振穩定性和幅值的影響。發現當平衡點為穩定的焦點時，機翼會進行衰減振動，滿足設計要求。但是，在平衡點的附近也有可能出現穩定的極限環顫振運動，甚至是不穩定的極限環顫振運動。在工程上，對于小振幅穩定極限環顫振，只要不引起疲勞破壞就可以滿足要求；而對于不穩定的極限環顫振運動，必須避免。設計中可以通過改變系統的參數來抑制振幅和臨界顫振速度的大小，從而避免不穩定極限環顫振運動的產生。趙海等[4-5]求解某彈性系統的顫振速度并進行了數值驗證，分析了系統各參數對顫振速度的影響，并討論了非線性剛度系數對系統沉浮和俯仰自由度顫振極限環幅值的影響。結果表明，忽略非線性阻尼項會對超聲速流二元機翼的顫振計算造成誤差。其研究還顯示，雖然系統運動微分方程具有對稱性，但兩個自由度上的結構參數對系統臨界顫振速度的影響卻不具有對稱性。張瑜等[6-7]利用變步長四階龍格庫塔法對二元機翼系統的運動微分方程進行數值模擬，對不同流速區間內平衡點、極限環個數以及穩定性進行研究。同時還對極限環分叉問題進行了數值模擬，結果表明系統的極限環也有可能產生和平衡點一樣的分岔現象，且隨著無量綱流速的增大，極限環振動會經過倍周期分叉轉化為混沌運動。

Hamdani等[8-9]將有限元方法應用于翼型的動力特性分析中，一方面可以得到翼型表面的壓力分布，另一方面也可以得到翼型附近的流體分布。翼型在低雷諾數下同樣會出現運動速度突然增大或減小的情況，并且在此期間會產生一個較大的升力，導致翼形的進行明顯的沉浮和俯仰運動。在這種不穩定的運動發生期間，會在翼型的上下表面不斷出現一層層范圍更大的漩渦。Münch等[10]利用CFX模擬了NACA0009翼型在湍流中的動力特性。翼型的初始攻角為2°，流速為5～15 m/s，計算結果與實驗數據吻合，誤差小于2%。Gnesin等[11]利用一種考慮旋轉葉片振動因素的三維流固耦合求解方法研究了渦輪轉子的振動問題。發現在葉片振動的幅值-頻率譜上存在高頻和低頻，但是這些頻率并不是轉子的倍頻。Amiralaei等[12]為分析NACA0012翼型在低雷諾數流體中的振動特性，分別選取攻角2°～10°，使用有限體積法求解Navier-Stokes方程，并將計算結果中的瞬時升力系數與利用Theodorsen理論求解的結果進行比較，發現雷諾數不僅對翼型的升阻力系數有影響，對翼型的振動特征也存在明顯的影響。Akcabay等[13]利用CFX求解非定常雷諾平均Navier-Stokes方程分析空泡對NACA66水翼水彈性穩定性的影響。隨著空泡量的增加，翼型的平均升力隨之減小而平均阻力增大，同時其平均變形幅度也隨之減小。此外，由流動引起的振動幅值和阻力隨著翼型剛度的減小而增大。劉胡濤等[14]利用流固耦合方法對二元水翼的運動進行數值模擬，分析了初始攻角、剛心位置以及來流速度對水翼振動特性的影響。使用大渦模擬方法計算高雷諾數下水翼繞流場以及流場力作用于兩自由度剛體上所導致的周期性沉浮和俯仰運動。采用龍格庫塔法求解剛體水翼的運動方程，位移參數作為下一時間步流場計算的邊界條件。結果表明，水翼振動狀態對系統的初始值具有依賴性，在沒有擾動的情況下，水翼在零攻角時始終保持微幅振動。隨著初始攻角的增大，振動平衡位置越偏離初始位置，且沉浮運動幅值衰減較明顯。來流速度對水翼振動影響較為顯著，并且當來流速度增大到一定程度時，水翼會出現顫振現象。此外，還發現水翼顫振產生的條件比較苛刻。盡管如此，隨著流速的增大，系統產生的振動加劇現象仍不容忽視。

若舵具有較小展弦比，利用兩自由度系統動力方程分析振動特性時需要考慮其三維效應，所以有必要對理論模型進行相應的參數修正。此外，鑒于低質量比(μ<1)舵系統顫振相關研究較少，通過數值模擬方法分析顫振發生時舵周圍的流場變化也很有必要。

1 模型和計算方法

1.1 小展弦比兩自由度系統動力模型

舵系統由一個截面為NACA0017的舵、一根剛性軸以及與其連接的彈性裝置組成，剛性軸與舵固接。其中，舵的展長l=0.39 m，兩端弦長不同，舵根側面(較大的側面)的弦長為0.3 m，舵稍側面(較小的側面)的弦長為0.2 m，則中截面的弦長c=0.25 m。可以將舵系統簡化為經典兩自由度舵模型，如圖1所示。其主要結構參數如表1所示。

圖1 有限展長舵的兩自由度振動理論計算模型Fig.1 Two-degree-freedom vibration mechanical model for a foil with a finite span

m/kgμAOA/(°)kh/(N·m-1)kα/(N·m·rad-1)xα/(N·m-1)a7.60.451.4×1062820.04-0.38

當舵做簡諧運動時，會受到升力L(向上為正)和俯仰力矩Tα(迎流抬頭為正)的共同作用，并在流體動力的激勵下發生兩個自由度的運動，一個是舵-軸部件的沉浮運動，位移為h，向下為正；另一個是舵隨剛性軸轉動而產生的俯仰運動，俯仰角度為α，迎流抬頭為正。

(1)

式中：V∞為流速；C(k)為Theodorsen函數； 式中與C(k)無關的項來自慣性效應。

考慮展弦比對結構振動的影響，對該式進行修正。代表慣性效應的項添加附加質量修正系數ε，而環量項添加環量修正系數δ，得到

(2)

其中，

(3)

式中：AR=l/(ccosAOA)，AOA為初始攻角；τ為舵的形狀參數。

通過Fourier變換和引入Wagner函數[15]，結合舵的兩自由度任意運動微分方程

(4)

可以得到

(5)

其中，

(6)

將式(5)寫為矩陣形式，有

(7)

(8)

從而得到兩自由度系統任意運動的無量綱動力方程為

(9)

1.2 流體模型

使用CFD(Computational fluid dynamics)商業軟件中的彈簧構件生成系統特定的支撐剛度，并利用UDF(User Defined Function)通過添加與俯仰運動方向相反的俯仰力矩來模擬相應的扭轉剛度。對于流體域，左側以及前后壁面為速度入口，上下為固壁面，右側為壓力出口，如圖2所示。

圖2 數值模擬幾何圖Fig.2 Geometric model of numerical simulation

當低速航行時，舵系統周圍的流體可以視作不可壓流體。質量守恒方程和動量守恒方程組成了黏性不可壓流體的基本控制方程，而控制方程正是這些守恒定律的數學描述。現實中的湍流流動是一個隨機變化的脈動量，由于很難在計算中完全反映這些脈動量，所以引入雷諾時均的概念，將任一流動變量的瞬時值進行分解，并看作平均值和脈動值之和。

三維瞬態不可壓流體質量守恒方程為

(10)

時均值的動量守恒方程，即雷諾時均納維爾—斯托克斯(Reynolds Average Navier-Stokes, RANS)為[16]

(11)

式中： (u,v,w)為流速在直角坐標系中的三個分量;ρ為流體的密度;p為壓力;μ為流體的動力黏性系數；g為重力加速度。

由Boussinesq假設，引入湍流動力黏性系數μt，可以得到雷諾應力張量和流體時均速度梯度之間的關系

只要選擇合適的湍流模型將μt確定，就可以使RANS封閉。本文選擇的SST(Shear Stress Transport)k-ω模型(經剪切應力修正的兩方程模型)結合了標準兩方程模型k-ω和k-ε的優點，在近壁面區域利用k-ω模型，而在遠離附面層的區域使用k-ε模型，可以精確地計算物體表面流體的逆壓梯度以得到物體附近流體的細微變化[17-18]。

1.3 流固耦合設置

把舵系統的運動看作是一種剛體運動，其重心位于剛體運動局部坐標系的原點。為了保持和實際情況一致，在剛心的垂向位置添加了彈性系數為kh的彈簧組構件，并通過UDF設置了舵系統繞剛性軸方向的俯仰力矩，其表達式為

Tα=-kαθ

(12)

式中：θ為舵轉動的弧度。

舵系統在運動過程中所受到的回復力F和俯仰力矩Tα為

(13)

(14)

式中： [τ]和p為作用的剪應力和壓力；r-rc為舵表面任意網格中心到舵重心的距離(用位置矢量表示)；n為舵表面的外法線方向。

得到舵系統的受力之后，需要進一步計算其運動姿態。引入舵系統的兩自由度運動方程可得

(15)

(16)

式中：v為舵系統沉浮運動的速度；M為慣性矩張量在繞軸方向上的分量；ω為舵系統俯仰運動角速度。由式(15)和式(16)可以求得舵系統運動的速度、角速度以及位移，進而確定其運動姿態。在求解過程中，每個時間步長內對離散方程進行多次迭代，當滿足收斂條件時輸出速度和角速度等參數，并即時更新網格節點的位置，從而實現舵系統在流體中運動的數值模擬。

1.4 網格劃分

在舵和流體耦合的過程中，網格移動是耦合界面數據傳遞的關鍵。當結構的運動形態存在較大變化時，網格移動速度和準確度會對數值模擬結果的效率和精度產生巨大影響。重疊網格技術(Overlapping Grid Method，OGM)的使結構網格的自動更新成為可能。OGM將覆蓋整個計算域的母網格和運動物體的子網格進行疊加來描述物體間的相對運動，各區域中的計算網格獨立生成，流場信息通過插值函數在重疊區邊界進行匹配和耦合，從而實現網格的自動更新。

依據舵模型的尺寸建立相應的幾何模型，并劃分面網格，網格數為23 834個，如圖3所示。

圖3 舵的概念、實物、幾何建模和網格圖Fig.3 Conceptual graph, photograph, geometric model and grid of the hydrofoil

圖4為體網格的劃分，網格數為122萬個。其中，對舵表面附近的網格進行了加密。在計算域內采用動網格技術，舵附近和外邊界分別設置了不同的網格密度，舵附近的網格最密，而越接近外邊界網格越稀疏。重疊網格實現了不同網格區域的數據傳遞，而通過設置疏密程度不同的網格，既可以精確反映舵表面附近的流場變化，又可以提高計算速度。

圖4 網格劃分Fig.4 Mesh generation

2 結果與分析

2.1 振動特性分析

圖5為舵系統振動理論值的V-g曲線，在給定的系統參數下存在一個臨界航速VF，當V∞VF時，沉浮分支所對應的g為負值，這時系統進行發散振動。

圖5 舵系統振動計算值V-g曲線圖Fig.5 Computed V-g curve for the vibration of the hydrofoil system

舵系統沉浮和俯仰振動幅度的實驗值變化曲線，如圖6所示。系統的兩種振動成份均隨著V∞增加而增大。沉浮振動幅度曲線在V∞=2.37 m/s時出現一個明顯的波動，并且此流速所對應的頻譜圖上在原有頻率成份的波峰1附近出現了明顯的伴隨峰值。可以認定此時系統發生顫振，但是這種顫振與高質量比(μ?1)系統的顫振不同，其振動V-g曲線向下凸。從實驗現象來看，一方面由于阻尼的存在，另一方面V∞較低并且系統的質量較小，所以并沒有出現高速水翼以及機翼經歷顫振時的結構損毀。

圖6 舵系統振動幅度變化和振動頻譜圖Fig.6 Amplitude growth of the vibrations and frequency spectrums of the hydrofoil system

圖7為V∞=1.0 m/s時舵系統振動時歷曲線(見圖7(a)和圖7(c))和相軌線(見圖7(b)和圖7(d))。可見，系統的運動中存在著沉浮和俯仰兩個成份，分別是系統對支撐彈簧回復力和扭轉彈簧回復力矩的響應。由于V∞=1.0 m/s沒有達到系統的臨界顫振速度VF，所以系統進行衰減振動。系統的動能會在運動的過程中發生耗散，表現為系統沉浮和俯仰運動的幅度逐漸變小，直至達到穩定狀態。從系統運動相軌線可以看出，沉浮成份表現為圍繞原點做周期運動，由于振幅單調減小，所以相軌線不斷持續向原點靠近，最終在原點達到穩定狀態。而俯仰成份的周期運動圍繞著一個飄移的中心進行。從計算結果來看，舵系統流致運動的沉浮成份是一種簡諧振動，而俯仰成份為兩個頻率的耦合振動。

圖7 V∞=1.0 m/s時，系統振動時歷曲線和相軌線Fig.7 Time history cures and phase orbits of the hydrofoil system at V∞=1.0 m/s

經計算，舵系統的臨界顫振速度為VF=2.37 m/s，圖8(a)和圖8(c)為系統發生顫振時的時歷曲線圖。可見，沉浮成份做等幅簡諧振動，即運動的過程中不會發生能量耗散；而俯仰成份在經歷了短暫的飄移振動之后，也進行和沉浮成份相似的等幅振動。這是顫振理論中典型的振動形態，在忽略阻尼的作用時，系統將進行無消減的等幅振動，如果航行體在高速運動中發生這種振動有可能造成結構的扭曲、損壞，還有可能產生較大的噪聲導致降低航行體的性能。系統此時相應的運動相軌線，如圖8(b)和圖8(d)所示。沉浮振動和俯仰振動最終都演化為典型的極限環振動。因為系統在運動中不存在能量的耗散，所以能量在動能和勢能之間不斷轉化。

圖8 V∞=2.37 m/s時，系統振動時歷曲線和相軌線Fig.8 Time history cures and phase orbits of the hydrofoil system at V∞=2.37 m/s

圖9為利用未經參數修正的運動方程求解所得的非等截面系統發生顫振時的時歷曲線圖和相軌線。可見，系統運動的沉浮成份和俯仰成份都進行發散振動。沉浮成份以原點為中心做周期運動，但是與低流速時不同，運動軌跡不斷向外擴散；而俯仰成份的運動軌跡在出現短暫的飄移后，也轉化以原點為中心向外不斷做周期性的擴散。理論上，出現這種情況是由于系統在運動的過程中不但沒有發生能量耗散，還持續不斷地從流體中吸收能量。這些能量由系統的動能轉化為勢能，又由勢能轉化為動能，如此循環往復便出現了系統運動過程中沉浮和俯仰成份振幅不斷擴大的現象。與圖8進行比較，可知在未進行參數修正的前提下，利用兩自由度系統任意運動方程求解會導致VF偏小。

當來流速度V∞=5.0 m/s時，非等截面舵系統振動的時歷曲線和相軌線，如圖10所示。可見，系統振動的沉浮成份和俯仰成份都呈現發散狀態，并且幅值迅速擴大。圖10(b)和圖10(d)為系統運動的相軌線圖，沉浮成份和俯仰成份都以原點為中心做周期運動，運動軌跡不斷向外擴散，并且擴散的速度不斷增大。如果不考慮阻尼的影響，現實中的舵系統若在此流速下進行發散振動，則很可能因為沉浮和俯仰振幅的激增而遭到損壞。

圖9 V∞=2.37 m/s時，利用未修正方程求解的時歷曲線和相軌線Fig.9 Time history cures and phase orbits of the hydrofoil system via unmodified equation for the vibrations at V∞=2.37 m/s

圖10 V∞=5.0 m/s時，系統振動時歷曲線和相軌線Fig.10 Time history cures and phase orbits of the hydrofoil system at V∞=5.0 m/s

2.2 繞流場分析

圖11為V∞=1.0 m/s時舵附近的流場變化圖。在運動最初(見圖11(a))，舵的前緣出現零星的細小渦，但是并不明顯，舵的后緣附近出現了細小的尾渦。t=0.05 s(見圖11(b))時，在舵的前緣出現了明顯的渦量變化，同時舵的后緣出現明顯脫落渦。t=0.1 s(見圖11(c))時，在舵上表面從前緣到剛心(約1/4弦長處)之間出現較大的渦量，而后3/4的上表面有形成多個細小的渦。t=0.15 s時(見圖11(d))，在舵的上表面從剛心到后緣出現了明顯的渦層。t=0.5 s(見圖11(e))時，在舵的上表面，前緣和剛心之間的渦量繼續增大，剛心和后緣之間的渦層不斷擴大；在舵的下表面，也有出現渦層的趨勢。由于系統的運動發生衰減，所以t=4 s(見圖11(f))時舵趨于穩定，俯仰運動幅度銳減導致其尾部脫落渦消失，而沉浮運動幅度銳減導致渦集中在舵的后半部附近。

圖11 V∞=1.0 m/s時，舵附近流場以及渦量云圖Fig.11 Vorticity contours of the hydrofoil at V∞=1.0 m/s

圖12為V∞=VF=2.37 m/s時舵附近的流場變化圖。可見，在運動最初(見圖12(a))，和衰減運動時相似，舵的后緣附近出現細小的尾渦，同時前緣也出現了零星的渦。t=0.05 s(見圖12(b))時，舵的上表面前緣和半弦長之間出現了明顯的渦量變化，并且剛心位置的渦最為突出，同時在尾部出現一系列明顯的脫落渦。t=0.1 s(見圖12(c))時，在舵的上表面，前緣和剛心之間的渦不斷向后緣擴散，并在剛心和后緣之間形成了明顯的渦層。此時舵的上表面形成一個完整的渦層，這與衰減運動時不同；在舵的下表面，在剛心和后緣之間也有形成渦層的趨勢。t=0.15 s(見圖12(d))時，舵上表面的渦層繼續向后擴散，并與尾渦融合。t=2 s(見圖12(e))時，舵出現大角度的俯仰運動，此時前緣和剛心之間出現加大的渦量，而剛心和后緣之間的渦層已經完全消失。t=4 s(見圖12(f))時，只存在前緣附近的一個巨大渦，舵的上下表面幾乎都沒有渦形成。由于已經忽略了阻尼的作用，所以舵的俯仰運動繼續加劇，而非經典的等幅運動。這不僅會大大降低航行體的操控性能，甚至有可能造成舵系統的結構損壞。

圖12 V∞=2.37 m/s時，舵附近流場以及渦量云圖Fig.12 Vorticity contours of the hydrofoil at V∞=2.37 m/s

3 結 論

本文通過對經典兩自由度系統任意運動的無量綱動力方程進行參數修正，在計算中引入了展弦比對流體動力的影響。兩自由度舵系統的流致振動實際上是動能和勢能相互轉化的過程，隨著V∞的增加相軌線的形狀越來越趨于扁平化，這意味著系統的動能可以轉化為更大的勢能。這是因為低流速時，系統的動能會產生耗散，致使系統進行衰減振動；當流速達到VF時，系統的動能可以全部轉化為勢能，勢能亦可全部轉化為動能，從而出現極限環形式的相軌線；而在高流速下，系統的動能不但不會發生耗散，還會在振動過程中不斷從流體中吸收能量并將其轉化為勢能，使得系統振動相軌線不斷向外擴散。

低質量比舵系統發生顫振時，系統的振動形態從穩態轉化為非穩態。在對應的頻譜圖上，系統的沉浮頻率或俯仰頻率附近會出現一個新的頻率，表明在這個過渡期間系統的運動受到了一定程度的擾動，但是并沒有像高質量比系統那樣形成一個突出的頻率，這是低質量比系統和高質量比系統顫振現象的主要區別。

通過數值模擬，發現在舵系統運動最初，舵的后緣附近會出現一系列細小的尾渦。并且，在舵的前緣會出現了明顯的流體分離，這種分離不斷向后緣擴散，同時會在舵的尾部出現一系列脫落渦。當系統進行衰減運動時，則尾渦會隨著系統的靜止而消失，舵上表面的渦層主要出現在1/4弦長處和后緣之間。如果系統進行發散運動，則舵會發生劇烈運動，其前緣附近產生巨大的渦，渦層遍布舵的整個上表面。并且在1/4弦長處會出現明顯層流分離現象，這與衰減運動不同。