魏子天,洪 亮,劉新月,南 栩
(1.中國船舶科學(xué)研究中心國家船舶噪聲與振動重點實驗室,江蘇 無錫 214000;2.南京理工大學(xué) 能源與動力工程學(xué)院, 南京 210000)
水翼是船舶上調(diào)整航向的結(jié)構(gòu),水翼處在密度較大的水流中,水流會對水翼產(chǎn)生非定常激勵從而導(dǎo)致水翼振動。近年來造船技術(shù)飛速發(fā)展,船舶航行速度越來越高。在高流速下,一種特殊的振動“顫振”愈發(fā)引起研究者的注意,顫振是結(jié)構(gòu)阻尼無法消耗流場的輸入能而導(dǎo)致振幅不斷擴大的振動。顫振最早發(fā)現(xiàn)于航空器,早期航空器多采用展弦比較大的機翼,顫振現(xiàn)象頻發(fā),對于機翼顫振的研究也較多,現(xiàn)階段水翼顫振的相關(guān)研究也多由空氣動力學(xué)領(lǐng)域的相關(guān)理論改進(jìn)而來。
目前常用的顫振預(yù)測方法有ONERA失速氣動力模型,展志煥等[1]基于該模型,建立二元翼型的非線性顫振運動的方程,獲得非線性顫振時域響應(yīng)曲線,由二分法尋找翼型的臨界顫振狀態(tài)。在發(fā)生顫振時機翼氣動力具有非線性,趙永輝等[2-3]使用該模型對三維彎扭耦合梁進(jìn)行了失速攻角下的顫振狀態(tài)分析。此外,Theodorsen[4]提出的基于非定常氣動力線性化精確解的Theodorsen非定常氣動力模型應(yīng)用也較多,劉胡濤[5-6]由該理論的理論解求解水翼算例的水彈性不穩(wěn)定邊界條件,并使用龍格-庫塔法進(jìn)行了驗證。李景奎等[7]使用V-g法對氣動彈性運動方程的特征值進(jìn)行求解,獲得機翼的臨界顫振速度。在試驗方面,Song[8]采用實驗的方法,對三維機翼上的非定常水動力及其水動彈性特性進(jìn)行了實驗測定。并通過水洞試驗對不同機翼的顫振速度進(jìn)行預(yù)測。Jewell[9]在美國的戴維-泰勒拖曳水池,通過設(shè)置2個自由度彈簧系統(tǒng),以結(jié)構(gòu)阻尼比為變量做了大量關(guān)于水翼顫振的試驗。
對于水翼顫振的研究理論仍較為匱乏,本文引入橋梁振動領(lǐng)域應(yīng)用較多的“Scanlan顫振導(dǎo)數(shù)理論”,使用CFD模擬的方法,根據(jù)理論求取水翼的顫振導(dǎo)數(shù),進(jìn)而求得水翼的臨界顫振狀態(tài),最后采用非線性振動求解方法紐馬克-β法對求得的結(jié)果進(jìn)行驗證與分析。
采用二維模型進(jìn)行計算,選取單位長度的水翼節(jié)段,設(shè)定水翼有沿著翼展方向“豎彎”方向運動自由度,以及沿旋轉(zhuǎn)軸轉(zhuǎn)動的“扭轉(zhuǎn)”自由度,將這2個自由度方向上的運動進(jìn)行結(jié)合,來表達(dá)二維水翼的振動情況。2個自由度上的運動由彈簧-阻尼系統(tǒng)進(jìn)行控制,如圖1所示。
使用NACA0015翼型,水翼弦長b=0.35 m,攻角為0°。水翼的豎彎及扭轉(zhuǎn)2個方向上水翼的固有頻率分別為ωh=4.37 Hz及ωα=2.95 Hz。水翼每延米的質(zhì)量為m=206 kg,扭轉(zhuǎn)慣性矩為Im=12.11 kg·m2/m。剛心設(shè)置在距離水翼前緣0.42倍弦長處,與水翼的重心位置重合。
圖1 2個自由度模型
應(yīng)用Fluent軟件進(jìn)行計算。水翼的振動涉及到邊界的運動,需要用到“動網(wǎng)格”,本文采用Fluent前處理軟件ICEM進(jìn)行建模和網(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,采用四邊形結(jié)構(gòu)化網(wǎng)格劃分,以保證水翼運動時計算的收斂性。B區(qū)域為動網(wǎng)格區(qū)域,采用了三角形非結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分。該區(qū)域的網(wǎng)格更新策略為“Smoothing”和“Remeshing”。C區(qū)域為外流場區(qū)域,流場的尺寸為20b×26b。流場設(shè)置的較大,可以避免流場邊界與水翼附近流場相互干涉而影響計算的準(zhǔn)確性。
圖3 流場網(wǎng)格
圖4 水翼前、后緣網(wǎng)格
Scanlan顫振導(dǎo)數(shù)理論內(nèi)容較多,具體可參考文獻(xiàn)[10],本文僅對關(guān)鍵點進(jìn)行介紹。
對于均勻水平來流中攻角為0°的理想薄平板作頻率為ω微小簡諧運動時,薄平板的微幅振動會對平板上下表面的氣流產(chǎn)生擾動。氣流的運動反過來也會對平板產(chǎn)生作用力。由于平板自身運動所產(chǎn)生的的力稱為自激力,自激力的理論解可表達(dá)為:
(1)
(2)
式(1)—(2)中:L、M分別為平板所受到的升力及扭矩;ρ為標(biāo)準(zhǔn)溫度下流體密度;b為薄平板半寬,2b=板寬B;U為流體流速;h;α為平板的豎向位移和扭轉(zhuǎn)角;K為無量綱折算頻率,K=ωB/U,ω為振動圓頻率。
式(1)—(2)是基于薄平板推導(dǎo)得出,Scanlan認(rèn)為無論是鈍體還是流線體都滿足該式,并對式(1)—(2)進(jìn)行了改編,引入一組僅與截面形狀有關(guān)的無量綱參數(shù)“顫振導(dǎo)數(shù)”,以實現(xiàn)模型力和原型力的轉(zhuǎn)換,此時自激力公式可改寫為:
(3)
(4)
根據(jù)所求得的水動力數(shù)據(jù),采用最小二乘法即可求解顫振導(dǎo)數(shù),最終求得水翼的顫振導(dǎo)數(shù)值如表1所示。
表1 顫振導(dǎo)數(shù)值Table 1 Flutter derivatives
求解顫振導(dǎo)數(shù)的目的是預(yù)測水翼的臨界顫振狀態(tài),本文使用Scanlan于1951年提出的計算二維截面臨界顫振速度的方法[10]。
(5)
將Scanlan自激力式(3)和式(4)代入2個自由度截面運動方程(5),并引入無量綱時間概念:s=tU/B,對截面運動方程進(jìn)行無量綱化,可得:
(6)
(7)
定義未知函數(shù)X=ω/ωh,代入上述方程進(jìn)行整理可得一個關(guān)于X的4次多項式。假定在顫振狀態(tài)下X總為實數(shù),可以得到實部和虛部2個方程為:
A4RX4+A3RX3+A2RX2+A1RX+A0R=0
A3lX3+A2lX2+A1lX+A0l=0
(8)
式(8)中各項系數(shù)值為:
(9)
式(8)為一元高次方程,可對其進(jìn)行分別求解,實部方程和虛部方程分別可求得4個解和3個解。舍去負(fù)解和不合理的解,繪出實部解XR和虛部解Xi隨折算速度Vr變化的曲線,它們的交點(XC,VrC)即代表臨界顫振狀態(tài),臨界顫振速度計算式為:
UC=BωhXCVrC
(10)
對于水翼算例,由求得的顫振導(dǎo)數(shù),結(jié)合式(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
利用高斯消元法對上式進(jìn)行求解:實部方程解為:-1.042 1;-0.612 5;1.045 7;0.611 3。虛部方程解為:0.728 3;-0.720 3 -0.055 1,其他折算速度下求解同理。實、虛部方程的解舍去負(fù)解并進(jìn)行二次多項式擬合,實、虛部解擬合曲線如圖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。
使用紐馬克-β法對第4節(jié)求得的結(jié)果進(jìn)行驗證。紐馬克-β法是求解非線性振動的一種方法,本文根據(jù)其原理,編寫UDF嵌入Fluent軟件中實現(xiàn)水翼的流固耦合計算。判斷水翼是否發(fā)生顫振的依據(jù)如圖6所示。
圖6 臨界顫振狀態(tài)判斷依據(jù)
給予水翼一個很小的位移,使水翼處于不平衡狀態(tài),待流場穩(wěn)定后釋放,水翼若在自激力的作用下,發(fā)生振幅越來越大的振動,則水翼進(jìn)入顫振狀態(tài)。編寫Profile文件使水翼在前一秒內(nèi)兩自由度皆做頻率為0.75 Hz的小幅正弦運動,待1 s時,豎彎和扭轉(zhuǎn)2個自由度的振幅皆達(dá)到最值時,將處于不平衡狀態(tài)的水翼釋放,啟動紐馬克-β算法進(jìn)行計算。
圖7—圖11給出了水翼在各流速下的豎向振動及扭轉(zhuǎn)角時程圖。
圖7 扭轉(zhuǎn)、豎彎振動時程圖(U=3.0 m/s)
圖8 扭轉(zhuǎn)、豎彎振動時程圖(U=7.0 m/s)
圖9 扭轉(zhuǎn)、豎彎振動時程圖(U=7.3 m/s)
圖11 扭轉(zhuǎn)、豎彎振動時程圖(U=8.0 m/s)
由圖7—圖11可知,前1 s內(nèi)水翼豎彎及扭轉(zhuǎn)方向皆由Profile驅(qū)動做規(guī)律的正弦運動。在流速為3 m/s時,水翼的振幅迅速衰減,而后做微幅振動以抵消流場輸入的能量。當(dāng)流速提升至7 m/s時,水翼被釋放瞬間豎彎及扭轉(zhuǎn)方向振幅明顯要比U=3 m/s要大,振幅衰減速度比U=3 m/s時要慢。流速為7.3 m/s時,水翼豎彎及扭轉(zhuǎn)方向皆做等幅簡諧振動。當(dāng)流速提升至7.5 m/s,2自由度方向的振動總體呈現(xiàn)緩慢發(fā)散態(tài)勢。U=8 m/s時,發(fā)散速度要比U=7.5 m/s時快的多,結(jié)構(gòu)已經(jīng)發(fā)生顫振,在計算中由于振幅發(fā)散過快,水翼運動至動網(wǎng)格區(qū)域邊緣導(dǎo)致計算報錯。根據(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ī)律與之相似??梢姰?dāng)結(jié)構(gòu)進(jìn)入顫振狀態(tài)時,受到的升力及力矩會隨時間推移而越來越大,相應(yīng)的振動振幅也會越來越大,最終導(dǎo)致結(jié)構(gòu)失穩(wěn)。
圖15—圖17給出了水翼在流速為7、7.3以及7.5 m/s時的扭轉(zhuǎn)、豎彎振動相軌線圖。
圖12 升力及力矩時程圖(U=7.0 m/s)
圖13 升力及力矩時程圖(U=7.3 m/s)
圖14 升力及力矩時程圖(U=7.5 m/s)
圖15 扭轉(zhuǎn)、豎彎振動相軌線圖(U=7.0 m/s)
圖16 扭轉(zhuǎn)、豎彎振動相軌線圖(U=7.3 m/s)
圖17 扭轉(zhuǎn)、豎彎振動相軌線圖(U=7.5 m/s)
由圖15—圖17可知,在U=7.0 m/s時,水翼沒有達(dá)到臨界顫振速度,水翼的動能及勢能在運動過程中逐漸被消散,水翼在豎彎及扭轉(zhuǎn)自由度上的振幅及振動速度都在不斷減小,相軌線不斷向系統(tǒng)平衡點接近,直至達(dá)到穩(wěn)定狀態(tài)。當(dāng)U=7.3 m/s時,豎彎及扭轉(zhuǎn)振動演化為典型的極限環(huán)振動,水翼的能量在動能和勢能之間不斷轉(zhuǎn)化。當(dāng)流速提升至7.5 m/s時,水翼已達(dá)到顫振速度,系統(tǒng)的振幅及振動速度在流場不斷做正功的原因下不斷增大,直至結(jié)構(gòu)破壞。
對水翼在2個自由度上的振動頻率進(jìn)行了統(tǒng)計,如圖18所示。
圖18 扭轉(zhuǎn)、豎彎振動頻率
由圖18可知,2自由度的振動頻率都隨著流速的增加而逐漸減小。流速較低時,2自由度的振動頻率差距較大,隨著流速增加,振動頻率趨向一致;當(dāng)流速大于或等于臨界顫振速度時,2自由度的振動頻率相等。印證了顫振是由具有2個自由度以上的結(jié)構(gòu)物以同一頻率耦合振動的現(xiàn)象。
經(jīng)過以上分析,由紐馬克-β法求得的水翼臨界顫振速度為7.3 m/s,與Scanlan顫振導(dǎo)數(shù)理論求得的結(jié)果基本一致。
1)引入Scanlan顫振導(dǎo)數(shù)理論,由該理論成功求解水翼的臨界顫振速度并進(jìn)行了驗證,證明Scanlan顫振導(dǎo)數(shù)理論可以應(yīng)用在水翼的臨界顫振狀態(tài)求解。
2) 由紐馬克-β法對水翼的振動狀態(tài)進(jìn)行了分析,水翼一旦進(jìn)入顫振狀態(tài),振幅會不斷擴大直至破壞。升力和力矩的變化情況和振幅一致。
3) 水翼進(jìn)入臨界顫振狀態(tài)時,動能和勢能總和不變,并在不斷轉(zhuǎn)化中。2自由度的頻率一致,印證了顫振是由具有2個自由度以上的結(jié)構(gòu)物以同一頻率耦合振動的現(xiàn)象。