陳濤， 徐敏， 謝丹， 安效民

西北工業(yè)大學(xué) 航天學(xué)院， 西安 710072

低速流場中柔性懸臂板的后顫振響應(yīng)

陳濤， 徐敏*， 謝丹， 安效民

西北工業(yè)大學(xué) 航天學(xué)院， 西安 710072

建立了一個(gè)新的非線性氣動(dòng)彈性模型,對低速流場中柔性懸臂板的后顫振響應(yīng)特性進(jìn)行了分析。建模中考慮了結(jié)構(gòu)幾何非線性、氣動(dòng)力非線性以及兩者之間的強(qiáng)耦合效應(yīng)。通過實(shí)驗(yàn)數(shù)據(jù)對所建立的氣動(dòng)彈性模型進(jìn)行了驗(yàn)證。發(fā)現(xiàn)在低速流場中柔性懸臂板可能會(huì)以周期加倍的方式進(jìn)入混沌運(yùn)動(dòng)。結(jié)構(gòu)幾何非線性效應(yīng)和翼尖渦引起的非定常氣動(dòng)力效應(yīng)對柔性懸臂板的結(jié)構(gòu)響應(yīng)有顯著影響，而尾渦變形引起的非定常氣動(dòng)力對結(jié)構(gòu)運(yùn)動(dòng)的影響較小。還研究了不同耦合算法的差異，給出了小展弦比大柔性結(jié)構(gòu)非線性氣動(dòng)彈性數(shù)值仿真時(shí)耦合策略的選擇依據(jù)。

低速流場； 柔性懸臂板； 后顫振； 混沌運(yùn)動(dòng)； 強(qiáng)耦合

低速流場中懸臂板的氣動(dòng)彈性行為可以用來解釋現(xiàn)實(shí)中的許多現(xiàn)象，如機(jī)翼顫振、旗幟飄動(dòng)和打鼾等。最近， Dunnmon[1]和Doaré[2]等甚至開始研究利用柔性懸壁板的顫振來收集能量用于發(fā)電。因此，低速流場中懸臂板的非線性氣動(dòng)彈性特性研究是一項(xiàng)很有意義的工作。

關(guān)于低速流場中的非線性氣動(dòng)彈性問題，早期的研究可以追溯到二維情形，Lee等[3]對此做過詳細(xì)的綜述，二維問題中極限環(huán)振蕩(LCO)、混沌運(yùn)動(dòng)等現(xiàn)象都已通過數(shù)值計(jì)算或?qū)嶒?yàn)等手段發(fā)現(xiàn)。對于三維問題，非線性效應(yīng)比二維情況要復(fù)雜的多。Tang[4]和Attar[5]等通過耦合非定常渦格法和計(jì)算結(jié)構(gòu)動(dòng)力學(xué)的方法建立了一個(gè)非線性氣動(dòng)彈性模型，并研究了切尖三角翼和小展弦比長方形機(jī)翼在低速流場中的極限環(huán)振蕩特性。研究發(fā)現(xiàn)這些非線性氣動(dòng)彈性現(xiàn)象主要受結(jié)構(gòu)幾何非線性效應(yīng)的主導(dǎo)。Cavallaro等[6]采用有限元和勢流理論方法研究了柔性連接翼的非線性氣動(dòng)彈性特性。Arena等[7]采用勢流理論方法研究了大柔性機(jī)翼的后顫振響應(yīng)特性，他們發(fā)現(xiàn)流場的非定常效應(yīng)會(huì)對機(jī)翼運(yùn)動(dòng)產(chǎn)生重要影響。此外，Relvas和Suleman[8]基于結(jié)構(gòu)有限元方法和非定常渦格法也建立了一個(gè)非線性氣動(dòng)彈性模型并研究了柔性懸臂板的極限環(huán)振蕩現(xiàn)象。國內(nèi)尹維龍和田東奎[9]基于勢流理論方法研究了二維柔性翼型的氣動(dòng)彈性特性。謝長川等[10]基于勢流理論和結(jié)構(gòu)有限元方法研究了柔性機(jī)翼的低速顫振特性，結(jié)果發(fā)現(xiàn)結(jié)構(gòu)幾何非線性效應(yīng)對柔性翼的顫振邊界有重要影響。張健和向錦武[11]基于ONERA失速氣動(dòng)力模型和非線性梁模型研究了柔性飛機(jī)的非線性氣動(dòng)彈性特性。

以上研究表明，低速流場中三維柔性結(jié)構(gòu)的氣動(dòng)彈性建模需要考慮3種非線性效應(yīng)：① 由于結(jié)構(gòu)大變形引起的結(jié)構(gòu)幾何非線性效應(yīng)；② 由于流場渦等引起的非定常氣動(dòng)力效應(yīng)；③ 這兩者的耦合效應(yīng)。首先，對于結(jié)構(gòu)幾何非線性問題，Gordnier和Visbal[12]研究表明，氣動(dòng)彈性建模中遇到的主要問題是像von Karmann板理論等方法在求解幾何大變形問題時(shí)會(huì)產(chǎn)生“鎖住”效應(yīng)。其次，對于低速氣動(dòng)彈性問題，目前國際上廣泛采用勢流理論方法進(jìn)行非定常氣動(dòng)力建模。雖然有更精確的方法如CFD等，但在低速問題中柔性結(jié)構(gòu)響應(yīng)周期往往較長，采用CFD方法會(huì)消耗大量計(jì)算時(shí)間。Murua等[13]的研究表明非定常渦格法在計(jì)算準(zhǔn)確性和計(jì)算效率2個(gè)方面都能達(dá)到較滿意的效果，國內(nèi)許多學(xué)者也采用非定常渦格法研究低速飛行器的非定常氣動(dòng)力建模問題[14-15]，該方法還常用于氣動(dòng)力非定常效應(yīng)很強(qiáng)的撲翼氣動(dòng)彈性現(xiàn)象研究[16-17]。最后，關(guān)于結(jié)構(gòu)幾何非線性與非定常氣動(dòng)力的耦合效應(yīng)，Piperno和Farhat[18]的研究表明，氣動(dòng)彈性建模中采用松耦合策略會(huì)產(chǎn)生一定的數(shù)值誤差。對于柔性翼這樣非線性特性很強(qiáng)的氣動(dòng)彈性問題，松耦合策略可能會(huì)引起較大的誤差。

如前所述，Attar和Relvas等只研究了柔性壁板的周期性極限環(huán)振蕩現(xiàn)象，雖然他們都考慮了結(jié)構(gòu)和氣動(dòng)力的非線性效應(yīng)，但對于耦合策略的影響仍然有待研究。本文的工作是對Attar[5]和Relvas[8]等研究工作的進(jìn)一步擴(kuò)展。在本文的氣動(dòng)彈性建模中不僅詳細(xì)考慮了結(jié)構(gòu)幾何非線性效應(yīng)帶來的“鎖住”問題和非定常氣動(dòng)力問題，還考慮了結(jié)構(gòu)與氣動(dòng)力之間的強(qiáng)耦合問題。其中，結(jié)構(gòu)建模采用一種多變量有限元方法，非定常氣動(dòng)力建模采用非定常渦格法，采用強(qiáng)耦合策略實(shí)現(xiàn)結(jié)構(gòu)和氣動(dòng)力求解器的耦合。本文的結(jié)構(gòu)安排如下，首先給出新的非線性氣動(dòng)彈性建模過程，然后通過柔性三角翼極限環(huán)振蕩的仿真對所提出的氣動(dòng)彈性模型進(jìn)行驗(yàn)證，最后對柔性懸臂板后顫振響應(yīng)以及相關(guān)的非線性影響因素給出分析結(jié)果。

1 結(jié)構(gòu)幾何非線性建模

柔性懸臂板的結(jié)構(gòu)建模采用多變量實(shí)體殼單元。這種單元的有限元列式過程簡要介紹如下。

實(shí)體殼單元是一種多變量的有限元單元，這種單元的位移場、應(yīng)變場和應(yīng)力場是相互獨(dú)立的，可以分別假設(shè)獨(dú)立的插值模式。其有限元列式基于三場FHW變分原理[19]：

(1)

對于上述處理方法，筆者早期研究工作已經(jīng)證明這是一種解決“鎖住”效應(yīng)的有效手段[20]。當(dāng)給定了位移、應(yīng)變和應(yīng)力的插值形函數(shù)以后，就可以對式(1)進(jìn)行離散。本文的算法中采用一種無條件穩(wěn)定的能量-動(dòng)量守恒算法作為時(shí)間推進(jìn)算法，采用Newton-Raphson算法對式(1)進(jìn)行線性化：

(2)

(3)

2 非定常氣動(dòng)力建模

與Attar[5]和Relvas[8]等相同，本文采用非定常渦格法進(jìn)行非定常氣動(dòng)力建模。該方法假設(shè)結(jié)構(gòu)周圍的流動(dòng)是無旋附著流，并且不考慮黏性效應(yīng)。對于本文的懸臂板問題，非定常渦格法的非線性特性主要體現(xiàn)在以下3個(gè)方面：① 結(jié)構(gòu)表面的渦格跟隨結(jié)構(gòu)變形而非線性地改變形狀和位置。② 機(jī)翼尾部拖出的尾渦可以自由變形。③ 隨著結(jié)構(gòu)的變形，可能會(huì)有翼尖渦不斷拖出并卷起。在本文的研究中，翼尖渦的拖出方式與尾渦的拖出方式相同。這種翼尖渦模型在三角翼的相關(guān)研究中已被證明是一種有效的翼尖渦模擬方法[21]。在后面的數(shù)值仿真中會(huì)發(fā)現(xiàn)，翼尖渦效應(yīng)其實(shí)對結(jié)構(gòu)的氣動(dòng)彈性響應(yīng)有很大影響。但是在Relvas等[8]的研究中并沒有考慮翼尖渦效應(yīng)，因此，本文是對其研究工作的拓展。下面簡要介紹非定常渦格法的求解流程。

定常渦格法的求解過程主要是為了得到結(jié)構(gòu)表面的渦量強(qiáng)度分布，這些渦量強(qiáng)度用來計(jì)算結(jié)構(gòu)表面的氣動(dòng)力載荷。首先，將結(jié)構(gòu)表面用一系列渦格離散，在每一個(gè)渦格中，渦的強(qiáng)度為常量。每一個(gè)渦格片段的渦強(qiáng)度由共享該片段的相鄰2個(gè)渦格所決定。通過Biot-Savart定律可以求解每一個(gè)渦格片段的誘導(dǎo)速度場。隨著時(shí)間的推移，結(jié)構(gòu)表面的渦會(huì)從結(jié)構(gòu)后緣拖出形成尾渦。剛拖出的尾渦強(qiáng)度與相鄰的結(jié)構(gòu)表面渦強(qiáng)度相同，當(dāng)結(jié)構(gòu)表面的渦拖出形成尾渦后，渦的強(qiáng)度保持不變。非定常渦格法通過滿足物面無穿透邊界條件計(jì)算得到物面的渦強(qiáng)分布：

[V∞+Vbound(r,t)+Vwake(r,t)]·n=0

AΓbound=Rhs

(4)

式中：Vbound和Vwake分別為結(jié)構(gòu)表面和尾渦中的誘導(dǎo)速度；r為位置向量；t為時(shí)間；n為結(jié)構(gòu)表面的單位法向量；A為氣動(dòng)力影響系數(shù)矩陣；Γbound為待求解的未知渦量向量；Rhs為殘差向量。當(dāng)物面渦量求得以后，結(jié)構(gòu)表面的空氣壓強(qiáng)可由Bernoulli方程求解得到：

(5)

式中：Φ為速度勢函數(shù)；p為壓強(qiáng)；V∞和p∞分別為遠(yuǎn)場處的速度和壓強(qiáng)。

結(jié)構(gòu)上下表面的壓強(qiáng)差表達(dá)式為

Δp=pu-pl=ρ∞[Vu·Vu-Vl·Vl+

(6)

式中：Δp為壓強(qiáng)差；pu和pl分別為結(jié)構(gòu)上下表面的壓強(qiáng)；Vu和Vl分別為結(jié)構(gòu)上下表面的速度；ru和rl分別為結(jié)構(gòu)上下表面的位置向量。

由壓強(qiáng)差便可求出對應(yīng)的氣動(dòng)力為

Fi=ΔpiAini

(7)

3 非線性氣動(dòng)彈性建模

基于上述非線性有限元求解方法和非定常氣動(dòng)力求解方法可以建立柔性壁板的氣動(dòng)彈性模型。下面簡要介紹建模中所采用的耦合策略和結(jié)構(gòu)/流場求解器之間的數(shù)據(jù)交換方法。

3.1 耦合策略

眾所周知，傳統(tǒng)的松耦合策略會(huì)產(chǎn)生時(shí)間滯后效應(yīng)。這會(huì)導(dǎo)致在每一個(gè)時(shí)間步計(jì)算得到的結(jié)構(gòu)/流體交界面的位置不準(zhǔn)確。然而，在引言中所提到的所有關(guān)于柔性懸臂板的研究都是采用松耦合算法。本文的氣動(dòng)彈性建模采用耦合精度更高的強(qiáng)耦合策略(見圖1)，在每一個(gè)時(shí)間步內(nèi)引入預(yù)估-校正方法來提高結(jié)構(gòu)/流體交界面位置的計(jì)算精度。圖1中UVLM表示非定常渦格法氣動(dòng)力求解器，F(xiàn)EM表示有限元求解器。

內(nèi)迭代的收斂標(biāo)準(zhǔn)給定為相鄰兩次求得的結(jié)構(gòu)位移向量的相對誤差小于1%。此外，在實(shí)際的數(shù)值算法實(shí)現(xiàn)中還應(yīng)給定最大內(nèi)迭代步數(shù)，以防止出現(xiàn)永遠(yuǎn)不收斂的情況，不過在本文的研究中沒有出現(xiàn)這種不收斂的問題。關(guān)于上述的強(qiáng)耦合算法需要指出的是當(dāng)最大內(nèi)迭代步數(shù)給定為1時(shí)，強(qiáng)耦合算法將自動(dòng)退化為松耦合算法。本文強(qiáng)耦合算法的最大內(nèi)迭代步數(shù)給定為10，實(shí)際中發(fā)現(xiàn)所有內(nèi)迭代計(jì)算均會(huì)在10步以內(nèi)達(dá)到收斂。

3.2 數(shù)據(jù)交換方法

在非線性氣動(dòng)彈性建模中，結(jié)構(gòu)和流體求解器之間的數(shù)據(jù)交換方法不僅影響計(jì)算速度，更影響計(jì)算精度，特別是當(dāng)結(jié)構(gòu)變形很大時(shí)(例如極限環(huán)振蕩等)，這一問題更加突出。本文采用徑向基函數(shù)方法實(shí)現(xiàn)氣動(dòng)網(wǎng)格的變形[20]，這一方法非常適合于大變形過程的動(dòng)網(wǎng)格插值。氣動(dòng)力的插值與結(jié)構(gòu)位移場插值不同，它需要保證插值前后力場的等效，即總的力和力矩需要保持守恒。本文采用有限元形狀函數(shù)插值方法實(shí)現(xiàn)氣動(dòng)力的插值。

4 數(shù)值計(jì)算

本節(jié)首先對上述所建立的非線性氣動(dòng)彈性模型進(jìn)行驗(yàn)證，在此基礎(chǔ)上對柔性懸臂板的后顫振特性進(jìn)行數(shù)值仿真并給出分析結(jié)果。

4.1 非線性氣動(dòng)彈性模型的驗(yàn)證

通過一個(gè)三角翼極限環(huán)振蕩(LCO)的數(shù)值仿真與實(shí)驗(yàn)結(jié)果的對比對所提出的非線性氣動(dòng)彈性模型進(jìn)行驗(yàn)證。該模型是一個(gè)45° 的直角三角形平板，厚度為1.6 mm，翼根中部60%的部位受到固支約束。彈性模量為3.3 GPa，密度為1 299 kg/m3[6]。圖2為該三角翼模型的示意圖，V∞為來流速度。

經(jīng)過網(wǎng)格收斂性分析后，氣動(dòng)彈性仿真的結(jié)構(gòu)有限元網(wǎng)格共有819個(gè)實(shí)體殼單元。氣動(dòng)網(wǎng)格的大小為60×20(弦向×展向)。尾渦的長度取為弦長的4倍。

圖3給出了在1° 迎角時(shí)機(jī)翼翼尖面外方向極限環(huán)振蕩的最大速度對比結(jié)果。計(jì)算所得該三角翼的顫振速度為23.0 m/s，實(shí)驗(yàn)值為24.0 m/s。

圖3表明在來流速度小于26.0 m/s時(shí)，本文的計(jì)算結(jié)果比Attar等[5]的計(jì)算結(jié)果要大，兩者之間的差別隨著流速的增大而逐漸減小。在較高的來流速度時(shí)，計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果非常接近，并且與Attar等[5]采用高精度有限元方法的計(jì)算結(jié)果也非常接近。在較低的速度范圍時(shí)，計(jì)算的顫振點(diǎn)和極限環(huán)振蕩結(jié)果與Cavallaro等[6]的結(jié)果比較接近。需要指出的是在較低的來流速度范圍時(shí)，本文計(jì)算結(jié)果與其他兩個(gè)計(jì)算結(jié)果的差別主要是由于對結(jié)構(gòu)阻尼的不同處理導(dǎo)致的。Attar[5]和Cavallaro[6]等分別采用了不同的結(jié)構(gòu)阻尼形式，由于非線性氣動(dòng)彈性計(jì)算中結(jié)構(gòu)阻尼的影響是一個(gè)比較復(fù)雜的問題，目前國際上仍沒有定論，因此本文沒有考慮結(jié)構(gòu)阻尼。從計(jì)算結(jié)果的對比來看，本文所提出的非線性氣動(dòng)彈性模型能較好地模擬大幅值的極限環(huán)運(yùn)動(dòng)。

需要指出的圖3中本文的計(jì)算結(jié)果采用的是強(qiáng)耦合策略，這主要是為了驗(yàn)證前文所建立的強(qiáng)耦合策略的正確性。強(qiáng)耦合策略由于帶有內(nèi)迭代過程而比較耗時(shí)，對于這里的三角翼模型，傳統(tǒng)松耦合策略與強(qiáng)耦合策略的計(jì)算結(jié)果也比較接近，見圖4。

由于本文的氣動(dòng)彈性建模借鑒并改進(jìn)了Relvas等[8]的方法，因此有必要將本文方法與其進(jìn)行對比。文獻(xiàn)[8]中研究的方形懸臂板模型展長(弦長)為0.3 m，板厚度為1 mm，材料為鋁合金，計(jì)算初始迎角為0°。圖5給出了采用Nastran線性分析方法得到的顫振V-g圖，圖6給出了本文的計(jì)算結(jié)果。由圖5可知，Nastran計(jì)算的顫振速度為53.7 m/s，由圖6可知，本文的計(jì)算結(jié)果為57.5 m/s。文獻(xiàn)[8]的計(jì)算結(jié)果為42.0 m/s，對比可知，本文結(jié)果與Nastran的相對誤差為7.08%，文獻(xiàn)[8]的結(jié)果與Nastran的相對誤差為21.79%。計(jì)算結(jié)果表明，本文的方法用于方形懸臂板后顫振響應(yīng)分析是可行的。

4.2 柔性懸臂板的非線性氣動(dòng)彈性特性

經(jīng)過實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證以后，將所建立的氣動(dòng)彈性模型用于柔性懸臂板(圖7)的后顫振響應(yīng)分析。板的厚度為0.5 mm，材料屬性如下：彈性模量為71 GPa，密度為2 800 kg/m3，泊松比為 0.33。需要說明的是文獻(xiàn)[8]中的方形懸臂板模型偏硬，不適合于較大幅值的極限環(huán)振蕩研究，而本文第2個(gè)算例主要是為了研究較大幅度的極限環(huán)振蕩，因此本文研究的懸臂板模型尺寸在文獻(xiàn)[8]的基礎(chǔ)上進(jìn)行了調(diào)整。在數(shù)值仿真之前，通過網(wǎng)格收斂性分析確定采用的計(jì)算網(wǎng)格。最終使用的結(jié)構(gòu)網(wǎng)格數(shù)量為30×30(弦向×展向)，對應(yīng)的氣動(dòng)力計(jì)算網(wǎng)格數(shù)量為40×40。計(jì)算過程中尾渦的長度為弦長的4倍。

對柔性懸臂板進(jìn)行氣動(dòng)彈性數(shù)值仿真，初始迎角為0°。表1給出了本文計(jì)算的結(jié)構(gòu)固有頻率和解析解[22](雙向梁函數(shù)組合級數(shù)逼近法)以及Nastran計(jì)算結(jié)果的對比。圖8給出了Nastran的顫振分析結(jié)果。圖9給出了來流速度分別為38.3、46、54、58、62、66 m/s時(shí)前緣翼尖的結(jié)構(gòu)位移響應(yīng)。

由表1可知，本文采用實(shí)體殼單元計(jì)算的結(jié)構(gòu)固有頻率與解析解和Nastran的計(jì)算結(jié)果非常接近，與解析解的最大誤差僅為1.49%(二階)。由圖8可知，Nastran分析的結(jié)構(gòu)顫振邊界為35.7 m/s，顫振發(fā)生的機(jī)理為二階扭轉(zhuǎn)模態(tài)發(fā)散。圖9(a)表明本文基于有限元和渦格法計(jì)算的顫振邊界為38.3 m/s，兩者之間的相對誤差為6.79%。

Method1st2nd3rd4th5thAnalysisresult10.53326.04464.56782.32194.088Nastran10.59925.61364.61082.70293.531Currentwork10.60525.65564.84583.11593.946

結(jié)構(gòu)響應(yīng)曲線圖9表明，當(dāng)來流速度高于顫振速度時(shí)，結(jié)構(gòu)會(huì)發(fā)生極限環(huán)振蕩。當(dāng)來流速度高于54 m/s時(shí)，結(jié)構(gòu)響應(yīng)變得非常復(fù)雜。為了更好地分析結(jié)構(gòu)響應(yīng)隨來流的變化趨勢，圖10給出了相應(yīng)的相平面圖。

相平面圖表明，當(dāng)來流速度為46 m/s時(shí)，結(jié)構(gòu)響應(yīng)形式為周期性的極限環(huán)振蕩。當(dāng)來流速度為54 m/s時(shí)，結(jié)構(gòu)響應(yīng)形式為準(zhǔn)周期的極限環(huán)振蕩。當(dāng)來流速度為58 m/s時(shí)，結(jié)構(gòu)響應(yīng)出現(xiàn)周期加倍現(xiàn)象。當(dāng)來流速度增大到62 m/s和66 m/s時(shí)，結(jié)構(gòu)響應(yīng)的規(guī)律變得非常復(fù)雜。筆者早期關(guān)于超聲速流中壁板顫振的研究表明這種運(yùn)動(dòng)形式的變化是一種進(jìn)入混沌運(yùn)動(dòng)的路徑[23]。因此，本文首次發(fā)現(xiàn)了低速流場中三維柔性懸臂板很有可能也會(huì)以周期加倍的方式進(jìn)入混沌運(yùn)動(dòng)。不過由于本文計(jì)算得到的時(shí)域響應(yīng)結(jié)果步數(shù)有限，尚不能斷定一定發(fā)生了混沌運(yùn)動(dòng)，定量的證明將在后續(xù)工作中采用降階的方法進(jìn)行研究。

從位移響應(yīng)曲線(圖9)還可以看出，在較大來流速度下，結(jié)構(gòu)首先經(jīng)歷一個(gè)較大幅值的瞬態(tài)響應(yīng)過程，然后收斂到穩(wěn)態(tài)響應(yīng)。需要指出的是，本文在氣動(dòng)彈性建模中沒有考慮結(jié)構(gòu)阻尼。另外，根據(jù)國際權(quán)威學(xué)者Hallissy和Cesnik[24]的研究，勢流理論與高精度的氣動(dòng)力計(jì)算方法(如CFD方法等)相比氣動(dòng)阻尼偏小，這可能導(dǎo)致結(jié)構(gòu)的瞬態(tài)響應(yīng)幅值偏大。在Cavallaro等[6]關(guān)于柔性連接翼非線性氣動(dòng)彈性特性的研究中也發(fā)現(xiàn)了類似的瞬態(tài)響應(yīng)幅值很大的情況，可見這是采用渦格法計(jì)算時(shí)的一種普遍現(xiàn)象。由于本文主要研究的是結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)特性，因此有必要考察結(jié)構(gòu)瞬態(tài)響應(yīng)對穩(wěn)態(tài)響應(yīng)的影響。

圖11給出了分別采用結(jié)構(gòu)一階模態(tài)和二階模態(tài)作為有限元計(jì)算的初始速度擾動(dòng)時(shí)結(jié)構(gòu)位移響應(yīng)的對比。由圖可知，結(jié)構(gòu)在不同初始擾動(dòng)激勵(lì)下經(jīng)歷了不同的瞬態(tài)響應(yīng)過程，但最終的穩(wěn)態(tài)響應(yīng)非常接近。因此說明結(jié)構(gòu)瞬態(tài)響應(yīng)不會(huì)對穩(wěn)態(tài)響應(yīng)造成很大影響。

以上研究表明低速流場中柔性懸臂板的后顫振運(yùn)動(dòng)存在很強(qiáng)的非線性特性，因此有必要對其中的非線性因素進(jìn)行深入研究。下面將從結(jié)構(gòu)幾何非線性、氣動(dòng)力非線性以及耦合方式的影響3個(gè)方面進(jìn)行分析。

圖12給出了懸臂板在來流速度為46 m/s時(shí)考慮和不考慮結(jié)構(gòu)幾何非線性效應(yīng)的計(jì)算結(jié)果對比，由圖可知在懸臂板的后顫振響應(yīng)中結(jié)構(gòu)幾何非線性效應(yīng)的影響非常大。

為了進(jìn)一步分析結(jié)構(gòu)幾何非線性效應(yīng)對后顫振響應(yīng)的影響，表2給出了在穩(wěn)態(tài)響應(yīng)階段結(jié)構(gòu)最大位移處計(jì)算的模態(tài)頻率?？梢钥闯?，隨著來流速度的增加，結(jié)構(gòu)在最大位移處的一階模態(tài)頻率變化較小，二階及高階模態(tài)頻率變化較大。在來流速度為54 m/s時(shí)，結(jié)構(gòu)極限環(huán)振蕩最大位移達(dá)到了板厚度的4倍左右，此時(shí)結(jié)構(gòu)二階模態(tài)頻率增大了約10%。圖12中采用線性結(jié)構(gòu)模型的時(shí)域響應(yīng)在t=0.3 s時(shí)結(jié)構(gòu)最大位移已達(dá)到板厚度的5倍左右，此時(shí)結(jié)構(gòu)幾何非線性效應(yīng)已經(jīng)非常明顯。

幾何非線性效應(yīng)主要是由結(jié)構(gòu)內(nèi)部應(yīng)力對結(jié)構(gòu)剛度的影響造成的，圖13給出了表2中各來流速度對應(yīng)的結(jié)構(gòu)Von Mises應(yīng)力分布。由圖可知，當(dāng)來流速度從38.3 m/s增大到46 m/s時(shí)，結(jié)構(gòu)內(nèi)部應(yīng)力最大值由0.06 MPa增大到了2.4 MPa，此時(shí)結(jié)構(gòu)剛度產(chǎn)生了較大變化，幾何非線性效應(yīng)已非常明顯。

除了結(jié)構(gòu)幾何非線性，氣動(dòng)力非線性可能也會(huì)對柔性懸臂板的后顫振響應(yīng)產(chǎn)生重要影響。下面主要考察尾渦和翼尖渦對非線性氣動(dòng)力的貢獻(xiàn)。圖14給出了考慮翼尖渦和不考慮翼尖渦的計(jì)算結(jié)果對比。圖15給出了在來流速度為66 m/s不考慮翼尖渦時(shí)結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)的變形情況。圖16 給出了分別采用自由變形尾渦模型和平面尾渦模型的計(jì)算結(jié)果。需要指出的是，此時(shí)為了單獨(dú)考察尾渦變形的影響，數(shù)值計(jì)算時(shí)沒有考慮翼尖渦模型。

當(dāng)不考慮翼尖渦時(shí)(圖14虛線)，結(jié)構(gòu)的穩(wěn)態(tài)響應(yīng)形式為動(dòng)態(tài)屈曲。雖然初始迎角為0°，但在瞬態(tài)響應(yīng)過程中，結(jié)構(gòu)產(chǎn)生了扭轉(zhuǎn)變形(圖15)，翼面上當(dāng)?shù)赜遣粸?°，所產(chǎn)生的升力使得結(jié)構(gòu)最終在某一平衡位置附近做幅值非常小的動(dòng)態(tài)屈曲運(yùn)動(dòng)，這種現(xiàn)象在筆者早期關(guān)于超聲速氣流中懸臂板的混沌運(yùn)動(dòng)研究中發(fā)現(xiàn)過(見文獻(xiàn)[23])。但考慮翼尖渦后，結(jié)構(gòu)響應(yīng)形式完全不同，圖14表明是否考慮翼尖渦會(huì)對計(jì)算結(jié)果產(chǎn)生重要影響。圖16表明是否考慮尾渦的變形對結(jié)構(gòu)響應(yīng)幾乎沒有影響。由以上分析可知，翼尖渦對非線性氣動(dòng)力的貢獻(xiàn)較大，而尾渦的自由變形引起的氣動(dòng)力非線性效應(yīng)較小。

表2 變形最大位移處結(jié)構(gòu)模態(tài)頻率

Table 2 Structural modal frequencies at position of maximum deformation Hz

Inflowvelocity/(m·s-1)1st2nd3rd4th5th010.60525.65564.84583.11593.94638.310.60625.66364.92983.19294.0564610.62626.98265.33083.54095.0025410.65428.25965.77484.80295.8795811.13441.06678.044108.315119.752

圖17給出了在來流速度為54 m/s和66 m/s時(shí)分別采用松耦合、改進(jìn)的松耦合以及強(qiáng)耦合策略的計(jì)算結(jié)果對比。這里的改進(jìn)松耦合策略為Relvas和Suleman[8]在研究相似柔性懸臂板非線性氣動(dòng)彈性問題時(shí)所采用的耦合方式。由圖可知，在較低的來流速度54 m/s時(shí)，改進(jìn)的松耦合方式與強(qiáng)耦合方式的計(jì)算結(jié)果非常接近，而松耦合算法與強(qiáng)耦合算法的差異隨著時(shí)間的推移會(huì)逐漸增大。在較高的來流速度66 m/s時(shí)，不管是松耦合還是改進(jìn)的松耦合方式與強(qiáng)耦合方式的差別都是非常大的。

仔細(xì)分析還發(fā)現(xiàn)，在數(shù)值仿真過程中前述柔性三角翼與柔性懸臂板在最大來流速度時(shí)結(jié)構(gòu)的位移響應(yīng)均達(dá)到了板厚度的20倍左右(圖4和圖9(f))，因此，2個(gè)模型的幾何非線性效應(yīng)都是很強(qiáng)的。但在柔性三角翼的研究中，當(dāng)來流速度較高時(shí)，不同的耦合策略計(jì)算得到的位移響應(yīng)差別不大。而對于柔性懸臂板模型，不同的耦合策略計(jì)算得到的位移響應(yīng)相差卻很大。2個(gè)模型的主要區(qū)別在于氣動(dòng)力非線性效應(yīng)的強(qiáng)弱，對于柔性三角翼，由于不存在翼尖渦引起的非線性氣動(dòng)力，因此，其氣動(dòng)力非線性效應(yīng)要弱于柔性懸臂板。據(jù)此可知，在小展弦比大柔性結(jié)構(gòu)的非線性氣動(dòng)彈性數(shù)值仿真時(shí)，當(dāng)只有結(jié)構(gòu)幾何非線性效應(yīng)占主導(dǎo)地位時(shí)，不同耦合策略的計(jì)算結(jié)果差別可能不是很大。當(dāng)結(jié)構(gòu)幾何非線性與氣動(dòng)力非線性效應(yīng)均較強(qiáng)時(shí)松耦合策略和改進(jìn)的松耦合策略都會(huì)引起很大的誤差。

5 結(jié) 論

建立一個(gè)新的非線性氣動(dòng)彈性模型對低速流場中柔性懸臂板的后顫振特性進(jìn)行了數(shù)值仿真。結(jié)構(gòu)幾何非線性問題采用多變量有限元方法求解，非定常氣動(dòng)力采用非定常渦格法求解，實(shí)現(xiàn)了結(jié)構(gòu)與氣動(dòng)力求解器之間的強(qiáng)耦合算法。采用柔性三角翼模型對所提出的氣動(dòng)彈性模型進(jìn)行了驗(yàn)證，結(jié)果表明，對于大柔性結(jié)構(gòu)的非線性氣動(dòng)彈性問題，所提出的算法可以給出較好的仿真結(jié)果。

首次發(fā)現(xiàn)了在低速流場中三維柔性懸臂板可能會(huì)以位移響應(yīng)周期加倍的形式進(jìn)入混沌運(yùn)動(dòng)。數(shù)值仿真結(jié)果表明，在柔性懸臂板的后顫振響應(yīng)中，結(jié)構(gòu)幾何非線性效應(yīng)和氣動(dòng)力非線性效應(yīng)對結(jié)構(gòu)運(yùn)動(dòng)都有很大影響。其中，氣動(dòng)力非線性效應(yīng)主要由翼尖渦引起，而尾渦的變形對非線性氣動(dòng)力的貢獻(xiàn)相對較小。在小展弦比大柔性結(jié)構(gòu)的非線性氣動(dòng)彈性數(shù)值仿真時(shí)，當(dāng)只有結(jié)構(gòu)幾何非線性效應(yīng)占主導(dǎo)地位時(shí)，松耦合策略與強(qiáng)耦合策略的計(jì)算結(jié)果可能差別不是很大，此時(shí)可以采用松耦合策略或改進(jìn)的松耦合策略。當(dāng)結(jié)構(gòu)幾何非線性與氣動(dòng)力非線性效應(yīng)均較強(qiáng)時(shí)必須采用強(qiáng)耦合策略。

需要說明的是本文目前只是初步研究了結(jié)構(gòu)幾何非線性、氣動(dòng)力非線性和耦合策略對柔性結(jié)構(gòu)非線性氣動(dòng)彈性仿真的影響。下一步的工作將對其他非線性因素如結(jié)構(gòu)阻尼效應(yīng)等進(jìn)行研究。另外，還將開展相關(guān)的風(fēng)洞試驗(yàn)研究。

[1] DUNNMON J A, STANTON S C, MANN B P, et al. Power extraction from aeroelastic limit cycle oscillations[J]. Journal of Fluids and Structures, 2011, 27(8): 1182-1198.

[2] DOARé O, MICHELIN S. Piezoelectric coupling in energy-harvesting fluttering flexible plates: Linear stability analysis and conversion efficiency[J]. Journal of Fluids and Structures, 2011, 27(8): 1357-1375.

[3] LEE B H K, PRICE S J, WONG Y S. Nonlinear aeroelastic analysis of airfoils: Bifurcation and chaos[J]. Progress in Aerospace Sciences, 1999, 35(3): 205-334.

[4] TANG D, DOWELL E H, HALL K C. Limit cycle oscillations of a cantilevered wing in low subsonic flow[J]. AIAA Journal, 1999, 37(3): 364-371.

[5] ATTAR P J, DOWELL E H, WHITE J R. Modeling the LCO of a delta wing using a high fidelity structural model：AIAA-2004-1692[R]. Reston: AIAA, 2004.

[6] CAVALLARO R, IANNELLI A, DEMASI L, et al. Phenomenology of nonlinear aeroelastic responses of highly deformable joined-wings configurations: AIAA-2014-1199[R]. Reston: AIAA, 2014.

[7] ARENA A, LACARBONARA W, MARZOCCA P. Nonlinear aeroelastic formulation and postflutter analysis of flexible high aspect ratio wings[J]. Journal of Aircraft, 2013, 50(6): 1748-1764.

[8] RELVAS A, SULEMAN A. Fluid-structure interaction modelling of nonlinear aeroelastic structures using the finite element corotational theory[J]. Journal of Fluids and Structures, 2006, 22(1): 59-75.

[9] 尹維龍, 田東奎. 柔性翼型的氣動(dòng)彈性建模與顫振特性分析[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào), 2012, 44(9): 69-72. YIN W L, TIAN D K. Aeroelastic modeling and flutter characteristics of flexible aerofoil[J]. Journal of Harbin Institute of Technology, 2012, 44(9): 69-72 (in Chinese).

[10] 謝長川, 吳志剛, 楊超. 大展弦比柔性機(jī)翼的氣動(dòng)彈性分析[J]. 北京航空航天大學(xué)學(xué)報(bào), 2003, 29(12): 1087-1090. XIE C C, WU Z G, YANG C. Aeroelastic analysis of flexible large aspect ratio wing[J]. Journal of Beijing University of Aeronautics and Astronautics, 2003, 29(12): 1087-1090 (in Chinese).

[11] 張健, 向錦武. 柔性飛機(jī)非線性氣動(dòng)彈性與飛行動(dòng)力學(xué)耦合靜、動(dòng)態(tài)特性[J]. 航空學(xué)報(bào), 2011, 32(9): 1569-1582. ZHANG J, XIANG J W. Static and dynamic characteristics of coupled nonlinear aeroelasticity and flight dynamics of flexible aircraft[J]. Acta Aeronautica et Astronautica Sinica, 2011, 32(9): 1569-1582 (in Chinese).

[12] GORDNIER R E, VISBAL M R. Development of a three-dimensional viscous aeroelastic solver for nonlinear panel flutter[J]. Journal of Fluids and Structures, 2002, 16(4): 497-527.

[13] MURUA J, PALACIOS R, GRAHAM J M R. Applications of the unsteady vortex-lattice method in aircraft aeroelasticity and flight dynamics[J]. Progress in Aerospace Sciences, 2012, 55(5): 46-72.

[14] 劉燚, 謝長川, 王立波, 等. 柔性飛機(jī)大變形曲面氣動(dòng)力計(jì)算及配平分析[J]. 工程力學(xué), 2015, 32(10): 239-249. LIU Y, XIE C C, WANG L B, et al. Nonplanar aerodynamic computation and trim analysis under large deflection of flexible aircraft[J]. Engineering Mechanics, 2015, 32(10):239-249 (in Chinese).

[15] 宋磊, 楊華, 解靜峰, 等. 基于改進(jìn)渦格法的飛翼布局飛機(jī)穩(wěn)定性導(dǎo)數(shù)計(jì)算[J]. 南京航空航天大學(xué)學(xué)報(bào), 2014, 46(3): 457-462. SONG L, YANG H, XIE J F, et al. Predicting stability derivatives of flying wing aircraft based on improved vortex lattice method[J]. Journal of Nanjing University of Aeronautics and Astronautics, 2014, 46(3): 457-462 (in Chinese).

[16] ROCCIA B A, PREIDIKMAN S, MASSA J C, et al. Modified unsteady vortex-lattice method to study flapping wings in hover flight[J]. AIAA Journal, 2013, 51(11): 2628-2642.

[17] 賀紅林, 周翔. 柔性撲翼非定常渦格法氣動(dòng)力計(jì)算的改進(jìn)與實(shí)現(xiàn)[J]. 航空學(xué)報(bào), 2010, 31(6): 1121-1126. HE H L, ZHOU X. Implementation of an improved unsteady vortex lattice method for flexible flapping-wing aerodynamic computation[J]. Acta Aeronautica et Astronautica Sinica, 2010, 31(6): 1121-1126 (in Chinese).

[18] PIPERNO S, FARHAT C. Partitioned procedures for the transient solution of coupled aeroelastic problems-Part II: Energy transfer analysis and three-dimensional applications[J]. Computer Methods in Applied Mechanics and Engineering, 2001, 190(24-25): 3147-3170.

[19] VU-QUOC L, TAN X G. Optimal solid shells for non-linear analysis of multilayer composites. 1. Static[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(9-10): 975-1016.

[20] CHEN T, XU M, XIE L. Aeroelastic modeling using geometrically nonlinear solid-shell element[J]. AIAA Journal, 2014, 52(9): 1980-1993.

[21] LEVIN D, KATZ J. A vortex-lattice method for the calculation of the nonsteady separated flow over delta wings: AIAA-1980-1803[R]. Reston: AIAA, 1980.

[22] 曹志遠(yuǎn). 板殼振動(dòng)理論[M]. 北京：中國鐵道出版社, 1989: 455-457. CAO Z Y. Vibration theory of plates and shells[M]. Beijing: China Railway Press, 1989: 455-457 (in Chinese).

[23] XIE D, XU M, DAI H, et al. Observation and evolution of chaos for a cantilever plate in supersonic flow[J]. Journal of Fluids and Structures, 2014, 50(6): 271-291.

[24] HALLISSY P B, CESNIK C E S. High-fidelity aeroelastic analysis of very flexible aircraft: AIAA-2011-1914[R]. Reston: AIAA, 2011.

(責(zé)任編輯：李明敏)

＊Corresponding author. E-mail: xumin@nwpu.edu.cn

Post flutter response of flexible cantilever plate inlow speed flow field

CHEN Tao, XU Min*, XIE Dan, AN Xiaomin

SchoolofAstronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China

A new nonlinear aeroelastic model is established for the analysis of post flutter behavior of the flexible cantilever plate in low speed flows. Geometric nonlinearity, aerodynamic nonlinearity as well as the strong coupling effect between geometric and aerodynamic nonlinearities are considered in the modeling. The aeroelastic model is verified with experimental data. The study shows that the flexible cantilever plate in low speed flows can enter into chaotic movement by means of periodic doubling in structural response. It is found that geometric nonlinearity and the unsteady aerodynamic effect caused by tip vortexes have significant impact on the structural response, while the unsteady aerodynamic force caused by deformation of wake vortexes has slight impact on the structural response. Different coupling strategies are also investigated to provide a reference for the selection of coupling strategy in nonlinear aeroelastic simulations of low-aspect-ratio flexible structures in low speed flows.

low speed flow field; flexible cantilever plate; post flutter; chaotic movement; strong coupling

2016-03-11； Revised：2016-10-10； Accepted：2016-11-14； Published online：2016-11-22 10:56

URL：www.cnki.net/kcms/detail/11.1929.V.20161122.1056.002.html

s：National Natural Science Foundation of China (11202165, 11402212)

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0296

2016-03-11； 退修日期：2016-10-10； 錄用日期：2016-11-14； 網(wǎng)絡(luò)出版時(shí)間：2016-11-22 10:56

www.cnki.net/kcms/detail/11.1929.V.20161122.1056.002.html

國家自然科學(xué)基金 (11202165, 11402212)

＊通訊作者.E-mail: xumin@nwpu.edu.cn

陳濤， 徐敏， 謝丹， 等． 低速流場中柔性懸臂板的后顫振響應(yīng)[J]． 航空學(xué)報(bào), 2017, 38(3): 120215. CHEN T, XU M, XIE D, et al. Post flutter response of flexible cantilever plate in low speed flow field[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(3): 120215.

V211.47

A

1000-6893(2017)03-120215-14