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

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

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

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

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

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

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

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

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

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

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

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

圖1　折疊機翼的有限元模型和坐標系　(每個圓點代表鉸鏈所處位置)

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

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

(1)

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

(2)

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

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

(3)

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

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

(4)

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

(5)

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

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

(6)

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

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

(7)

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

(8)

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

(9)

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

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

(10)

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

(11)

(12)

(13)

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

(14)

V2ΠTA0Π。

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

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

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

(15)

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

XF+FTX-XGGTX=0

(16)

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

(17)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.2顫振邊界的分析

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

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

隨著折疊角的增加,顫振速度并不是折疊角的單一函數(shù)。在某些折疊角度范圍內(nèi),如45°≤θ≤90°,顫振速度隨著折疊角的增加而增加;在某些折疊角度范圍內(nèi),如30°<θ<45°,顫振速度隨著折疊角的增加而減小。折疊角在接近30°時,顫振速度達到了最大值,即48.36 m/s;折疊角在接近105°時,顫振速度達到了最小值,即29.59 m/s。即在整個折疊過程中,顫振速度的變化幅度很大。顯然,折疊角對顫振速度有很大的影響。這主要是由于折疊角的變化,使得結(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é)論

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

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

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

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

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

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

參考文獻：

[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]. 航空學報, 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

文獻標志碼：A

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

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

收稿日期：2015-04-09