一種基于CR理論的大柔性機(jī)翼非線性氣動彈性求解方法

王偉1，周洲1，祝小平2，王睿1

(1.西北工業(yè)大學(xué)航空學(xué)院，西安710072；2. 西北工業(yè)大學(xué)無人機(jī)研究所，西安710065)

摘要：大展弦比大柔性機(jī)翼在氣動載荷的作用下，產(chǎn)生較大的彈性變形，其慣性特性、剛度特性、動氣動彈性特性等亦發(fā)生較大改變，常規(guī)的線性氣動彈性分析方法不再適用?；贑o-rotational(CR)理論，推導(dǎo)了機(jī)翼變形后的切線剛度矩陣和質(zhì)量矩陣，建立了考慮幾何非線性效應(yīng)的大柔性機(jī)翼結(jié)構(gòu)動力學(xué)模型；耦合改進(jìn)的ONERA非線性非定常氣動力模型，提出了一種適用于大柔性機(jī)翼的非線性氣動彈性求解方法。采用Newmark直接數(shù)值積分法及松耦合技術(shù)在時域內(nèi)對氣動彈性運(yùn)動方程進(jìn)行求解，對所提出的非線性氣動彈性求解方法的正確性和精度進(jìn)行了驗證，并研究了大柔性機(jī)翼的極限環(huán)顫振特性。研究表明：適用于大柔性機(jī)翼完整的非線性氣動彈性建模需要考慮機(jī)翼結(jié)構(gòu)大變形和非定常氣動力動態(tài)失速等非線性因素；彎曲變形可降低臨界極限環(huán)顫振速度的15%以上，而前移彈性軸能夠有效的提高臨界極限環(huán)顫振速度；所提出的非線性氣動彈性求解方法具有較好的精度和效率，滿足大柔性機(jī)翼非線性氣動彈性的求解需求。

關(guān)鍵詞：非線性氣動彈性；極限環(huán)顫振；CR理論；非定常氣動力；動態(tài)失速；Newmark積分法

中圖分類號:V211.47

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

DOI:10.13465/j.cnki.jvs.2015.19.010

Abstract:Very flexible wings under aerodynamic loads tend to produce larger deformation, it results in significant changes in inertial and stiffness characteristics, and dynamic aeroelastic features, the linear aeroelastic analysis method is no longer applicable. Here, based on the co-rotational(CR) theory, the tangent stiffness matrix and mass matrix of a wing after deformation were derived, the structural dynamic model of very flexible wings considering geometric nonlinearity was then established. Coupled with ONERA dynamic stall model, an efficient method to solve nonlinear aeroelasticity of very flexible wings was proposed. Using Newmark direct integration method and loose coupled algorithms, a numerical procedure for solving nonlinear aeroelastic dynamic equations was presented, and the efficiency and precision of the method were verified through tests. The results showed that structural and aerodynamic nonlinearities should be considered for complete nonlinear dynamic aeroelastic simulations of very flexible wings; the wing’s critical limit cycle oscillation speed decreases 15% or more due to its bending deformation, but it increases through shifting forward the wing’s elastic axis; the proposed method has a good precision and efficiency, and meets requirements of nonlinear aeroelastic analysis of very flexible wings.

基金項目：國家自然科學(xué)基金(51475387)

收稿日期：2014-06-12修改稿收到日期：2014-09-30

A CR theory-based approach for solving nonlinear aeroelasticity of very flexible wings

WANGWei1,ZHOUZhou1,ZHUXiao-ping2,WANGRui1(1. College of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China；2. UAV Research Institute, Northwestern Polytechnical University, Xi’an 710065, China)

Key words:nonlinear aeroelasticity; limit cycle oscillation; CR theory; unsteady aerodynamics loads; dynamic stall; Newmark integration method

高空超長航時太陽能無人機(jī)在飛行高度上能夠填補(bǔ)低軌道衛(wèi)星和常規(guī)能源飛機(jī)的空白，執(zhí)行環(huán)境監(jiān)測、通信中繼等任務(wù)，具有覆蓋區(qū)域廣、持續(xù)時間長和花費(fèi)低等特點，由于高空長航時的性能要求，這類飛機(jī)一般具有較小的結(jié)構(gòu)重量和超大展弦比機(jī)翼，在氣動載荷的作用下，機(jī)翼彎曲變形可達(dá)半翼展的50%，機(jī)翼的結(jié)構(gòu)動力學(xué)特性和氣動彈性特性隨變形的增加而產(chǎn)生較大改變。建立高精度的大柔性機(jī)翼結(jié)構(gòu)動力學(xué)模型進(jìn)而盡可能高精度、高效率的預(yù)測大柔性大變形機(jī)翼的非線性氣動彈性特征是這類飛機(jī)的主要研究內(nèi)容之一。由于面內(nèi)、面外彎曲變形和扭轉(zhuǎn)變形與結(jié)構(gòu)當(dāng)前的幾何位形有關(guān)，并且機(jī)翼局部剖面的有效攻角可能很大，這類機(jī)翼的氣動彈性建模需要考慮結(jié)構(gòu)幾何非線性和大攻角時非定常氣動力動態(tài)失速的影響。

目前考慮幾何非線性效應(yīng)的非線性氣動彈性運(yùn)動方程的求解方法主要有直接數(shù)值積分法和線性化后的頻域法[1-6]。直接數(shù)值積分法可以考慮結(jié)構(gòu)大振幅運(yùn)動以及氣動力的動態(tài)失速效應(yīng)等非線性因素；線性化的頻域法一般假設(shè)結(jié)構(gòu)在靜平衡位置附近作小振幅振動，然后耦合線性頻域非定常氣動力，采用經(jīng)典的P-K法或V-g法求解而無法考慮氣動非線性的影響[5]。Patil等[4]基于經(jīng)典ONERA動態(tài)失速模型建立了可以描述二維翼段準(zhǔn)靜態(tài)失速的非線性非定常氣動力模型，并耦合混合變分形式的彎-彎-扭耦合本征梁模型建立了非線性氣動彈性運(yùn)動方程，開展了針對大柔性機(jī)翼的非線性氣動彈性的研究。Shams等[6]采用Duhamel積分形式的Wagner函數(shù)描述線性時域非定常氣動力并耦合非線性Euler-Bernouli理論梁模型建立了描述大柔性機(jī)翼的非線性氣動彈性運(yùn)動方程，沒有考慮動態(tài)失速效應(yīng)引起的非線性效應(yīng)。非線性理論梁模型固然具有高精度、高效率等優(yōu)點，但物理變量不直接，通用性受到一定的限制。隨著非線性有限元技術(shù)的發(fā)展，謝長川等[5]基于線性化的UL有限元理論對復(fù)雜結(jié)構(gòu)機(jī)翼的非線性顫振問題開展了大量研究，氣動彈性運(yùn)動方程在頻域內(nèi)求解，而無法預(yù)測臨界速度前后的結(jié)構(gòu)運(yùn)動行為；在時域內(nèi)采用TL法或UL法求解結(jié)構(gòu)動力學(xué)方程時為了獲得較好的精度一般選擇非常小的時間步長，從而顯著地增加求解時間；隨著CR有限元技術(shù)的發(fā)展，其靜力學(xué)求解允許適度大的載荷步、動力學(xué)方程求解同樣允許適度大的時間步，從而引起了許多學(xué)者基于CR列式法對有限旋轉(zhuǎn)問題的研究[7-13]。作者在大柔性太陽能無人機(jī)的靜氣動彈性研究中，引入了CR非線性有限元梁模型描述機(jī)翼的幾何大變形，減少了流固耦合迭代次數(shù)，取得了較好的效果[14]。

本文在所編寫的大柔性結(jié)構(gòu)幾何非線性靜力學(xué)求解程序的基礎(chǔ)上，進(jìn)一步推導(dǎo)了空間CR梁單元的結(jié)構(gòu)動力學(xué)方程，首先以懸臂梁模型為例，研究了大柔性結(jié)構(gòu)在外載荷作用下的強(qiáng)迫振動，驗證了基于CR有限元理論所建立的結(jié)構(gòu)動力學(xué)模型的精度、效率及所編寫程序的正確性。其次，以簡諧振蕩的NACA0012翼型為例，采用改進(jìn)的ONERA動態(tài)失速模型對典型翼段的非線性非定常氣動力進(jìn)行求解，驗證了非定常非線性氣動力模型對動態(tài)失速效應(yīng)模擬的精度。最后，耦合上述結(jié)構(gòu)動力學(xué)模型和非線性非定常氣動力模型，基于松耦合求解技術(shù)，提出了一種適用于幾何大變形大柔性機(jī)翼的非線性氣動彈性求解方法，并研究了大柔性機(jī)翼的極限環(huán)顫振特性以及彈性軸站位對極限環(huán)顫振的影響等。

1大柔性機(jī)翼結(jié)構(gòu)動力學(xué)建模

大柔性機(jī)翼在氣動載荷的作用下，變形前后的幾何位形具有明顯的差異，常規(guī)的線彈性小變形假設(shè)不再成立，而局部結(jié)構(gòu)的材料應(yīng)力應(yīng)變關(guān)系仍處于線彈性范圍內(nèi)。CR有限元法的主要思想是將結(jié)構(gòu)的變形分為剛體的旋轉(zhuǎn)及平移和單元坐標(biāo)系內(nèi)的線彈性變形，平衡方程在單元坐標(biāo)系內(nèi)建立，通過坐標(biāo)轉(zhuǎn)換矩陣及其微分形式構(gòu)造總體坐標(biāo)系下的切線剛度矩陣，因而其核心是坐標(biāo)轉(zhuǎn)換矩陣的建立及其微分形式的推導(dǎo)，其計算精度主要取決于剛度矩陣和坐標(biāo)系轉(zhuǎn)換矩陣對幾何非線性描述的精確性。

1.1CR梁單元靜力學(xué)平衡方程

CR有限元理論的平衡方程建立在單元坐標(biāo)系內(nèi)，通過總體坐標(biāo)系Ixyz、節(jié)點坐標(biāo)系ui和單元平均構(gòu)形坐標(biāo)系um建立單元坐標(biāo)系ue，見圖1。

結(jié)構(gòu)位形確定后，單元節(jié)點的三個歐拉角和現(xiàn)時坐標(biāo)隨之確定，按照式(1)旋轉(zhuǎn)總體坐標(biāo)系得到節(jié)點坐標(biāo)系：

ui(x)=Ti(x)Ixyz

(1)

那么參照文獻(xiàn)[9,14]，CR空間梁單元的單元坐標(biāo)系定義為：

ue1=(d2-d1)/ln

ue2=um2-((uTm2ue1)/2)(ue1+um1)

ue3=um3-((uTm3ue1)/2)(ue1+um1)

(2)

式中，d1與d2分別是節(jié)點i,j的現(xiàn)時坐標(biāo)ln=[(d2-d1)T(d2-d1)]1/2。

圖1　柔性結(jié)構(gòu)變形示意圖 Fig.1 Sketch of a deformed very flexible structure

可見節(jié)點坐標(biāo)系和單元坐標(biāo)系都是結(jié)構(gòu)位形的函數(shù)，在求解過程中需要隨變形的增加而不斷更新。

CR有限元法充分利用了幾何非線性問題的小應(yīng)變、大位移特征，彈性變形在單元坐標(biāo)系內(nèi)描述，位移-應(yīng)變關(guān)系仍然是線性的。利用方程2建立的單元坐標(biāo)系與方程1建立的節(jié)點坐標(biāo)系，得到單元坐標(biāo)系下的節(jié)點位移：

ul=ln-lo

2θl1,1=-uTi3ue2+uTi2ue3

2θl1,2=-uTi1ue3+uTi3ue1

2θl1,3=-uTi2ue1+uTi1ue2

2θl2,1=-uTi3ue2+uTi2ue3

2θl2,2=-uTi1ue3+uTi3ue1

2θl2,3=-uTi2ue1+uTi1ue2

(3)

對式(3)進(jìn)行微分得到單元坐標(biāo)系下位移增量δpl到總體坐標(biāo)系下位移增量δp的坐標(biāo)轉(zhuǎn)換矩陣T：

δpl=Tδp

(4)

總體坐標(biāo)系下廣義節(jié)點力向量f為：

f=TTfl=TTKlpl

(5)

δf=δTTKlpl+TTKlδpl=

[KTσ(fl)+TTKlT]δp

(6)

增量平衡方程為：

KTΔp=ΔFe

(7)

KT=KTσ(fl)+TTKlT為所求的切線剛度矩陣，其中TTKlT為材料剛度矩陣，KTσ(fl)為幾何剛度矩陣，是單元坐標(biāo)系下單元節(jié)點力fl的函數(shù)。

當(dāng)大柔性機(jī)翼上布置有發(fā)動機(jī)等推進(jìn)系統(tǒng)時，將會引入隨動力進(jìn)而影響機(jī)翼的顫振穩(wěn)定性[17]。假設(shè)該發(fā)動機(jī)掛載在第i個節(jié)點下，F(xiàn)s(t)表示總體坐標(biāo)系下隨動載荷等效到節(jié)點i處的外載荷。那么，第i個節(jié)點的總外載荷可以描述為：

Fi,e(t)=Fi,a(t)+Fs(t)

(8)

當(dāng)隨動載荷描述在局部坐標(biāo)系下時，需要把Fs,i(t)轉(zhuǎn)換到總體坐標(biāo)系下，轉(zhuǎn)換表達(dá)式為：

Fs(t)=Ui,IFs,i(t)

(9)

具體求解時需要注意，隨結(jié)構(gòu)運(yùn)動需要實時更新轉(zhuǎn)換矩陣Ui,I。

1.2CR梁單元動力學(xué)平衡方程

大柔性機(jī)翼在變形的過程中可以認(rèn)為其局部翼剖面沒有發(fā)生翹曲[4，6]，則剖面內(nèi)所有點具有相同的角速度，截面線動量與角動量為：

(10)

那么，截面慣性力和慣性力矩可以表示為：

(11)

截面剛心處的線加速度、角速度和角加速度依次通過節(jié)點線加速度、角速度和角加速度插值得到：

(12)

同樣可以插值得到截面虛位移，那么慣性力在虛位移上所做的虛功為：

(13)

根據(jù)虛功原理，等效節(jié)點慣性力所做的虛功為：

δW=f Tmass(δuTi,δθTi,δuTj,δθTj)T

(14)

把方程13代入式(14)得：

由方程(15)可以看出，線加速度是定義在總體坐標(biāo)系下的，而角加速度是定義在單元坐標(biāo)系下的，并且慣性力與時間和結(jié)構(gòu)位形有關(guān)。

根據(jù)方程(15)所得到的慣性力表達(dá)式，忽略阻尼后考慮幾何非線性效應(yīng)的大柔性結(jié)構(gòu)動力學(xué)平衡方程可以表示為：

Fe-fmass-fi=0

(16)

式中，F(xiàn)e表示結(jié)構(gòu)t時刻所受的外載荷，fmass表示結(jié)構(gòu)t時刻的慣性力，fi表示結(jié)構(gòu)t時刻的內(nèi)力，它們均與時間t有關(guān)。

1.3非線性動力學(xué)平衡方程的時域推進(jìn)求解

采用隱式Newmark時域積分格式求解大柔性機(jī)翼的結(jié)構(gòu)動力學(xué)問題時，需要考慮系統(tǒng)在tn+1時的平衡，并采用Newton-Raphson法迭代求解控制平衡方程。Newmark積分法可以認(rèn)為是對線性加速度法的推廣，其插值假設(shè)是：

(17)

式中，ΔΨ=UTnΔψ，為隨動坐標(biāo)系下描述節(jié)點坐標(biāo)系Un旋轉(zhuǎn)到Un+1所對應(yīng)的旋轉(zhuǎn)偽矢量；Δψ為總體坐標(biāo)系下描述節(jié)點坐標(biāo)系Un旋轉(zhuǎn)到Un+1所對應(yīng)的旋轉(zhuǎn)偽矢量。

根據(jù)上述假設(shè)，對方程15進(jìn)行微分，得：

δfmass=δ([Ue]fmass,in)=MT,massδp

(18)

根據(jù)方程(18)得到慣性力的切線表達(dá)形式，可以采用預(yù)估校正技術(shù)求解方程(16)。預(yù)估的本質(zhì)是用tn時刻位移及增量斜率預(yù)測tn+1時刻的位移。第一次位移增量Δp0可以通過下式預(yù)估求解得到：

(19)

將預(yù)估求解得到的位移增量代入方程式(16)，計算tn+1時刻第k次迭代后的非平衡力：

(20)

校正步內(nèi)采用下式迭代求解，直到滿足收斂標(biāo)準(zhǔn)：

gk=(KkT,tn+1+MkT,tn+1)Δpk

(21)

與靜力學(xué)增量平衡方程的求解基本一致，采用預(yù)估校正推進(jìn)技術(shù)求解方程式(16)時，在校正步內(nèi)一般采用Newton-Raphson或修正的Newton-Raphson迭代求解，只是動力學(xué)求解中殘差gk隨所求解的時間步內(nèi)位移增量的更新而不斷更新。

2非線性非定常氣動力模型

采用ONERA動態(tài)失速模型描述非線性非定常氣動力，在不同文獻(xiàn)中，其表達(dá)式形式也不盡相同；這里采用改進(jìn)的ONEAR(EDLin)動態(tài)失速模型[18]。非定常升力、力矩作用在四分之一弦長處，表達(dá)式為：

(22)

(23)

(24)

把式(23)和(24)寫成緊湊矩陣形式：

(25)

(26)

把式(26)寫成緊湊矩陣形式：

采用四級四階Runge-Kutta公式求解改進(jìn)的ONERA動態(tài)失速模型中的微分方程，求解格式為：

(27)

3非線性氣動彈性運(yùn)動方程

第二節(jié)中的動力學(xué)方程求解是以總體坐標(biāo)系下的節(jié)點線位移、角位移、線速度、線加速度和單元坐標(biāo)系下角速度、角加速度為未知量的；非線性非定常氣動力的求解輸入是局部氣流坐標(biāo)系下的片條沉浮線位移、線速度、線加速度、角位移、角速度、角加速度等，因此需要對動力學(xué)模型的解向量進(jìn)行轉(zhuǎn)換和插值。

首先繞Iy旋轉(zhuǎn)90°，再繞現(xiàn)在的x軸旋轉(zhuǎn)90°得到總體氣流坐標(biāo)系，再繞總體氣流坐標(biāo)系的負(fù)x軸旋轉(zhuǎn)θ得到局部氣流坐標(biāo)系：

(28)

那么總體坐標(biāo)系下的線位移、角位移、線加速度、角加速度和單元坐標(biāo)系下的角速度及角加速度到局部氣流坐標(biāo)系下的轉(zhuǎn)換關(guān)系式為：

(29)

經(jīng)過插值得到每個片條的振動及沉浮、扭轉(zhuǎn)運(yùn)動：

(30)

根據(jù)式(19)和(23)計算得到局部氣流坐標(biāo)系下的片條非定常氣動力Qs，然后根據(jù)功互等原理，將片條在局部氣流坐標(biāo)系下的非定常氣動載荷按照等效節(jié)點力的形式插值到相鄰節(jié)點上，計算公式為：

(31)

采用下式把局部氣流坐標(biāo)系下的節(jié)點力轉(zhuǎn)換到總體坐標(biāo)系下：

(32)

對式(32)得到的節(jié)點氣動力組裝后可以得到總體坐標(biāo)系下的非定常氣動力Qe，代入方程(16)代替外載荷向量Fe，可以得到機(jī)翼的氣動彈性運(yùn)動方程：

Qe-fmass-fi=0

(33)

在時域內(nèi)對式(33)進(jìn)行直接數(shù)值積分求解時，與結(jié)構(gòu)動力學(xué)求解中tn+1時刻的外載荷已知不同，tn+1時刻的非定常氣動力Qe與tn+1時刻的結(jié)構(gòu)運(yùn)動狀態(tài)有關(guān)，因而需要預(yù)估tn+1時刻的非定常氣動力，然后根據(jù)解得的tn+1時刻的結(jié)構(gòu)狀態(tài)更新非定常氣動力，根據(jù)Newton等距向前插值多項式進(jìn)行外推：

Qe,tn+1=3Qe,tn-3Qe,tn-1+Qe,tn-2

(34)

松耦合的求解思想是：首先預(yù)估下一時刻的非定常氣動力，然后帶入氣動彈性運(yùn)動方程中求解結(jié)構(gòu)運(yùn)動狀態(tài)，以解得的下一時刻結(jié)構(gòu)狀態(tài)重新計算非定常氣動力代替預(yù)估值，但不參與迭代求解，而是直接推進(jìn)到下一時間步的求解；若參與迭代求解，上述計算方法則變成緊耦合推進(jìn)求解算法。緊耦合迭代求解算法具有精度高的優(yōu)點但求解易發(fā)散，一般很少在非線性氣動彈性系統(tǒng)的求解中使用。綜上所述，松耦合推進(jìn)求解方程式(33)的流程見圖2。

圖2　非線性氣動彈性系統(tǒng)松耦合求解流程 Fig.2 Calculation process of the nonlinear aeroelastic system based on loosely coupled technology

4算例及驗證

大柔性機(jī)翼非線性氣動彈性求解所需的幾何參數(shù)、剛度參數(shù)、質(zhì)量參數(shù)及大氣參數(shù)見表1。

4.1大柔性懸臂梁幾何非線性動力學(xué)求解

首先把大柔性機(jī)翼作為懸臂梁，不考慮氣動力的作用，在翼尖施加集中載荷歷程見圖3(b)，分別采用經(jīng)典UL法(時間步長為1E-4s)和CR法(時間步長為5E-2s)進(jìn)行求解，以驗證基于CR有限元法采用直接數(shù)值積分技術(shù)(Newmark法)求解大柔性結(jié)構(gòu)動力學(xué)響應(yīng)的精度等。

表1　大柔性機(jī)翼模型參數(shù) [4]

圖3　大柔性懸臂梁及其端部加載示意圖 Fig.3 Sketch of a very flexible cantilever beam and the tip loading history

圖4　不同步長下CR法與UL法求解結(jié)果對比 Fig.4 Results comparing between CR method and UL method under variant time steps

由圖4中的計算結(jié)果可以看出，采用CR法求解大柔性結(jié)構(gòu)的動力學(xué)響應(yīng)問題時允許采用較大的時間步長，且具有較好的求解精度，滿足大柔性機(jī)翼結(jié)構(gòu)動力學(xué)求解對幾何非線性因素的描述精度要求。

4.2非線性非定常氣動力模型驗證

以NACA0012翼型為例，采用ONERA動態(tài)失速模型及其改進(jìn)形式對非定常非線性氣動力進(jìn)行求解，結(jié)果見圖5。

圖5　α=15°+10° cos(0.1τ)的非定常氣動力 Fig.5 Unsteady aerodynamics atα=15°+10° cos(0.1τ)

大迎角區(qū)域內(nèi)，經(jīng)典的ONERA(EDLin)模型求解精度較好，而改進(jìn)后的ONERA模型雖然略微降低了對大迎角區(qū)的模擬精度但是顯著的改善了小迎角區(qū)的求解精度。

仍然以NACA0012翼型為例，平均迎角分別取6°、12°、17°，h=0.1sin(0.1τ)m，τ=Vb/t；采用改進(jìn)的ONERA動態(tài)失速模型研究完全線性段、部分線性段部分非線性段及完全非線性段內(nèi)的非定常氣動力特性，得到翼段純沉浮運(yùn)動下的非定常氣動力見圖6。

圖6　h=0.1sin(0.1τ)m的非定常氣動力 Fig.6 Unsteady aerodynamics at h=0.1 sin (0.1τ)m

可見采用改進(jìn)的ONERA模型無論是線性段還是非線性段均可以較好的模擬二維機(jī)翼的非定常氣動力，特別是非線性段的動態(tài)失速效應(yīng)。

4.3大柔性機(jī)翼非線性氣動彈性響應(yīng)求解

以表1中所給出的大柔性機(jī)翼為例，機(jī)翼剖面為NACA0012翼型，初始翼尖彎曲位移為4m時，本文與文獻(xiàn)4及文獻(xiàn)6得到的臨界極限環(huán)顫振速度均約為28m/s，翼尖位移響應(yīng)歷程與文獻(xiàn)值的對比見圖7示。

圖7　V=28m/s，y 0=4.0m時翼尖彎曲位移響應(yīng) Fig.7 Wingtip displacement response at V=28m/s，y 0=4.0m

從圖7可以看出，采用所提出的非線性氣動彈性求解方法解得的位移響應(yīng)與文獻(xiàn)值6吻合的更好，可能是因為文獻(xiàn)4中沒有考慮非定常氣動力的動態(tài)失速效應(yīng)。翼尖位移響應(yīng)的相位圖見圖8。

圖8 V=28m/s,y0=4.0m翼尖相位圖Fig.8Phase-planeplotsofthewingtipresponseatV=28m/s,y0=4.0m圖9 y0=4.0m時不同來流速度下的翼尖彎曲位移響應(yīng)Fig.9Wingtipdisplacementresponseaty0=4.0mandvariantairspeeds

圖10　y 0=4.0m時不同來流速度下的翼尖響應(yīng)相位圖 Fig.10 Phase-plane plots of the wingtip response at y 0=4.0m and variant air speeds

由計算結(jié)果可以看出，采用所提出的非線性氣動彈性求解方法可以較好的預(yù)測大柔性機(jī)翼的極限環(huán)顫振。進(jìn)一步增加來流速度，研究來流速度大于臨界速度時的大柔性機(jī)翼非線性氣動彈性響應(yīng)特性；初始彎曲位移為4.0m。

由計算結(jié)果可以看出，當(dāng)來流速度略大于28m/s時，大柔性機(jī)翼的氣動彈性翼尖位移響應(yīng)仍具有周期性，當(dāng)來流速度增加到32m/s時，翼尖位移響應(yīng)不再具有周期性。隨來流速度的增加，翼尖彎曲位移響應(yīng)及扭轉(zhuǎn)響應(yīng)幅值隨之增加，且極限環(huán)顫振變得越來越復(fù)雜。來流速度為28m/s時，極限環(huán)顫振每個周期的振動幅值一致，但來流速度增加后，翼尖位移響應(yīng)需要經(jīng)過一個幅值較小的周期和一個幅值較大的周期才能夠回到原位置。另外，從圖10中還可以看出隨來流速度的增加，局部動態(tài)失速效應(yīng)越來越明顯。

4.4彈性軸站位對大柔性機(jī)翼極限環(huán)顫振的影響

給定初始條件下，誘發(fā)非線性氣動彈性極限環(huán)顫振的最小速度稱為臨界極限環(huán)顫振速度[4]。結(jié)構(gòu)參數(shù)仍以表一中給出的參數(shù)作為基準(zhǔn)不變，研究彈性軸和剖面質(zhì)心一起向前緣移動時對臨界極限環(huán)顫振速度的影響。初始翼尖彎曲位移分別選為4m、3m、2m、1m，得到的臨界極限環(huán)顫振速度隨彈性軸在弦向站位的不同而增減的趨勢見圖11。

圖11　彈性軸站位對極限環(huán)顫振的影響 Fig.11 The effects on limit cycle oscillation due to variant elastic axis positions

從圖11中可以看出，翼尖彎曲位移的增加將降低機(jī)翼的臨界極限環(huán)顫振速度；彈性軸在50%的弦長處時，初始翼尖彎曲位移為4m時解得的臨界極限環(huán)顫振速度與初始翼尖彎曲位移為1m時相比降低了15.36%；另外，彈性軸的前移可以有效的提高臨界極限環(huán)顫振速度；因此，對于這類大柔性機(jī)翼，通過合理設(shè)計彈性軸的弦向站位，可以有效的改善彎曲變形對機(jī)翼氣動彈性所引起的不利效應(yīng)。

4.5積分步長的選擇及收斂性

直接數(shù)值積分的收斂性和精度可以通過不同時間步長的選取以證明，在上述研究中時間步長均為0.005s，這里進(jìn)一步補(bǔ)充來流速度為28m/s時，初始彎曲位移為4m，時間步長為0.01s和0.0001s的計算結(jié)果。

圖12　不同時間步長的積分結(jié)果 Fig.12 Comparing of direct integration results between variant time steps

從圖12中可以看出，時間步長為0.005s和0.0001s的積分結(jié)果基本一致，因此時間步長選為0.005s滿足計算精度的要求，并且與時間步長為0.0001s時相比大大減少了求解時間。針對具體的求解問題，對比不同的積分步長對求解精度和穩(wěn)定性的影響，指導(dǎo)選擇合適的積分步長，可以有效地提高求解效率；具體執(zhí)行時，可選擇下一次積分步長為本次積分步長的1/N(N為大于等于2的整數(shù))，兩次求解得到的最大響應(yīng)幅值相差小于幅值的5%時，即可把本次試算積分步長作為采用直接數(shù)值積分求解目標(biāo)算例的較為合理的時間積分步長。

5結(jié)論

本文基于CR有限元非線性動力學(xué)求解技術(shù)和改進(jìn)的ONERA動態(tài)失速非線性非定常氣動力模型提出了一種適用于大柔性大變形機(jī)翼的非線性氣動彈性求解方法，進(jìn)而對大柔性機(jī)翼的非線性氣動彈性響應(yīng)特性進(jìn)行了較為深入的研究，主要得到以下結(jié)論：

(1)本文所提出的非線性氣動彈性響應(yīng)計算方法，以總體坐標(biāo)系下的位移等為變量，物理意義簡單直接，具有較好的求解精度和計算效率，能夠滿足大柔性機(jī)翼非線性氣動彈性響應(yīng)求解的需求，并可以較好的預(yù)測大柔性機(jī)翼的極限環(huán)顫振。

(2)較小的時間步長可以提高求解的精度但計算花費(fèi)隨之增加，在具體的求解中可以預(yù)先對比不同積分步長對求解精度的影響，指導(dǎo)選取合適的積分步長，以滿足求解的需求。

(3)來流速度大于臨界極限環(huán)顫振速度時，大柔性機(jī)翼的非線性氣動彈性響應(yīng)是非常復(fù)雜的，非定常氣動力的動態(tài)失速效應(yīng)雖然限制了振幅的增加，但也會帶來結(jié)構(gòu)疲勞等問題，在這類機(jī)翼的結(jié)構(gòu)設(shè)計時需要特別注意。

(4)前移彈性軸可以有效地提高臨界極限環(huán)顫振速度，增加周期性極限環(huán)顫振的速度區(qū)間，顯著地改善大柔性機(jī)翼的非線性氣動彈性特性，可以用于指導(dǎo)大柔性機(jī)翼彈性軸在弦向的站位設(shè)計。進(jìn)一步的研究工作中可以引入高效的預(yù)估求解算子，減少動力學(xué)響應(yīng)求解時的迭代次數(shù)，提高求解效率。

參考文獻(xiàn)

[1]Patil M J, Hodges D H. On the importance of aerodynamic and structural geometrical nonlinearities in aeroelastic behavior of high-aspect-ration wings[J].Journal of Fluids and Structures,2004,19:905-915.

[2]Patil M J, Hodges D H, Cesnik C E S. Nonlinearaeroelasticity and flight dynamics of high-altitude long-endurance Aircraft [R]. AIAA-99-1470.

[3]Patil M J, Hodges D H. Flightdynamics of highly flexible flying wings [J]. Journal of Aircraft，2006，43(6):1790-1798.

[4]Patil M J, Hodges D J, Cesnik C E S. Limit cycle oscillations in high-aspect-ratio wings[J].Journal of Fluid and Structures,2001,15:107-132.

[5]Xie Chang-chuan, Yang Chao. Linearization methods of nonlinear aeroelastic stability for complete aircraft with high-aspect-ratio wings[J]. Sci China Tech Sci,2011, 54:403-411.

[6]Shams S, Sadr M H, Haddadpour H.An efficient method for nonlinear aeroelasticy of slender wings[J].Nonlinear Dynamics,2012,67:659-681.

[7]周凌遠(yuǎn)，李喬.基于UL法的CR列式三維梁單元計算方法[J]. 西南交通大學(xué)學(xué)報，2006，41(6):690-695.

ZHOU Ling-yuan, LI Qiao. Updated lagrangian co-rotational formulation for geometrically nonlinear FE analysis of 3D beam element[J].Journal of Southwest Jiaotong Unoversity,2006,41(6):690-695.

[8]Belytschko T, Schwer L. Large displacement, transient analysis of space frames[J]. International Journal for Numerical Methods in Engineering, 1977, 11:65-84.

[9]Crisfield M A. A consistent Co-rotational formulation for non-linear, three-dimensional, beam element[J]. Computer Methods In Applied Mechanics And Engineering, 1990,81:131-150.

[10]Crisfield M A. Non-linear finite element analysis of solids and structures, Volume 2: Advanced topics[M]. John Wiley & Sons, Chichester, New York, Weinheim, Brisbane, Singapore, Toronto,2000.

[11]Crisfield M A,Galvanetto U, Jelenicé G. Dynamics of 3-D co-rotational beams[J]. Computational Mechanics, 1997,20：507-519.

[12]Ghosh S, Roy D. Consistent quaternion interpolation for objective finite element approximation of geometrically exact beams[J]. Computer Methods In Applied Mechanics And Engineering,2008, 198:555-571.

[13]Le T N, Battini J M, Hjiaj M. Dynamics of 3d beam elements in a corotational context: a comparative study of established and new formulation[J]. Finite Elements in Analysis and Design,2012,61:97-111.

[14]王偉,周洲,祝小平,等.考慮幾何非線性效應(yīng)的大柔性太陽能無人機(jī)靜氣動彈性分析[J]. 西北工業(yè)大學(xué)學(xué)報, 2014, 32(4):499-504.

WANG Wei, ZHOU Zhou, ZHU Xiao-ping,et al. Static aeroelastic characteristics analysis of a very flexible solar powered UAV with geometrical nonlinear effect considered[J]. Journal of Northwestern Polytechnical University, 2014, 32(4):499-504.

[15]Palacios R, Murua J, Cook R. Structural and aerodynamic models in nonlinear flight dynamics of very flexible aircraft [J]. AIAA Journal, 2010,48 (11):2648-2659.

[16]吳國榮,顏桂云，陳福全.柔性梁大柔度動力響應(yīng)分析的多體系統(tǒng)方法[J].振動與沖擊,2007,26(3):87-89.

WU Guo-rong,YAN Gui-yun,CHEN Fu-quan.Large deflection dynamic response analysis of flexible beams bi multibody system method[J].Journal of Vibration and Shock, 2007,26(3):87-89.

[17]張健,向錦武.側(cè)向隨動力作用下大展弦比柔性機(jī)翼的穩(wěn)定性[J].航空學(xué)報,2010,31(11):2115-2123.

ZHANG Jian, XIANG Jin-wu. Stability of high-aspect-ratio fiexible wings loaded by a lateral follower force[J].Chinese Journal of Aeronautics,2010,31(11):2115-2123.

[18]Laxman V,Venkatesan C.Chaotic response of an airfoil due to aeroelastic coupling and dynamic stall[J].AIAA Jourmal, 2007,45(1):271-280.

第一作者何洪陽男，碩士，1988年生

通信作者陳春俊男，博士，教授，1967年生