張忠源,段靜波,路 平,葉 廷

(1.陸軍工程大學(xué)石家莊校區(qū) 無(wú)人機(jī)工程系， 石家莊 050003； 2.石家莊鐵道大學(xué) 工程力學(xué)系， 石家莊 050003；3.中國(guó)人民解放軍71939部隊(duì)， 濟(jì)南 250000)

當(dāng)機(jī)翼攻角大于失速攻角時(shí)，機(jī)翼升力會(huì)迅速減小，對(duì)飛行器的安全造成極大威脅[1]。所以實(shí)際情況中，機(jī)翼所受氣動(dòng)力是非線性的，尤其對(duì)于大柔性機(jī)翼，受氣動(dòng)力作用，機(jī)翼變形更加明顯[2]，極大提高了機(jī)翼產(chǎn)生失速的可能性。為解決這一問(wèn)題，就要在計(jì)算大柔性機(jī)翼顫振時(shí)引入非線性氣動(dòng)力模型，本文將ONERA非線性氣動(dòng)力模型引入大展弦比柔性機(jī)翼顫振問(wèn)題，同時(shí)基于半解析半數(shù)值的傳遞函數(shù)方法[3]，該方法求解過(guò)程簡(jiǎn)潔和統(tǒng)一，邊界條件處理規(guī)范和方便，可以更加快捷、有效、準(zhǔn)確地計(jì)算出機(jī)翼的顫振特性。

由于大柔性機(jī)翼扭轉(zhuǎn)會(huì)產(chǎn)生非線性氣動(dòng)力，CC Xie等[4]針對(duì)這個(gè)問(wèn)題進(jìn)行了研究，用平面雙點(diǎn)陣方法計(jì)算頻域內(nèi)的非定常氣動(dòng)力，忽略偏轉(zhuǎn)翼的彎曲效應(yīng)。在給定的載荷條件下，對(duì)系統(tǒng)進(jìn)行氣彈性穩(wěn)定性分析。同樣針對(duì)大展弦比大柔性機(jī)翼，密歇根大學(xué)C Cesnik[5]團(tuán)隊(duì)也進(jìn)行了深入研究，并搭建了收集幾何非線性氣動(dòng)彈性響應(yīng)的數(shù)據(jù)實(shí)驗(yàn)平臺(tái)，為飛機(jī)提供可在飛行中測(cè)量的特定氣動(dòng)彈性特征，例如，耦合的剛性、彈性體不穩(wěn)定性，陣風(fēng)期間的大的機(jī)翼偏轉(zhuǎn)等。D Tang[6]將柔性機(jī)翼的氣動(dòng)彈性分析與風(fēng)洞試驗(yàn)相結(jié)合，介紹了一種彈性載荷作用下柔性大展弦比翼型氣動(dòng)彈性模型的理論氣動(dòng)彈性模型。MJ Patil等[7,8]提出使用完整飛機(jī)模型的氣動(dòng)彈性特性來(lái)獲得結(jié)果，由于機(jī)翼具有大柔性，飛機(jī)整體的飛行動(dòng)態(tài)特性也會(huì)發(fā)生變化，并用嚴(yán)格的非線性氣動(dòng)彈性來(lái)分析。

本研究對(duì)大展弦比柔性機(jī)翼進(jìn)行非線性顫振分析，求解思路如下：先將柔性機(jī)翼視為曲梁，結(jié)合非線性O(shè)NERA氣動(dòng)力模型得到機(jī)翼顫振微分方程，為傳遞函數(shù)法求解做好準(zhǔn)備。然后，將顫振微分方程轉(zhuǎn)換為狀態(tài)空間方程形式，通過(guò)求解復(fù)特征值問(wèn)題，獲得了大展弦比柔性機(jī)翼在非線性氣動(dòng)力作用下的顫振速度和顫振頻率。

1 柔性機(jī)翼非線性顫振分析模型

1.1 曲梁?jiǎn)卧駝?dòng)微分方程

將柔性機(jī)翼視為曲梁，根據(jù)文獻(xiàn)[9]，考慮曲梁質(zhì)心軸與彈性軸不重合的情形，可獲得曲梁六自由度振動(dòng)方程為：

(1)

如圖1所示，曲梁曲率為常數(shù)R，u,w,v分別為曲梁?jiǎn)卧刈鴺?biāo)軸ξ,η,ζ三個(gè)方向的位移，ψζ,ψη,φξ分別為曲梁?jiǎn)卧@ζ,η,ξ軸的扭轉(zhuǎn)角。ρ為曲梁的密度，A為曲梁?jiǎn)卧孛娣e。E,G為拉伸彈性模量和剪切彈性模量，k為剪切剛度修正系數(shù)。ζa為彈性軸到質(zhì)心軸的距離。Iζ,Iη,J分別為繞坐標(biāo)軸ζ,η,ξ三個(gè)方向的慣性矩。t為時(shí)間。qξ(ξ,t)、qη(ξ,t)、qζ(ξ,t)為沿坐標(biāo)軸ξ,η,ζ三個(gè)方向的分布力，mξ(ξ,t)、mη(ξ,t)、mζ(ξ,t)為繞坐標(biāo)軸ξ,η,ζ三個(gè)方向的分布彎矩。

1.2 ONEAR非線性氣動(dòng)力模型

在忽略機(jī)翼重力影響條件下，機(jī)翼顫振時(shí)的外力為氣動(dòng)力產(chǎn)生的分布升力以及分布扭矩。本文采用片條理論進(jìn)行非定常氣動(dòng)力計(jì)算。對(duì)于大展弦比大柔性機(jī)翼，在正常飛行期間局部迎角可能非常大，這將導(dǎo)致機(jī)翼失速，即強(qiáng)烈的非線性不穩(wěn)定失速現(xiàn)象。在本章中，使用了ONEAR[10]動(dòng)態(tài)失速模型，并且給出了作用在四分之一弦長(zhǎng)處的升力和力矩的公式：

(2)

式中：L表示機(jī)翼升力；V是來(lái)流速度；α為瞬時(shí)攻角；φξ為有效攻角，相當(dāng)于機(jī)翼實(shí)際扭轉(zhuǎn)角；b為機(jī)翼半弦長(zhǎng)；Cz表示空氣動(dòng)力系數(shù)；Cza表示線性氣動(dòng)力部分對(duì)應(yīng)的氣動(dòng)力函數(shù)；Czb表示非線性氣動(dòng)力部分對(duì)應(yīng)的氣動(dòng)力系數(shù)。參數(shù)下標(biāo)z=L時(shí)，表示與升力有關(guān)的系數(shù)；當(dāng)下標(biāo)z=M時(shí)，表示與力矩有關(guān)的系數(shù)。由于ONEAR線性部分與經(jīng)典Theodorson理論一致，所以其線性部分參數(shù)選取如表1所示。

表1 ONERA模型線性部分參數(shù)值

非線性氣動(dòng)力部分的有關(guān)參數(shù)有r1z，r2z，r3z。

為簡(jiǎn)化計(jì)算，通常用折線來(lái)近似代替靜態(tài)空氣動(dòng)力曲線，如圖2所示。

ΔCz可以表示為：

(3)

(4)

這里w1/4表示四分之一弦長(zhǎng)處的位移，根據(jù)幾何關(guān)系，w1/4可以表示為：

(5)

將式(5)代入式(4)并整理得：

(6)

將式(6)代入式(2)，可得：

z=LorM

(7)

單位展長(zhǎng)上的升力和相應(yīng)的俯仰力矩可以表示為：

(8)

1.3 機(jī)翼單元顫振微分方程

將式(8)代入式(1)，即可得到大變形柔性機(jī)翼的顫振微分方程。

(9)

2 傳遞函數(shù)法求解

對(duì)機(jī)翼的顫振微分方程(9)進(jìn)行Fourier變換，并整理得：

(10)

其中，氣動(dòng)力系數(shù)CL和CM的Fourier變換可以表示為：

(11)

根據(jù)傳遞函數(shù)方法，定義狀態(tài)向量：

將式(10)改寫為狀態(tài)空間形式方程：

(12)

其中，ge(ξ,ω)=0，轉(zhuǎn)移矩陣Fe(ω,V)為12×12的方陣，其非零元素為：

Fe(8,10)=-1,Fe(9,10)=1

(13)

邊界條件為：

Mbeηe(ξ=0,ω)+Nbeηe(ξ=1,ω)=γe(ω)

(14)

其中，γe(ω)為由位移和力組成的列向量，Mbe、Nbe分別為單元邊界條件選擇矩陣，γe(ω)、Mbe、Nbe表達(dá)式為：

式(12)的傳遞函數(shù)解為：

ηe(ξ,ω)=He(ξ,ω,V)γe(ω)

(15)

其中，

He(ξ,ω,V)=eFe(ω,V)ξ[Mb+NbeFe(ω,V)]-1

(16)

由于大變形后的機(jī)翼劃分為若干常曲率曲梁?jiǎn)卧獊?lái)描述整個(gè)機(jī)翼，其求解方程可借鑒有限元方法的思想進(jìn)行組集。曲梁?jiǎn)卧孛嫔系膬?nèi)力為：

(17)

其中：Mη、Mζ分別為曲梁?jiǎn)卧@η、ζ軸彎矩；Tξ分別為曲梁?jiǎn)卧@ξ軸的扭矩；Nξ為沿ξ軸的軸力；Qη和Qζ分別為沿η、ζ軸的剪力。

將曲梁?jiǎn)卧獌?nèi)力寫成矩陣形式：

Qe(ξ)=Qeη(ξ)ηe(ξ,ω)

(18)

(19)

將式(15)代入式(18)，可的達(dá)到曲梁?jiǎn)卧它c(diǎn)處的內(nèi)力：

(20)

按照有限元組集方法對(duì)各節(jié)點(diǎn)統(tǒng)一編號(hào)，并進(jìn)行組集拼接，可得到機(jī)翼整體平衡方程為：

K(ω,V)γ(ω)=f

(21)

其中，

(22)

其中：K(ω,V)可視為整體剛度矩陣；γ(ω)可視為整體節(jié)點(diǎn)位移向量；f為各單元節(jié)點(diǎn)內(nèi)力拼裝成的向量。由于本文中將機(jī)翼的氣動(dòng)力與機(jī)翼作為一個(gè)完整的系統(tǒng)來(lái)考慮，除此之外，機(jī)翼沒(méi)有受到其他外力作用，因而根據(jù)單元節(jié)點(diǎn)內(nèi)力與外載荷平衡，可得出

f=0

(23)

當(dāng)機(jī)翼顫振時(shí)，γ(ω)應(yīng)有非零解，此時(shí)須滿足條件

det[K(ω,V)]=0

(24)

由于K(ω,V)為復(fù)矩陣，其行列式值等于零的必要條件為矩陣行列式值的實(shí)部與虛部均為零，即

(25)

上式兩個(gè)方程包含兩個(gè)未知變量V和ω，因而求解式(25)可得到V和ω的解。其中，V即為機(jī)翼顫振速度，ω即為機(jī)翼顫振頻率。通常滿足式(25)條件的解可能有多個(gè)，其中V最小的一組解即為機(jī)翼的顫振特性。

3 算例分析與討論

前兩節(jié)推導(dǎo)出了非線性氣動(dòng)力影響下的機(jī)翼顫振求解方法，為研究非線性氣動(dòng)力與線性氣動(dòng)力對(duì)機(jī)翼顫振造成的差異，選取某大柔性機(jī)翼模型進(jìn)行非線性氣動(dòng)力和線性氣動(dòng)力的顫振計(jì)算，機(jī)翼模型參數(shù)如表2所示，通過(guò)改變機(jī)翼線密度，展弦比，抗彎剛度和抗扭剛度等參數(shù)，對(duì)比線性與非線性氣動(dòng)力模型下的顫振結(jié)果的差異。

表2 大柔性機(jī)翼模型參數(shù)

這里取3種典型變形進(jìn)行分析，記變形一、變形二和變形三，對(duì)應(yīng)的翼尖位移取半展長(zhǎng)的3.125%、6.250%、12.50%，分別對(duì)3種變形使用線性和非線性顫振分析。計(jì)算得到3種變形下的線性和非線性的機(jī)翼顫振速度和顫振頻率，如圖3和圖4所示。

使用表2中的參數(shù)進(jìn)行非線性氣動(dòng)力條件下的顫振計(jì)算，通過(guò)比較線性氣動(dòng)力和非線性氣動(dòng)力求解的顫振計(jì)算結(jié)果，可以發(fā)現(xiàn)，運(yùn)用非線性氣動(dòng)力求解的顫振速度要低于線性氣動(dòng)力求解的顫振速度，是因?yàn)殡S著彈性變形的增加，機(jī)翼的氣動(dòng)力系數(shù)減小，其相應(yīng)的氣動(dòng)特性也降低，隨之顫振速度也降低。顫振頻率通過(guò)對(duì)比發(fā)現(xiàn)，隨著變形程度的加大，變化趨勢(shì)相同，但是非線性氣動(dòng)力計(jì)算得到的顫振頻率變化更加劇烈。

將氣動(dòng)彈性系統(tǒng)中的彎曲剛度、扭轉(zhuǎn)剛度、彈性軸位置、半展現(xiàn)比和機(jī)翼線密度等參數(shù)作為變量，以表2中的參數(shù)為基準(zhǔn)，以變形條件6.25%為例，對(duì)比線性氣動(dòng)力和非線性氣動(dòng)力求解的顫振計(jì)算結(jié)果的變化規(guī)律。定義垂直彎曲剛度與表2中垂直彎曲剛度值的比為彎曲剛度比：ηb=EIb/EIζ，扭轉(zhuǎn)剛度與表2中扭轉(zhuǎn)剛度值的比：ηt=GJt/GJ，改變垂直彎曲剛度和扭轉(zhuǎn)剛度，得到的線性與非線性顫振速度變化趨勢(shì)如圖5、圖6所示。隨著垂直彎曲剛度和扭轉(zhuǎn)剛度的增加，顫振速度也隨之提高，增加垂直彎曲剛度和扭轉(zhuǎn)剛度可以有效改善其氣動(dòng)彈性穩(wěn)定性，可見在進(jìn)行大柔性無(wú)人機(jī)設(shè)計(jì)時(shí)需要著重注意彎曲和扭轉(zhuǎn)剛度的設(shè)計(jì)。

密度和展弦比是大展弦比機(jī)翼的重要設(shè)計(jì)參數(shù)(見圖7、圖8)，以表2中的參數(shù)為基準(zhǔn)，以變形條件6.25%為例，機(jī)翼顫振速度隨機(jī)翼半展現(xiàn)比的變化比較簡(jiǎn)單，非線性和線性的結(jié)果都是單調(diào)遞減的，也就是說(shuō)展現(xiàn)比越大，柔性機(jī)翼的氣動(dòng)彈性穩(wěn)定性越差。顫振速度隨機(jī)翼密度的增加而減小，但是線性和非線性的計(jì)算誤差隨著密度增加而增大，這是由于隨著機(jī)翼的質(zhì)量增加，所受氣動(dòng)力增加，導(dǎo)致兩種氣動(dòng)力的計(jì)算結(jié)果誤差變大。

最后來(lái)研究彈性軸在機(jī)翼中的相對(duì)位置對(duì)大柔性無(wú)人機(jī)氣動(dòng)彈性穩(wěn)定的影響，大柔性無(wú)人機(jī)的彈性軸在弦向位置是非常重要的氣動(dòng)彈性設(shè)計(jì)量，如圖9所示，橫軸表示彈性軸在弦長(zhǎng)相對(duì)位置的百分比，縱軸為顫振速度，隨著彈性軸后移，顫振速度降低較快，即彈性軸后移嚴(yán)重削弱了大柔性無(wú)人機(jī)的氣動(dòng)彈性穩(wěn)定性，非線性顫振速度相比于線性計(jì)算結(jié)果降低的比較慢。

4 結(jié)論

1) 使用非線性氣動(dòng)力計(jì)算得到機(jī)翼顫振速度相較于使用線性氣動(dòng)力計(jì)算得到結(jié)果偏低，表明將氣動(dòng)非線性引入計(jì)算，機(jī)翼的氣動(dòng)特性降低，同時(shí)影響到機(jī)翼的顫振特性。

2) 不同參數(shù)變化對(duì)比線性氣動(dòng)力和非線性氣動(dòng)力計(jì)算得到的顫振結(jié)果相對(duì)誤差也不同，當(dāng)相對(duì)誤差較小時(shí)可以采用計(jì)算量較小的線性氣動(dòng)力模型進(jìn)行計(jì)算，當(dāng)誤差較大時(shí)，則需要采用非線性氣動(dòng)力計(jì)算才能反映更真實(shí)的機(jī)翼顫振特性。