倪迎鴿, 侯赤, 萬小朋, 趙美英

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安　710072)

折疊機(jī)翼的參數(shù)化氣動彈性建模與顫振分析

倪迎鴿, 侯赤, 萬小朋, 趙美英

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安710072)

摘要：針對折疊機(jī)翼的特點(diǎn)建立了顫振分析的參數(shù)化氣動彈性模型。參數(shù)化的結(jié)構(gòu)模型基于模態(tài)綜合法實(shí)現(xiàn);參數(shù)化的氣動力模型采用偶極子網(wǎng)格法建立;并且闡述了基于Gram矩陣范數(shù)對于氣動彈性系統(tǒng)顫振邊界的預(yù)測方法。以折疊機(jī)翼完全展開和完全折疊構(gòu)型為例,將文中方法獲得的顫振邊界與特征值法獲得的結(jié)果進(jìn)行了對比,驗(yàn)證了該方法的正確性。通過對不同折疊角度下的顫振邊界的分析可以得出:顫振邊界對折疊角度很敏感;隨著折疊角度的增加,顫振模態(tài)發(fā)生了變化;較高的結(jié)構(gòu)模態(tài)阻尼比可以提高顫振速度,推遲顫振現(xiàn)象的發(fā)生。

關(guān)鍵詞：折疊機(jī)翼；Gram矩陣；參數(shù)化模型；氣動彈性；顫振

變體飛機(jī)可以根據(jù)不同的飛行任務(wù)改變構(gòu)型,因此,受到了極大關(guān)注[1]。針對變體飛機(jī)的顫振問題,許多學(xué)者開展了廣泛的研究[2-4]。與固定翼的顫振分析不同,折疊機(jī)翼因?yàn)檎郫B角的變化需要進(jìn)行不同構(gòu)型下的顫振分析。如果不采用參數(shù)化的建模方法,其顫振分析需要建立不同的結(jié)構(gòu)模型與氣動模型,顯然,該過程是比較低效的。Zhao等[5]提出了參數(shù)化的氣動彈性模型的建立方法,在MATLAB中實(shí)現(xiàn)了顫振分析。但在參數(shù)化過程中會遇到2個(gè)問題:①在不同的折疊角下,整個(gè)折疊機(jī)翼的質(zhì)量矩陣和剛度矩陣會變?yōu)椴粚ΨQ的稀疏矩陣,此時(shí)采用MATLAB中eigs函數(shù)求解模態(tài)振型會出現(xiàn)求解問題,因?yàn)閑igs函數(shù)要求質(zhì)量矩陣為對稱矩陣;若將稀疏矩陣轉(zhuǎn)化為全矩陣,采用eig函數(shù)求解,又會出現(xiàn)數(shù)據(jù)的存儲問題;②采用經(jīng)典的顫振分析時(shí),需要求解顫振特征方程的特征值,倘若對顫振特征值跟綜失敗,會出現(xiàn)“竄支”問題,從而得到錯(cuò)誤的結(jié)果[6]。

隨著以狀態(tài)空間模型為基礎(chǔ)的現(xiàn)代控制理論的發(fā)展,系統(tǒng)的穩(wěn)定性也具有了成熟的理論可以借鑒。當(dāng)氣動彈性方程表示成時(shí)域的狀態(tài)空間時(shí),根據(jù)線性系統(tǒng)理論,當(dāng)系統(tǒng)狀態(tài)矩陣的全部特征值均位于復(fù)平面的左半部時(shí),系統(tǒng)是穩(wěn)定的,此結(jié)論為顫振發(fā)生的判定提供了依據(jù)。但特征值法的缺點(diǎn)是需要對不同速度下的系統(tǒng)特征值進(jìn)行重復(fù)求解,當(dāng)系統(tǒng)的維度較高時(shí),特征值的求解比較費(fèi)時(shí),不便于工程應(yīng)用。在現(xiàn)代控制理論中,Gram 矩陣具有廣泛的應(yīng)用[7-9]。最近,Bueno等[10]提出了氣動彈性系統(tǒng)穩(wěn)定性分析的新方法,該方法以Gram矩陣范數(shù)為基礎(chǔ),從能量的觀點(diǎn)進(jìn)一步闡述了顫振產(chǎn)生的機(jī)理,同時(shí)避免了對高維氣動彈性系統(tǒng)特征值的求解,提高了計(jì)算效率。

本文采用模態(tài)綜合法建立了折疊機(jī)翼的參數(shù)化結(jié)構(gòu)模型,結(jié)合偶極子網(wǎng)格法建立了參數(shù)化氣動力模型,并構(gòu)建了氣動彈性狀態(tài)方程,采用Gram矩陣范數(shù)獲得了折疊機(jī)翼的顫振邊界;將折疊機(jī)翼完全展開和完全折疊2種構(gòu)型的顫振邊界與特征值法獲得的結(jié)果進(jìn)行了對比,驗(yàn)證了本文方法的正確性。本文方法不需要重復(fù)建立折疊機(jī)翼的結(jié)構(gòu)和氣動力模型,可以定性地確定顫振耦合的模態(tài)信息,提高了折疊機(jī)翼顫振分析的效率。

1參數(shù)化的氣動彈性模型

折疊機(jī)翼結(jié)構(gòu)由機(jī)身、內(nèi)翼和外翼3個(gè)子結(jié)構(gòu)組成,分別記為子結(jié)構(gòu)A、B、C。本文將各個(gè)子結(jié)構(gòu)簡化為等厚度的鋁板,在MSC.NASTRAN中采用CQUAD4單元來模擬,該單元的每個(gè)節(jié)點(diǎn)有6個(gè)自由度,即3個(gè)平動自由度ux、uy、uz以及3個(gè)轉(zhuǎn)動自由度θx、θy、θz。對于連接處的每個(gè)界面節(jié)點(diǎn),扭轉(zhuǎn)彈簧耦合θx,MPC耦合其余5個(gè)自由度,實(shí)現(xiàn)了各個(gè)子結(jié)構(gòu)的界面節(jié)點(diǎn)之間的鉸鏈模擬。其中,子結(jié)構(gòu)A與子結(jié)構(gòu)B之間的扭轉(zhuǎn)彈簧剛度為KA,子結(jié)構(gòu)B與子結(jié)構(gòu)C之間的扭轉(zhuǎn)彈簧剛度為KB。折疊機(jī)翼在局部坐標(biāo)下的有限元模型如圖1所示。

圖1　折疊機(jī)翼的有限元模型和坐標(biāo)系　(每個(gè)圓點(diǎn)代表鉸鏈所處位置)

1.1參數(shù)化的結(jié)構(gòu)模型

各個(gè)子結(jié)構(gòu)在局部坐標(biāo)系下的無阻尼的動力學(xué)方程如下所示:

(1)

當(dāng)采用圖1中的坐標(biāo)系時(shí),子結(jié)構(gòu)A與子結(jié)構(gòu)C在總體坐標(biāo)下與局部坐標(biāo)系的位移矢量相同,不需要進(jìn)行變換。對于子結(jié)構(gòu)B,其在總體坐標(biāo)系下的位移矢量與在局部坐標(biāo)系下的位移矢量有如下關(guān)系:

(2)

式中,T為位移轉(zhuǎn)換矩陣。同時(shí)對子結(jié)構(gòu)B的動力學(xué)方程的兩邊左乘以T,就完成了子結(jié)構(gòu)B的動力學(xué)方程從局部坐標(biāo)系向總體坐標(biāo)系的過渡。

為了方便模態(tài)綜合法的應(yīng)用,將各個(gè)子結(jié)構(gòu)在總體坐標(biāo)系下的位移分為內(nèi)部節(jié)點(diǎn)位移和界面節(jié)點(diǎn)位移,對質(zhì)量矩陣和剛度矩陣進(jìn)行分塊。以子結(jié)構(gòu)A為例,其動力學(xué)方程為:

(3)

式中,i代表內(nèi)部節(jié)點(diǎn)位移,m代表MPC耦合的節(jié)點(diǎn)位移,s代表扭轉(zhuǎn)彈簧耦合的節(jié)點(diǎn)位移。

對于子結(jié)構(gòu)A進(jìn)行模態(tài)分析,引入模態(tài)坐標(biāo)pA以及模態(tài)矩陣φA,則:

(4)

類似地,引入模態(tài)坐標(biāo)pB、pC以及模態(tài)矩陣φB、φC,對于子結(jié)構(gòu)B和子結(jié)構(gòu)C進(jìn)行模態(tài)分析。此時(shí),整體動力學(xué)方程可以寫成如下形式:

(5)

1.2參數(shù)化的氣動力模型

本文采用偶極子網(wǎng)格法獲得了折疊機(jī)翼的氣動力模型?？紤]到各個(gè)子結(jié)構(gòu)的升力面之間的相互干擾,每個(gè)升力面上的壓力分布Δp與氣動網(wǎng)格控制點(diǎn)的無量綱法洗幅值w滿足如下關(guān)系:

(6)

式中,D為非定常氣動力影響系數(shù)矩陣,ρ為密度,V為速度。

在簡諧振蕩的假設(shè)下,根據(jù)非定常氣動力的邊界條件,w與氣動網(wǎng)格節(jié)點(diǎn)的法向位移z有如下關(guān)系:

(7)

式中,k為減縮頻率,b為參考長度。由于氣動網(wǎng)格與結(jié)構(gòu)網(wǎng)格是獨(dú)立的,因此在進(jìn)行氣動彈性分析時(shí),需要進(jìn)行力和位移的轉(zhuǎn)換,此時(shí)可以采用樣條插值法來實(shí)現(xiàn),如無限板樣條,薄板樣條等。其位移轉(zhuǎn)換關(guān)系為:

(8)

式中,φ=blockdiag(φA,φB,φC),Gs1和Gs2為位移插值轉(zhuǎn)換矩陣。根據(jù)虛功原理,可以得到作用在結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)的等效力:

(9)

1.3氣動彈性系統(tǒng)的狀態(tài)方程

在整體動力學(xué)方程中引入作用在結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)的等效力:

(10)

式中,模態(tài)坐標(biāo)p中各個(gè)分量并不獨(dú)立,對各個(gè)子結(jié)構(gòu)進(jìn)行綜合,消去非獨(dú)立的模態(tài)坐標(biāo)便可以得到整體的氣動彈性模型。引入獨(dú)立模態(tài)坐標(biāo)變換矩陣Π,并左乘ΠT,則獨(dú)立模態(tài)坐標(biāo)下的氣動彈性方程為:

(11)

(12)

(13)

此時(shí)Q為離散形式的頻域氣動力。為了建立狀態(tài)空間的氣動彈性模型,利用有理函數(shù)近似,如Roger法,最小狀態(tài)法等可以獲得時(shí)域內(nèi)的氣動力。這里采用最小狀態(tài)法進(jìn)行頻域氣動力的有理近似,引入空氣狀態(tài)變量qa,則公式(13)表達(dá)為狀態(tài)空間的形式為:

(14)

V2ΠTA0Π。

此時(shí)A與折疊角θ,扭轉(zhuǎn)彈簧剛度KA、KB均有關(guān)。系統(tǒng)的穩(wěn)定性由A的特征值來確定:當(dāng)速度V小于顫振速度VF時(shí),全部特征值均位于復(fù)平面的左半部;當(dāng)速度V增加到一定值時(shí),出現(xiàn)一個(gè)系統(tǒng)特征值的實(shí)部由負(fù)值變?yōu)檎禃r(shí),系統(tǒng)會變?yōu)椴环€(wěn)定的。當(dāng)實(shí)部為零時(shí),此時(shí)的速度V即為顫振速度。盡管該方法可以預(yù)測系統(tǒng)的顫振速度,但是當(dāng)以簡化顫振邊界預(yù)測為目的時(shí),計(jì)算量較大。反之,采用Gram矩陣可以有效地避免特征值的計(jì)算,提高計(jì)算效率。

2基于Gram矩陣范數(shù)氣動彈性系統(tǒng)的穩(wěn)定性預(yù)測

在現(xiàn)代控制理論中,其反饋信息是由系統(tǒng)的狀態(tài)變量組合而成。但并不是所有的狀態(tài)變量都能測量到,于是提出了能否通過對輸出的測量獲得全部狀態(tài)變量的信息,便成了系統(tǒng)的可觀性問題。對于線性時(shí)不變的穩(wěn)定系統(tǒng),其可觀Gram矩陣可由下面的Lyapunov矩陣獲得[10]:

(15)

根據(jù)文獻(xiàn)[10],輸出的狀態(tài)變量所具有的能量,可以通過可觀Gram矩陣來體現(xiàn)。但是此時(shí),Wo矩陣不能采用公式(15)來求解。這是因?yàn)楣?15)只適用于穩(wěn)定系統(tǒng)。當(dāng)顫振發(fā)生時(shí),系統(tǒng)是不穩(wěn)定的。所以文獻(xiàn)[10]對其進(jìn)行了修正,使得Wo的求解不受系統(tǒng)穩(wěn)定性的影響。假設(shè)X是以下Riccati方程的解:

XF+FTX-XGGTX=0

(16)

(F+GMg)W+W(F+GMg)T+GGT=0

(17)

這里將Wo的矩陣范數(shù)σ作為參數(shù),此時(shí)不同速度V和C對應(yīng)一個(gè)σ,該值意味著機(jī)翼從來流中獲得的能量在結(jié)構(gòu)上的累積。當(dāng)C相同時(shí),隨著速度的增加,σ會達(dá)到最大值,該階模態(tài)不穩(wěn)定;倘若隨著速度的增加,σ沒有最大值出現(xiàn),則該階模態(tài)是穩(wěn)定的,進(jìn)而可以確定不穩(wěn)定模態(tài)以及顫振速度[10]。

3折疊機(jī)翼的顫振參數(shù)化分析

折疊機(jī)翼的各個(gè)子結(jié)構(gòu)的質(zhì)量矩陣和剛度矩陣通過DMAP語言獲得。采用DLM(doubletlatticemethod)獲得折疊機(jī)翼在頻域下的氣動力,通過最小狀態(tài)近似,最終獲得了折翼機(jī)翼的氣動彈性方程在狀態(tài)空間的表達(dá)形式。在公式(14)中,ρ=1.226kg/m3,M∞=0.0,采用前6階模態(tài)參與計(jì)算,KA=KB=1Nm/rad。最后利用文中的Gram矩陣的F范數(shù),確定了折疊機(jī)翼的顫振邊界,整個(gè)程序的編制均在MATLAB中完成。

3.1基于Gram矩陣范數(shù)顫振邊界預(yù)測的有效性驗(yàn)證

首先以機(jī)翼完全展開和完全折疊為例,即θ=0°和θ=120°,獲得了其顫振邊界。為了驗(yàn)證方法的有效性,將結(jié)果與采用特征值法的結(jié)果進(jìn)行了對比,在驗(yàn)證的過程中均沒有考慮結(jié)構(gòu)阻尼。

3.1.1機(jī)翼完全展開(θ=0°)

當(dāng)θ=0°時(shí),各階模態(tài)的范數(shù)σ隨速度的變化如圖2所示。

圖2　范數(shù)σ隨來流速度的變化(θ=0°)

從圖2中可以發(fā)現(xiàn),隨著速度的增加,第1階模態(tài)與第2階模態(tài)的σ,在接近42.7 m/s時(shí)達(dá)到了最大值,并且第2階模態(tài)σ的最大值明顯高于第1階,說明第2階模態(tài)不穩(wěn)定,其余各階模態(tài)的σ均沒有最大值出現(xiàn)。這一現(xiàn)象表明第1階模態(tài)與第2階模態(tài)相互耦合,使得結(jié)構(gòu)耗散的能量小于從來流中吸收的能量,達(dá)到顫振臨界點(diǎn)而發(fā)生顫振現(xiàn)象。與特征值法得到的結(jié)果對比如表1所示。

表1　顫振結(jié)果對比(θ=0°)

3.1.2機(jī)翼完全折疊(θ=120°)

當(dāng)θ=120°時(shí),各階模態(tài)的范數(shù)σ隨速度的變化如圖3所示。

圖3　范數(shù)σ隨來流速度的變化(θ=120°)

從圖3中可以看出,第3階模態(tài)和第4階模態(tài)的σ,在32.7 m/s附近,達(dá)到了最大值。隨著速度的增加,其余各階模態(tài)的σ沒有最大值出現(xiàn)。第3階模態(tài)的σ最大值明顯高于第4階,說明第3階模態(tài)不穩(wěn)定。與特征值法得到的結(jié)果對比如表2所示。

表2　顫振結(jié)果對比(θ=120°)

由于顫振發(fā)生時(shí),總是伴隨著模態(tài)之間的耦合,在某2個(gè)模態(tài)耦合較強(qiáng)的情況下,產(chǎn)生了能量轉(zhuǎn)換而發(fā)生顫振,所以在這2個(gè)例子中,均有2階模態(tài)的σ出現(xiàn)了最大值,意味著這2階模態(tài)的耦合作用較強(qiáng),使得結(jié)構(gòu)從氣流中吸收能量,達(dá)到顫振臨界點(diǎn)而發(fā)生顫振。這些與顫振發(fā)生的機(jī)理保持一致。從表1和表2 的對比結(jié)果可以看出:本文方法與特征值法獲得的顫振速度的相對誤差較小,對于不穩(wěn)定模態(tài)的預(yù)測也吻合。表明本文方法可以獲得與特征值法精度相當(dāng)?shù)慕Y(jié)果。

3.2顫振邊界的分析

顫振速度隨折疊角的變化情況如圖4 所示。

圖4　顫振速度隨折疊角的變化

隨著折疊角的增加,顫振速度并不是折疊角的單一函數(shù)。在某些折疊角度范圍內(nèi),如45°≤θ≤90°,顫振速度隨著折疊角的增加而增加;在某些折疊角度范圍內(nèi),如30°<θ<45°,顫振速度隨著折疊角的增加而減小。折疊角在接近30°時(shí),顫振速度達(dá)到了最大值,即48.36 m/s;折疊角在接近105°時(shí),顫振速度達(dá)到了最小值,即29.59 m/s。即在整個(gè)折疊過程中,顫振速度的變化幅度很大。顯然,折疊角對顫振速度有很大的影響。這主要是由于折疊角的變化,使得結(jié)構(gòu)的質(zhì)量矩陣,剛度矩陣以及氣動模型均發(fā)生了很大的變化。從表3 中,可以看出,隨著折疊角的增加,顫振不穩(wěn)定的模態(tài)也發(fā)生了變化,從第2階模態(tài)過渡到了第3階模態(tài)。

表3　不穩(wěn)定模態(tài)隨折疊角的變化

值得注意的是,在不同的折疊角度下,隨著模態(tài)阻尼比的增加,顫振速度均有所提高。這是由于較高的模態(tài)阻尼比能夠抑制不穩(wěn)定的模態(tài)之間的耦合,提高結(jié)構(gòu)的能量耗散,減小從來流中吸收的能量在結(jié)構(gòu)上的累積,推遲顫振現(xiàn)象的發(fā)生。

4結(jié)論

本文建立了折疊機(jī)翼顫振分析的參數(shù)化氣動彈性模型,該模型以折疊角和內(nèi)外鉸鏈剛度為參數(shù),可以有效地形成不同構(gòu)型下的折疊機(jī)翼的氣動彈性模型。根據(jù)控制理論中的Gram矩陣范數(shù)獲得了折疊機(jī)翼的顫振邊界,可以得出以下結(jié)論:

1) 提出一種適合折疊機(jī)翼的參數(shù)化的氣動彈性模型。該模型具有以下優(yōu)點(diǎn):適用于類似折疊機(jī)翼的多段機(jī)翼結(jié)構(gòu);結(jié)構(gòu)模型和氣動力模型均不需要重復(fù)建立。

2) 采用Gram矩陣范數(shù)預(yù)測氣動彈性系統(tǒng)的穩(wěn)定性,可以定性地分析參與顫振耦合的模態(tài)信息,避免高維氣動彈性系統(tǒng)的特征值求解,提高計(jì)算效率。

3) Gram矩陣的范數(shù),不僅僅局限于F范數(shù),如1范數(shù),無窮范數(shù)等。當(dāng)采用不同范數(shù)時(shí),得到的顫振速度基本不變。

4) 由于折疊角對機(jī)翼結(jié)構(gòu)以及氣動特性有很明顯的影響,使得顫振速度對折疊角很敏感;并且隨著折疊角的增加,出現(xiàn)了不穩(wěn)定模態(tài)之間的過渡。

5) 在一定的范圍內(nèi),結(jié)構(gòu)的模態(tài)阻尼比可以提高折疊機(jī)翼的顫振速度。

參考文獻(xiàn)：

[1]Wilson J R. Morphing UAVs Change the Shape of Warfare[J]. Aerospace America, 2004, 42(2): 28-29

[2]Tang D, Dowell E H. Theoretical and Experimental Aeroelastic Study for Folding Wing Strucutres[J]. Journal of Aircraft, 2008, 45(4):1136-1147

[3]Snyder M P, Sanders B, Eastep F E, Frank G J. Vibration and Flutter Characteristics of a Folding Wing [J]. Journal of Aircraft, 2009, 46(3): 791-799

[4]Wamg I, Gibbs S C, Dowell E H. Aeroelastic Model of Mutisegmented Folding Wing: Theory and Experiment[J]. Journal of Aircraft, 2012, 49(3): 911-921

[5]Zhao Youhui, Hu Haiyan. Paremeterized Aeroelastic Modeling and Flutter Analysis for a Folding Wing[J]. Journal of Sound and Vibration, 2012,331(2): 308-324

[6]楊智春, 谷迎松, 夏巍. 基于矩陣奇異值理論的顫振分析新方法[J]. 航空學(xué)報(bào), 2009, 30(6): 985-989

Yang Zhichun, Gu Yingsong, Xia Wei. New Type of Flutter Solution Based on Matrix Singularity[J]. Acta Aeronautica et Astronautica Sinica, 2009,30(6): 985-989 (in Chinese)

[7]Moore B C. Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduciton[J]. IEEE Trans on Automatic Control, 1981, 26(1): 17-32

[8]Hyounsurk Roh,Yongjin Park. Actuator and Exciter Placement for Flexible Structures[J]. Journal of Guidance, Control and Dynamics, 1997, 20(5): 850-856

[9]Rowley C W. Model Reduction for Fluids, Using Balanced Proper Orthogonal Decomposition[J]. International Journal of Bifurcation and Chaos, 2005, 15(3): 997-1013

[10] Bueno D D, Marqui C R, Goes L C S, Goncalves P J P. The Use of Gramian Matrices for Aeroelastic Stability Analysis[C]//Mathematical Problems in Engineering, 2013

Parametric Aeroelastic Modeling and Flutter Analysis for a Folding Wing

Ni Yingge, Hou Chi, Wan Xiaopeng, Zhao Meiying

(College of Aeronautics, Northwestern Polytechnical University, Xi′an 710072, China)

Abstract:As for the characteristics of a folding wing, a parametric aeroelastic model for flutter analysis is proposed in this paper. Parametric structural model is achieved based on modal synthesis method. Parametric aerodynamic model is established using doublet lattice method. And the flutter boundary prediction of an aeroelastic system is elaborated using Gram matrix norm. The fully extended and completely folded configurations are taken as examples. The results of this method are compared with those of eigenvalue method to verify the correctness. Some analyses of different folding angles are performed. Flutter boundary is very sensitive to the folding angle. With increasing folding angle, the flutter mode changes. Higher modal damping ratio can raise flutter speed and delay flutter phenomenon.

Key words:Folding wing, Gram matrix, parametric model, aeroelasticity; flutter, eigenvalues and eigenfunctions, finite element method, MATLAB, Riccati equation, stability

中圖分類號：V271.4

文獻(xiàn)標(biāo)志碼：A

文章編號:1000-2758(2015)05-0788-06

作者簡介：倪迎鴿(1987—),女,西北工業(yè)大學(xué)博士研究生,主要從事變體飛機(jī)的一體化設(shè)計(jì)的研究。

收稿日期：2015-04-09