大展弦比柔性機翼非線性顫振研究

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

(1.陸軍工程大學石家莊校區 無人機工程系， 石家莊 050003； 2.石家莊鐵道大學 工程力學系， 石家莊 050003；3.中國人民解放軍71939部隊， 濟南 250000)

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

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

本研究對大展弦比柔性機翼進行非線性顫振分析，求解思路如下：先將柔性機翼視為曲梁，結合非線性ONERA氣動力模型得到機翼顫振微分方程，為傳遞函數法求解做好準備。然后，將顫振微分方程轉換為狀態空間方程形式，通過求解復特征值問題，獲得了大展弦比柔性機翼在非線性氣動力作用下的顫振速度和顫振頻率。

1 柔性機翼非線性顫振分析模型

1.1 曲梁單元振動微分方程

將柔性機翼視為曲梁，根據文獻[9]，考慮曲梁質心軸與彈性軸不重合的情形，可獲得曲梁六自由度振動方程為：

(1)

如圖1所示，曲梁曲率為常數R，u,w,v分別為曲梁單元沿坐標軸ξ,η,ζ三個方向的位移，ψζ,ψη,φξ分別為曲梁單元繞ζ,η,ξ軸的扭轉角。ρ為曲梁的密度，A為曲梁單元截面積。E,G為拉伸彈性模量和剪切彈性模量，k為剪切剛度修正系數。ζa為彈性軸到質心軸的距離。Iζ,Iη,J分別為繞坐標軸ζ,η,ξ三個方向的慣性矩。t為時間。qξ(ξ,t)、qη(ξ,t)、qζ(ξ,t)為沿坐標軸ξ,η,ζ三個方向的分布力，mξ(ξ,t)、mη(ξ,t)、mζ(ξ,t)為繞坐標軸ξ,η,ζ三個方向的分布彎矩。

1.2 ONEAR非線性氣動力模型

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

(2)

式中：L表示機翼升力；V是來流速度；α為瞬時攻角；φξ為有效攻角，相當于機翼實際扭轉角；b為機翼半弦長；Cz表示空氣動力系數；Cza表示線性氣動力部分對應的氣動力函數；Czb表示非線性氣動力部分對應的氣動力系數。參數下標z=L時，表示與升力有關的系數；當下標z=M時，表示與力矩有關的系數。由于ONEAR線性部分與經典Theodorson理論一致，所以其線性部分參數選取如表1所示。

表1 ONERA模型線性部分參數值

非線性氣動力部分的有關參數有r1z，r2z，r3z。

為簡化計算，通常用折線來近似代替靜態空氣動力曲線，如圖2所示。

ΔCz可以表示為：

(3)

(4)

這里w1/4表示四分之一弦長處的位移，根據幾何關系，w1/4可以表示為：

(5)

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

(6)

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

z=LorM

(7)

單位展長上的升力和相應的俯仰力矩可以表示為：

(8)

1.3 機翼單元顫振微分方程

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

(9)

2 傳遞函數法求解

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

(10)

其中，氣動力系數CL和CM的Fourier變換可以表示為：

(11)

根據傳遞函數方法，定義狀態向量：

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

(12)

其中，ge(ξ,ω)=0，轉移矩陣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表達式為：

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

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

(15)

其中，

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

(16)

由于大變形后的機翼劃分為若干常曲率曲梁單元來描述整個機翼，其求解方程可借鑒有限元方法的思想進行組集。曲梁單元截面上的內力為：

(17)

其中：Mη、Mζ分別為曲梁單元繞η、ζ軸彎矩；Tξ分別為曲梁單元繞ξ軸的扭矩；Nξ為沿ξ軸的軸力；Qη和Qζ分別為沿η、ζ軸的剪力。

將曲梁單元內力寫成矩陣形式：

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

(18)

(19)

將式(15)代入式(18)，可的達到曲梁單元端點處的內力：

(20)

按照有限元組集方法對各節點統一編號，并進行組集拼接，可得到機翼整體平衡方程為：

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

(21)

其中，

(22)

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

f=0

(23)

當機翼顫振時，γ(ω)應有非零解，此時須滿足條件

det[K(ω,V)]=0

(24)

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

(25)

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

3 算例分析與討論

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

表2 大柔性機翼模型參數

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

使用表2中的參數進行非線性氣動力條件下的顫振計算，通過比較線性氣動力和非線性氣動力求解的顫振計算結果，可以發現，運用非線性氣動力求解的顫振速度要低于線性氣動力求解的顫振速度，是因為隨著彈性變形的增加，機翼的氣動力系數減小，其相應的氣動特性也降低，隨之顫振速度也降低。顫振頻率通過對比發現，隨著變形程度的加大，變化趨勢相同，但是非線性氣動力計算得到的顫振頻率變化更加劇烈。

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

密度和展弦比是大展弦比機翼的重要設計參數(見圖7、圖8)，以表2中的參數為基準，以變形條件6.25%為例，機翼顫振速度隨機翼半展現比的變化比較簡單，非線性和線性的結果都是單調遞減的，也就是說展現比越大，柔性機翼的氣動彈性穩定性越差。顫振速度隨機翼密度的增加而減小，但是線性和非線性的計算誤差隨著密度增加而增大，這是由于隨著機翼的質量增加，所受氣動力增加，導致兩種氣動力的計算結果誤差變大。

最后來研究彈性軸在機翼中的相對位置對大柔性無人機氣動彈性穩定的影響，大柔性無人機的彈性軸在弦向位置是非常重要的氣動彈性設計量，如圖9所示，橫軸表示彈性軸在弦長相對位置的百分比，縱軸為顫振速度，隨著彈性軸后移，顫振速度降低較快，即彈性軸后移嚴重削弱了大柔性無人機的氣動彈性穩定性，非線性顫振速度相比于線性計算結果降低的比較慢。

4 結論

1) 使用非線性氣動力計算得到機翼顫振速度相較于使用線性氣動力計算得到結果偏低，表明將氣動非線性引入計算，機翼的氣動特性降低，同時影響到機翼的顫振特性。

2) 不同參數變化對比線性氣動力和非線性氣動力計算得到的顫振結果相對誤差也不同，當相對誤差較小時可以采用計算量較小的線性氣動力模型進行計算，當誤差較大時，則需要采用非線性氣動力計算才能反映更真實的機翼顫振特性。