魏子天，洪 亮，劉新月，南 栩

(1.中國船舶科學研究中心國家船舶噪聲與振動重點實驗室，江蘇 無錫 214000；2.南京理工大學 能源與動力工程學院， 南京 210000)

1 引言

水翼是船舶上調整航向的結構，水翼處在密度較大的水流中，水流會對水翼產(chǎn)生非定常激勵從而導致水翼振動。近年來造船技術飛速發(fā)展，船舶航行速度越來越高。在高流速下，一種特殊的振動“顫振”愈發(fā)引起研究者的注意，顫振是結構阻尼無法消耗流場的輸入能而導致振幅不斷擴大的振動。顫振最早發(fā)現(xiàn)于航空器，早期航空器多采用展弦比較大的機翼，顫振現(xiàn)象頻發(fā)，對于機翼顫振的研究也較多，現(xiàn)階段水翼顫振的相關研究也多由空氣動力學領域的相關理論改進而來。

目前常用的顫振預測方法有ONERA失速氣動力模型，展志煥等[1]基于該模型，建立二元翼型的非線性顫振運動的方程，獲得非線性顫振時域響應曲線，由二分法尋找翼型的臨界顫振狀態(tài)。在發(fā)生顫振時機翼氣動力具有非線性，趙永輝等[2-3]使用該模型對三維彎扭耦合梁進行了失速攻角下的顫振狀態(tài)分析。此外，Theodorsen[4]提出的基于非定常氣動力線性化精確解的Theodorsen非定常氣動力模型應用也較多，劉胡濤[5-6]由該理論的理論解求解水翼算例的水彈性不穩(wěn)定邊界條件，并使用龍格-庫塔法進行了驗證。李景奎等[7]使用V-g法對氣動彈性運動方程的特征值進行求解，獲得機翼的臨界顫振速度。在試驗方面，Song[8]采用實驗的方法，對三維機翼上的非定常水動力及其水動彈性特性進行了實驗測定。并通過水洞試驗對不同機翼的顫振速度進行預測。Jewell[9]在美國的戴維-泰勒拖曳水池，通過設置2個自由度彈簧系統(tǒng)，以結構阻尼比為變量做了大量關于水翼顫振的試驗。

對于水翼顫振的研究理論仍較為匱乏，本文引入橋梁振動領域應用較多的“Scanlan顫振導數(shù)理論”，使用CFD模擬的方法，根據(jù)理論求取水翼的顫振導數(shù)，進而求得水翼的臨界顫振狀態(tài)，最后采用非線性振動求解方法紐馬克-β法對求得的結果進行驗證與分析。

2 水翼計算模型

采用二維模型進行計算，選取單位長度的水翼節(jié)段，設定水翼有沿著翼展方向“豎彎”方向運動自由度，以及沿旋轉軸轉動的“扭轉”自由度，將這2個自由度方向上的運動進行結合，來表達二維水翼的振動情況。2個自由度上的運動由彈簧-阻尼系統(tǒng)進行控制，如圖1所示。

使用NACA0015翼型，水翼弦長b=0.35 m，攻角為0°。水翼的豎彎及扭轉2個方向上水翼的固有頻率分別為ωh=4.37 Hz及ωα=2.95 Hz。水翼每延米的質量為m=206 kg，扭轉慣性矩為Im=12.11 kg·m2/m。剛心設置在距離水翼前緣0.42倍弦長處，與水翼的重心位置重合。

圖1 2個自由度模型

應用Fluent軟件進行計算。水翼的振動涉及到邊界的運動，需要用到“動網(wǎng)格”，本文采用Fluent前處理軟件ICEM進行建模和網(wǎng)格劃分，網(wǎng)格使用“剛性邊界層運動區(qū)域+動網(wǎng)格變形區(qū)域+外流場”的劃分策略。圖2給出了網(wǎng)格劃分的大致示意圖。流場網(wǎng)格如圖3所示，水翼前、后緣網(wǎng)格如圖4所示。

圖2 網(wǎng)格區(qū)域示意圖

圖2中，A區(qū)域為剛性運動區(qū)域，該區(qū)域與水翼一起運動。該區(qū)域的半徑為b，采用四邊形結構化網(wǎng)格劃分，以保證水翼運動時計算的收斂性。B區(qū)域為動網(wǎng)格區(qū)域，采用了三角形非結構化網(wǎng)格進行劃分。該區(qū)域的網(wǎng)格更新策略為“Smoothing”和“Remeshing”。C區(qū)域為外流場區(qū)域，流場的尺寸為20b×26b。流場設置的較大，可以避免流場邊界與水翼附近流場相互干涉而影響計算的準確性。

圖3 流場網(wǎng)格

圖4 水翼前、后緣網(wǎng)格

3 水翼顫振導數(shù)求解

Scanlan顫振導數(shù)理論內容較多，具體可參考文獻[10]，本文僅對關鍵點進行介紹。

對于均勻水平來流中攻角為0°的理想薄平板作頻率為ω微小簡諧運動時，薄平板的微幅振動會對平板上下表面的氣流產(chǎn)生擾動。氣流的運動反過來也會對平板產(chǎn)生作用力。由于平板自身運動所產(chǎn)生的的力稱為自激力，自激力的理論解可表達為：

(1)

(2)

式(1)—(2)中：L、M分別為平板所受到的升力及扭矩；ρ為標準溫度下流體密度；b為薄平板半寬，2b=板寬B；U為流體流速；h；α為平板的豎向位移和扭轉角；K為無量綱折算頻率，K=ωB/U，ω為振動圓頻率。

式(1)—(2)是基于薄平板推導得出，Scanlan認為無論是鈍體還是流線體都滿足該式，并對式(1)—(2)進行了改編，引入一組僅與截面形狀有關的無量綱參數(shù)“顫振導數(shù)”，以實現(xiàn)模型力和原型力的轉換，此時自激力公式可改寫為：

(3)

(4)

根據(jù)所求得的水動力數(shù)據(jù)，采用最小二乘法即可求解顫振導數(shù)，最終求得水翼的顫振導數(shù)值如表1所示。

表1 顫振導數(shù)值Table 1 Flutter derivatives

4 水翼臨界顫振狀態(tài)求解

求解顫振導數(shù)的目的是預測水翼的臨界顫振狀態(tài)，本文使用Scanlan于1951年提出的計算二維截面臨界顫振速度的方法[10]。

(5)

將Scanlan自激力式(3)和式(4)代入2個自由度截面運動方程(5)，并引入無量綱時間概念：s=tU/B，對截面運動方程進行無量綱化，可得：

(6)

(7)

定義未知函數(shù)X=ω/ωh，代入上述方程進行整理可得一個關于X的4次多項式。假定在顫振狀態(tài)下X總為實數(shù)，可以得到實部和虛部2個方程為：

A4RX4+A3RX3+A2RX2+A1RX+A0R=0

A3lX3+A2lX2+A1lX+A0l=0

(8)

式(8)中各項系數(shù)值為：

(9)

式(8)為一元高次方程，可對其進行分別求解，實部方程和虛部方程分別可求得4個解和3個解。舍去負解和不合理的解，繪出實部解XR和虛部解Xi隨折算速度Vr變化的曲線，它們的交點(XC,VrC)即代表臨界顫振狀態(tài)，臨界顫振速度計算式為：

UC=BωhXCVrC

(10)

對于水翼算例，由求得的顫振導數(shù)，結合式(9)可以求得式(8)的各項系數(shù)，如在折算速度為3.33時，實、虛部方程為：

1.118 6X4-0.002 7X3-1.637 8X2+0X+0.456 4=0

-0.393 0X3-0.018 4X2+0.205 3X+0.011 3=0

利用高斯消元法對上式進行求解：實部方程解為：-1.042 1；-0.612 5；1.045 7；0.611 3。虛部方程解為：0.728 3；-0.720 3 -0.055 1，其他折算速度下求解同理。實、虛部方程的解舍去負解并進行二次多項式擬合，實、虛部解擬合曲線如圖5所示。

圖5 實、虛部解擬合曲線

擬合曲線的交點為(12.72，0.38)，由式(10)求得臨界顫振速度為：

Uc=BωhXCXrC=0.35×4.37×12.72×0.38=7.39 m/s

即求得該水翼的臨界顫振速度為7.39 m/s。

5 驗證與分析

使用紐馬克-β法對第4節(jié)求得的結果進行驗證。紐馬克-β法是求解非線性振動的一種方法，本文根據(jù)其原理，編寫UDF嵌入Fluent軟件中實現(xiàn)水翼的流固耦合計算。判斷水翼是否發(fā)生顫振的依據(jù)如圖6所示。

圖6 臨界顫振狀態(tài)判斷依據(jù)

給予水翼一個很小的位移，使水翼處于不平衡狀態(tài)，待流場穩(wěn)定后釋放，水翼若在自激力的作用下，發(fā)生振幅越來越大的振動，則水翼進入顫振狀態(tài)。編寫Profile文件使水翼在前一秒內兩自由度皆做頻率為0.75 Hz的小幅正弦運動，待1 s時，豎彎和扭轉2個自由度的振幅皆達到最值時，將處于不平衡狀態(tài)的水翼釋放，啟動紐馬克-β算法進行計算。

圖7—圖11給出了水翼在各流速下的豎向振動及扭轉角時程圖。

圖7 扭轉、豎彎振動時程圖(U=3.0 m/s)

圖8 扭轉、豎彎振動時程圖(U=7.0 m/s)

圖9 扭轉、豎彎振動時程圖(U=7.3 m/s)

圖11 扭轉、豎彎振動時程圖(U=8.0 m/s)

由圖7—圖11可知，前1 s內水翼豎彎及扭轉方向皆由Profile驅動做規(guī)律的正弦運動。在流速為3 m/s時，水翼的振幅迅速衰減，而后做微幅振動以抵消流場輸入的能量。當流速提升至7 m/s時，水翼被釋放瞬間豎彎及扭轉方向振幅明顯要比U=3 m/s要大，振幅衰減速度比U=3 m/s時要慢。流速為7.3 m/s時，水翼豎彎及扭轉方向皆做等幅簡諧振動。當流速提升至7.5 m/s，2自由度方向的振動總體呈現(xiàn)緩慢發(fā)散態(tài)勢。U=8 m/s時，發(fā)散速度要比U=7.5 m/s時快的多，結構已經(jīng)發(fā)生顫振，在計算中由于振幅發(fā)散過快，水翼運動至動網(wǎng)格區(qū)域邊緣導致計算報錯。根據(jù)以上運動圖像及分析可以推斷，水翼的臨界顫振速度約為7.3 m/s。

圖12—圖14給出了水翼在流速為7、7.3以及7.5 m/s時的升力、力矩時程圖。

由圖12—圖14可知，水翼所受升力在1 s后的幅值差距不大，約為3 800 N左右，但之后流速為7 m/s時的水翼受到的升力越來越小，7.5 m/s時的水翼受到了越來越大的升力。流速為7.3 m/s時水翼受到了呈穩(wěn)定大小以正弦規(guī)律變化的升力。水翼所受力矩規(guī)律與之相似。可見當結構進入顫振狀態(tài)時，受到的升力及力矩會隨時間推移而越來越大，相應的振動振幅也會越來越大，最終導致結構失穩(wěn)。

圖15—圖17給出了水翼在流速為7、7.3以及7.5 m/s時的扭轉、豎彎振動相軌線圖。

圖12 升力及力矩時程圖(U=7.0 m/s)

圖13 升力及力矩時程圖(U=7.3 m/s)

圖14 升力及力矩時程圖(U=7.5 m/s)

圖15 扭轉、豎彎振動相軌線圖(U=7.0 m/s)

圖16 扭轉、豎彎振動相軌線圖(U=7.3 m/s)

圖17 扭轉、豎彎振動相軌線圖(U=7.5 m/s)

由圖15—圖17可知，在U=7.0 m/s時，水翼沒有達到臨界顫振速度，水翼的動能及勢能在運動過程中逐漸被消散，水翼在豎彎及扭轉自由度上的振幅及振動速度都在不斷減小，相軌線不斷向系統(tǒng)平衡點接近，直至達到穩(wěn)定狀態(tài)。當U=7.3 m/s時，豎彎及扭轉振動演化為典型的極限環(huán)振動，水翼的能量在動能和勢能之間不斷轉化。當流速提升至7.5 m/s時，水翼已達到顫振速度，系統(tǒng)的振幅及振動速度在流場不斷做正功的原因下不斷增大，直至結構破壞。

對水翼在2個自由度上的振動頻率進行了統(tǒng)計，如圖18所示。

圖18 扭轉、豎彎振動頻率

由圖18可知，2自由度的振動頻率都隨著流速的增加而逐漸減小。流速較低時，2自由度的振動頻率差距較大，隨著流速增加，振動頻率趨向一致；當流速大于或等于臨界顫振速度時，2自由度的振動頻率相等。印證了顫振是由具有2個自由度以上的結構物以同一頻率耦合振動的現(xiàn)象。

經(jīng)過以上分析，由紐馬克-β法求得的水翼臨界顫振速度為7.3 m/s，與Scanlan顫振導數(shù)理論求得的結果基本一致。

6 結論

1)引入Scanlan顫振導數(shù)理論，由該理論成功求解水翼的臨界顫振速度并進行了驗證，證明Scanlan顫振導數(shù)理論可以應用在水翼的臨界顫振狀態(tài)求解。

2) 由紐馬克-β法對水翼的振動狀態(tài)進行了分析，水翼一旦進入顫振狀態(tài)，振幅會不斷擴大直至破壞。升力和力矩的變化情況和振幅一致。

3) 水翼進入臨界顫振狀態(tài)時，動能和勢能總和不變，并在不斷轉化中。2自由度的頻率一致，印證了顫振是由具有2個自由度以上的結構物以同一頻率耦合振動的現(xiàn)象。