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

(1.中國船舶科學(xué)研究中心國家船舶噪聲與振動(dòng)重點(diǎn)實(shí)驗(yàn)室，江蘇 無錫 214000；2.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院， 南京 210000)

1 引言

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

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

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

2 水翼計(jì)算模型

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

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

圖1 2個(gè)自由度模型

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

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

圖2中，A區(qū)域?yàn)閯傂赃\(yùn)動(dòng)區(qū)域，該區(qū)域與水翼一起運(yùn)動(dòng)。該區(qū)域的半徑為b，采用四邊形結(jié)構(gòu)化網(wǎng)格劃分，以保證水翼運(yùn)動(dòng)時(shí)計(jì)算的收斂性。B區(qū)域?yàn)閯?dòng)網(wǎng)格區(qū)域，采用了三角形非結(jié)構(gòu)化網(wǎng)格進(jìn)行劃分。該區(qū)域的網(wǎng)格更新策略為“Smoothing”和“Remeshing”。C區(qū)域?yàn)橥饬鲌鰠^(qū)域，流場的尺寸為20b×26b。流場設(shè)置的較大，可以避免流場邊界與水翼附近流場相互干涉而影響計(jì)算的準(zhǔn)確性。

圖3 流場網(wǎng)格

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

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

Scanlan顫振導(dǎo)數(shù)理論內(nèi)容較多，具體可參考文獻(xiàn)[10]，本文僅對(duì)關(guān)鍵點(diǎn)進(jìn)行介紹。

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

(1)

(2)

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

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

(3)

(4)

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

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

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

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

(5)

將Scanlan自激力式(3)和式(4)代入2個(gè)自由度截面運(yùn)動(dòng)方程(5)，并引入無量綱時(shí)間概念：s=tU/B，對(duì)截面運(yùn)動(dòng)方程進(jìn)行無量綱化，可得：

(6)

(7)

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

A4RX4+A3RX3+A2RX2+A1RX+A0R=0

A3lX3+A2lX2+A1lX+A0l=0

(8)

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

(9)

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

UC=BωhXCVrC

(10)

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

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

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

圖5 實(shí)、虛部解擬合曲線

擬合曲線的交點(diǎn)為(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 驗(yàn)證與分析

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

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

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

圖7—圖11給出了水翼在各流速下的豎向振動(dòng)及扭轉(zhuǎn)角時(shí)程圖。

圖7 扭轉(zhuǎn)、豎彎振動(dòng)時(shí)程圖(U=3.0 m/s)

圖8 扭轉(zhuǎn)、豎彎振動(dòng)時(shí)程圖(U=7.0 m/s)

圖9 扭轉(zhuǎn)、豎彎振動(dòng)時(shí)程圖(U=7.3 m/s)

圖11 扭轉(zhuǎn)、豎彎振動(dòng)時(shí)程圖(U=8.0 m/s)

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

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

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

圖15—圖17給出了水翼在流速為7、7.3以及7.5 m/s時(shí)的扭轉(zhuǎn)、豎彎振動(dòng)相軌線圖。

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

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

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

圖15 扭轉(zhuǎn)、豎彎振動(dòng)相軌線圖(U=7.0 m/s)

圖16 扭轉(zhuǎn)、豎彎振動(dòng)相軌線圖(U=7.3 m/s)

圖17 扭轉(zhuǎn)、豎彎振動(dòng)相軌線圖(U=7.5 m/s)

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

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

圖18 扭轉(zhuǎn)、豎彎振動(dòng)頻率

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

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

6 結(jié)論

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

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

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