郭同彪，白俊強, *，孫智偉，王晨

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

考慮機翼幾何非線性的氣動彈性建模與分析

郭同彪1，白俊強1, *，孫智偉2，王晨1

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

大展弦比柔性機翼結(jié)構(gòu)重量輕、氣動效率高，廣泛應(yīng)用于高空長航時無人機(UAVs)。飛行過程中，這類機翼在氣動力作用下發(fā)生大變形，線性結(jié)構(gòu)模型不再適用，需要建立考慮幾何大變形的結(jié)構(gòu)模型。采用牛頓力學方法推導(dǎo)了考慮結(jié)構(gòu)幾何非線性的機翼結(jié)構(gòu)動力學模型，該方法推導(dǎo)過程簡潔、物理意義明確，可以與Hodges基于哈密頓原理的推導(dǎo)方法相互補充，相互驗證。為了能夠更準確地求解大展弦比柔性機翼的非定常氣動力，建立了能夠考慮機翼三維效應(yīng)且適用于機翼空間大變形的非定常氣動力模型?；诮⒌姆蔷€性結(jié)構(gòu)模型和非定常氣動力模型，采用松耦合方法建立了非線性氣動彈性模型，并通過算例驗證了氣彈模型的準確性。研究結(jié)果表明，大展弦比柔性機翼顫振速度對來流迎角和機翼的展長均較為敏感；當來流速度大于顫振速度時，由于幾何非線性，機翼振動并未發(fā)散而是形成穩(wěn)定的極限環(huán)振蕩(LCO)；隨著來流速度進一步增加，機翼再次穿過臨界穩(wěn)定點，由不穩(wěn)定系統(tǒng)變?yōu)榉€(wěn)定系統(tǒng)，直到隨著速度的增加系統(tǒng)再次達到臨界穩(wěn)定狀態(tài)。

非線性氣動彈性；幾何非線性；柔性機翼；顫振速度；極限環(huán)振蕩；時域響應(yīng)

高空長航時無人機在情報、偵查、監(jiān)視等軍用領(lǐng)域以及環(huán)境監(jiān)測、通信中繼、氣象等民用領(lǐng)域，有著廣闊的應(yīng)用前景，近年來在世界各國得到了廣泛的發(fā)展。為了減輕結(jié)構(gòu)重量、提高升阻比，高空長航時無人機普遍采用大展弦比柔性機翼[1]。在飛行過程中特別是受到大氣紊流擾動的情況下，這類機翼產(chǎn)生較大的彎曲和扭轉(zhuǎn)變形，對結(jié)構(gòu)固有模態(tài)和動力學特性有著較大的影響。文獻[2]指出，機翼發(fā)生大變形時其顫振速度和頻率相對于未變形狀態(tài)變化可高達50%，此時基于小變形假設(shè)的線性結(jié)構(gòu)模型無法滿足精度要求，需要建立考慮幾何非線性的結(jié)構(gòu)模型[2-4]。

Hodges建立的幾何非線性梁動力學模型[5]，能夠考慮任意大位移/轉(zhuǎn)動和橫向剪切變形，被廣泛應(yīng)用于大展弦比柔性機翼的結(jié)構(gòu)建模[2,6-8]。Hodges基于哈密頓原理推導(dǎo)梁動力學模型，該方法無需對梁進行受力分析，不足之處在于推導(dǎo)過程比較繁瑣，且物理意義不明確，這也是分析力學方法建立動力學方程的固有缺陷。本文采用經(jīng)典的牛頓方法推導(dǎo)了幾何非線性梁模型，該方法推導(dǎo)過程簡潔，物理意義明確，可以與Hodges的推導(dǎo)方法相互補充，具有一定的理論意義。

相比于傳統(tǒng)剛度較大的機翼，大展弦比柔性機翼的氣動彈性設(shè)計更具有挑戰(zhàn)性[1]，其在氣動力作用下易發(fā)生幾何大變形，并受大變形影響，發(fā)生氣動失速[4]。這類機翼的設(shè)計難點主要在于幾何大變形，以及由此引發(fā)的結(jié)構(gòu)幾何非線性和氣動力非線性。

在靜氣彈方面，由于幾何大變形，機翼氣動載荷分布[9-10]、配平迎角[11-12]等均發(fā)生了顯著的變化，在設(shè)計初期必須考慮大變形帶來的影響。

相比于靜氣彈，大展弦比柔性機翼的動氣動彈性問題更加復(fù)雜，近年來吸引了大批研究者。文獻[6]指出，對于大展弦比柔性機翼，即使來流速度小于顫振速度，在較大的初始擾動作用下，機翼振蕩也可能不收斂，表現(xiàn)為等幅振蕩。多個學者的研究結(jié)果表明[6-7,13-15]，當來流速度大于顫振速度時，大變形引起的幾何剛化效應(yīng)和氣動分離誘導(dǎo)機翼表現(xiàn)為有限幅值的極限環(huán)振蕩(LCO)，而非一直發(fā)散；而且當不考慮氣動非線性時，僅考慮幾何非線性也能導(dǎo)致機翼呈現(xiàn)極限環(huán)振蕩[16]。

此外，各國研究者研究了結(jié)構(gòu)幾何參數(shù)和來流迎角對大展弦比柔性機翼氣動彈性特性的影響。Patil[8]和楊智春[17]等研究了來流迎角對機翼顫振速度和顫振頻率的影響。Eskandary等[14]研究了質(zhì)量比和剛度比對大展弦比機翼顫振速度和發(fā)散速度的影響。Suleman等[18]研究了機翼弦長和梢根比對大展弦比機翼時域響應(yīng)的影響。以上研究中，非定常氣動力模型大多采用二元翼氣動力和片條理論進行建模，無法考慮氣動力的三維效應(yīng)。

建立能夠考慮三維效應(yīng)的非定常氣動力模型是柔性機翼非線性氣動彈性建模的一個難點。經(jīng)典氣彈模型中三維氣動力采用偶極子網(wǎng)格法進行求解；當結(jié)構(gòu)發(fā)生變形時，偶極子網(wǎng)格無法跟隨結(jié)構(gòu)運動形成曲面，不能準確求解機翼發(fā)生大變形時的氣動載荷。楊超等[19-20]采用曲面渦格法建立了柔性機翼和飛機的氣動彈性模型，它能夠求解任意曲面變形機翼的氣動力，且氣動網(wǎng)格隨著結(jié)構(gòu)變形進行更新，適用于大變形機翼的氣動彈性建模。

非定常渦格法[21](Unsteady Vortex-Lattice Method，UVLM)是一種時域非定常氣動力模型，氣動網(wǎng)格可布置成空間任意曲面形狀，且可實現(xiàn)氣動面隨結(jié)構(gòu)變形的自適應(yīng)更新，適用于大變形機翼的定常和非定常氣動力求解。本文耦合該氣動力模型和結(jié)構(gòu)模型建立了考慮結(jié)構(gòu)幾何非線性的機翼動氣動彈性模型，并基于該模型研究了大展弦比柔性機翼的動氣動彈性特性。

1 大展弦比機翼梁結(jié)構(gòu)動力學建模

大展弦比柔性機翼在氣動力作用下會發(fā)生大變形。機翼各個剖面不僅發(fā)生大的平移，還會產(chǎn)生較大的轉(zhuǎn)動。為了準確地對機翼梁結(jié)構(gòu)的變形進行描述，需引入4類坐標系(見圖1)：慣性坐標系I、未變形梁截面坐標系b、變形梁截面坐標系B以及機翼隨體坐標系A(chǔ)。

未變形梁、變形梁的參考線和橫截面如圖2所示，本文在推導(dǎo)動力學方程時以剛心軸作為梁的參考線。x1為沿未變形梁參考線的曲線坐標，s為沿變形梁參考線的曲線坐標，且s=s(x1)。未變形梁參考線上任一點在慣性系下的位置矢量用r(x1)表示，位移矢量用u(x1)表示，則變形梁參考線上該點在慣性系下的位置矢量R(x1)=r(x1)+u(x1)。根據(jù)坐標系的定義可知，未變形梁參考線的單位切線矢量為r′=b1，z′表示任意變量z對未變形梁曲線坐標x1求偏導(dǎo)；變形梁參考線的單位切線矢量為?R/?s=R′/s′,s′為基于未變形梁參考線的拉伸應(yīng)變。

在考慮大變形的預(yù)扭/彎梁模型中，未變形梁曲率矢量k和變形梁曲率矢量K是兩個關(guān)鍵的幾何參數(shù)，分別表示在b坐標系和B坐標系下，即

{ki=k·bi

Ki=K·Bi

(1)

圖1 本征梁模型相關(guān)的坐標系Fig.1 Coordinate system of intrinsic beam model

圖2 本征梁變形示意圖Fig.2 Schematic of deformation of intrinsic beam

式中：i=1,2,3；k1為未變形梁單位長度的扭轉(zhuǎn)(即預(yù)扭)；k2和k3分別為未變形梁在b2、b3方向的曲率；K1為變形梁單位長度的扭轉(zhuǎn)；K2和K3分別為變形梁在B2、B3方向的曲率。

梁截面坐標系與梁曲率矢量有如下關(guān)系[22]：

(2)

根據(jù)式(2)，任意矢量z在A坐標系和B坐標系下對x1的偏導(dǎo)數(shù)有如下關(guān)系：

(3)

下面針對梁截面微元采用動量定理和動量矩定理推導(dǎo)幾何非線性梁的動力學方程。圖3為梁截面微元O1O2在平面內(nèi)的放大示意圖，其中，微元長度為dx1,T坐標系為忽略剪切變形時梁截面的局部坐標系；B坐標系基矢量和T坐標系基矢量方向的差別反映了梁剖面的剪切應(yīng)變，2γ12和2γ13為剪切應(yīng)變，與拉伸應(yīng)變γ11構(gòu)成梁的廣義力應(yīng)變。下面基于該梁微元對梁進行受力分析，并建立梁的動力學方程。

圖3 梁微元變形放大示意圖Fig.3 Schematic of infinitesimal deformation of beam

1.1 基于動量定理推導(dǎo)幾何非線性梁的動力學方程

設(shè)微元左端受到的彈性力為F(x1)，則在右端受到的彈性力可表示為

(4)

則單位長度梁元受到的彈性合力為

(5)

令梁橫截面內(nèi)任意質(zhì)點P在局部坐標系B內(nèi)的位置矢量為

ξP=[0xP2xP3]

(6)

則該質(zhì)點的速度矢量為

VP=VB+Ω×ξP

(7)

式中：VB為梁截面B坐標系原點的速度；Ω為梁截面的角速度矢量。

單位長度梁元的線動量為

μ(VB+Ω×ξcg)

(8)

式中：ξcg為梁截面質(zhì)心的位置矢量；μ為單位長度梁元的質(zhì)量。

單位長度梁元的線動量對時間的導(dǎo)數(shù)為

I(dP/dt)=B(dP/dt)+Ω×P

(9)

根據(jù)動量定理，質(zhì)點系的動量對時間的導(dǎo)數(shù)等于作用在該質(zhì)點系的所有外力的矢量和，即

(10)

1.2 基于動量定理推導(dǎo)幾何非線性梁的力矩方程

設(shè)微元左端的彈性力矩為M(x1)，則在右端的彈性力矩可表示為

M(x1+dx1)=M(x1)+A(?M/?x1)dx1=

(11)

則單位長度梁元受到的彈性合力矩為

(12)

[(1 +γ11) 2γ122γ13]Tdx1=(e1+γB)dx1

(13)

單位長度梁元對慣性系原點OI的動量矩為

(14)

(15)

根據(jù)動量矩定理可知，質(zhì)點系對固定點的動量矩對時間的導(dǎo)數(shù)等于作用在該質(zhì)點系的所有外力對同一點的矩的矢量和，則

(16)

將式(10)代入式(16)，得到

(17)

式(17)即描述梁轉(zhuǎn)動的動力學方程。至此完成了幾何非線性梁動力學方程的推導(dǎo)。將平動和轉(zhuǎn)動動力學方程統(tǒng)一表示為

(18)

(19)

可以看出，以上推導(dǎo)具有以下幾個特點：① 推導(dǎo)過程十分簡潔、直觀；② 直接得到了動力學方程的強形式；③ 動力學方程中每一項物理意義明確。

1.3 結(jié)構(gòu)動力學方程的空間離散

幾何非線性梁的動力學方程是對時間和空間求導(dǎo)的偏微分方程組，求解困難。為了便于求解，采用有限元方法進行空間離散，將偏微分方程組轉(zhuǎn)化為時間坐標的常微分方程組。選取R(s)(見圖2)和Ψ(s)(從A坐標系到B坐標系的笛卡兒旋轉(zhuǎn)矢量)作為結(jié)構(gòu)狀態(tài)變量[23]。第n個有限元內(nèi)任意點的結(jié)構(gòu)狀態(tài)變量為

(20)

式中：Ni(s)為形函數(shù)，sn-1≤s≤sn。令η為所有有限元節(jié)點的位置和旋轉(zhuǎn)矢量，則對機翼結(jié)構(gòu)動力學方程進行有限元離散后表述為

(21)

式中：M為切線質(zhì)量矩陣(Tangent Mass Matrix)；Qgyr、Qstiff、Qext分別為回轉(zhuǎn)力、彈性力和外力，詳細推導(dǎo)過程見文獻[23]。

2 大展弦比柔性機翼氣動力求解

本文采用非定常渦格法(UVLM)計算大展弦比柔性機翼的氣動力。該算法在每個離散時間步處理翼面網(wǎng)格和尾跡網(wǎng)格的運動以及附著渦和尾跡渦的變化，適用于大變形機翼氣動力求解。非定常渦格法的求解流程如圖4所示。

圖4 非定常渦格法計算流程Fig.4 Flow chart for unsteady vortex-lattice method calculation

圖5 渦格劃分示意圖Fig.5 Schematics of vortex lattice discretization

圖5為渦格法機翼網(wǎng)格劃分示意圖。與偶極子格網(wǎng)法類似，翼面渦環(huán)的前緣線布置在翼面網(wǎng)格的1/4弦線；渦環(huán)前緣線中點為氣動力作用點；控制點布置在翼面網(wǎng)格3/4弦線中點，在該點滿足無穿透邊界條件。

設(shè)翼面渦環(huán)數(shù)為N=ncnb(nc和nb分別為沿弦向和展向的網(wǎng)格數(shù))，尾跡渦環(huán)數(shù)為M=nwnb(nw為沿弦向的網(wǎng)格數(shù))。對于強迫振動，各個網(wǎng)格點的位移和速度為時間的函數(shù)，可以直接求解；對于氣動彈性問題，網(wǎng)格點的位移和速度通過結(jié)構(gòu)動力學模型求解得到。通過在各個翼面網(wǎng)格的控制點(Collocation Point)處滿足無穿透邊界條件建立非定常渦格法的控制方程，則各控制點處法向洗流速度vn與翼面渦強度γb和尾渦γw之間的關(guān)系為

(22)

即

AΓb=vn-BΓw

(23)

式中：A和B分別為翼面渦環(huán)和尾跡渦環(huán)的氣動力影響矩陣；Γb和Γw分別為翼面渦強度和尾跡渦強度；aij為第j個翼面渦環(huán)的單位渦強度在第i個翼面渦環(huán)控制點的誘導(dǎo)速度；bik為第k個尾跡渦環(huán)的單位渦強度在第i個翼面渦環(huán)控制點的誘導(dǎo)速度。在每一步時間推進求解中，根據(jù)翼面網(wǎng)格的變形量和速度更新A和vn。 尾跡渦環(huán)的渦強Γw是各個時刻機翼后緣渦脫落而產(chǎn)生的，且根據(jù)卡爾文環(huán)量定理，t+1時刻翼面鄰近的尾跡渦環(huán)的強度與t時刻翼面后緣處渦環(huán)的強度相等。因此，Γw可表示為

(24)

式中：Cb和Cw為常數(shù)矩陣。求解矩陣B的關(guān)鍵在于計算尾跡渦環(huán)的位置，該過程計算較為復(fù)雜，詳細過程可參考文獻[22]。

3 氣動彈性模型

基于機翼結(jié)構(gòu)非線性動力學模型和非定常氣動力模型，采用松耦合的方法建立大展弦比柔性機翼的非線性氣動彈性模型。耦合的難點在于氣動和結(jié)構(gòu)之間的插值，即結(jié)構(gòu)位移和速度插值到翼面渦環(huán)的節(jié)點，氣動力插值到結(jié)構(gòu)有限元節(jié)點上，但對于梁結(jié)構(gòu)模型，該過程很容易實現(xiàn)。圖6為氣動網(wǎng)格和結(jié)構(gòu)有限元離散示意圖。

在大展弦比柔性機翼的結(jié)構(gòu)建模中，假設(shè)機翼剖面是剛性的，即同一機翼剖面各點的相對位置保持不變。通過結(jié)構(gòu)動力學方程，可以求解彈性軸上各點的位置和速度，進而通過式(6)和式(7)可求解翼面任意點的位移和速度。同時，作用在翼面的氣動載荷可以直接通過虛功原理轉(zhuǎn)換至結(jié)構(gòu)有限元節(jié)點上。

本文采用Newmark-β方法對非線性動力學方程進行時域推進求解，非線性結(jié)構(gòu)動力學方程采用Newton-Raphson方法迭代求解，且為了提高求解效率，在動力學方程求解中，不對氣動力進行更新，具體求解過程見圖7。為了彌補這種松耦合求解方法帶來的精度損失，本文在時域響應(yīng)推進求解中采用較小的時間步長。

圖6 氣動網(wǎng)格和結(jié)構(gòu)有限元離散示意圖Fig.6 Schematic of aerodynamic panels and structure finite element discretization

圖7 非線性氣動彈性時域響應(yīng)流程Fig.7 Flow chart for nonlinear aeroelasticity time marching response

4 模型驗證與算例

Patil和Hodges采用的機翼[2]被廣泛應(yīng)用于考慮幾何非線性的非線性氣動彈性特性研究中，機翼的幾何和結(jié)構(gòu)參數(shù)如表1所示，平面形狀如圖8所示。機翼結(jié)構(gòu)模型和氣動模型沿展向采用相同的劃分方式，有限元個數(shù)為24個，氣動網(wǎng)格沿弦向布置8個，其中大氣密度為0.088 9 kg/m3。

表1 機翼幾何和結(jié)構(gòu)參數(shù)Table 1 Geometry and structure parameters of wing

圖8 懸臂機翼示意圖Fig.8 Schematics of cantilever wing

4.1 大展弦比柔性機翼顫振速度求解

本文采用時域仿真的方法求解顫振速度。在給定迎角下，逐漸增加來流速度，在靜氣彈平衡位置的基礎(chǔ)上給機翼一個小擾動(0.02°迎角)，進行時域仿真，通過觀察機翼是否呈等幅振蕩判斷其是否達到臨界顫振速度。

圖9為翼根迎角為3°時，不同來流速度U∞下機翼翼梢扭轉(zhuǎn)角度的時域響應(yīng)對比圖，從圖中可以看出，來流速度U∞=23.7 m/s時，隨著時間的推進振蕩振幅減小，逐漸收斂；來流速度U∞=24.8 m/s時，隨著時間的推進振幅增加，逐漸發(fā)散；來流速度U∞=24.2 m/s時，為等幅振蕩。因此顫振速度為24.2 m/s(在求解時，進行了不同速度下大量的時域仿真計算，以保證該速度為呈現(xiàn)等幅振蕩的最小速度)。

圖9 3°翼根迎角時不同來流速度翼梢扭轉(zhuǎn)的時域響應(yīng)Fig. 9 Time history of tip twist at various free-stream speeds for root angle of attack=3°

圖10 不同翼根迎角下的顫振速度Fig.10 Flutter speed at various root angles of attack

圖10為根據(jù)上述準則求解的機翼顫振速度隨機翼迎角的變化曲線，從圖中可以看出，本文計算結(jié)果與文獻[17]基本相同，與文獻[8]趨勢一致，但有一個約2 m/s的差量，差量主要是因為文獻采用的氣動力模型是基于片條理論的二元機翼有限狀態(tài)氣動力模型，無法考慮機翼的三維效應(yīng)，計算的氣動力偏大，使得其計算的顫振速度偏小。

4.2 來流迎角和機翼展長對顫振速度的影響

展長是決定機翼柔性大小的關(guān)鍵參數(shù)，主導(dǎo)著飛行過程中機翼產(chǎn)生大變形的程度。從4.1節(jié)可以看出，來流迎角對顫振速度也有明顯的影響。本節(jié)將進一步探索這兩個參數(shù)對顫振速度的影響。其余各個參數(shù)與表1相同。

圖11為不同機翼展長在不同翼根迎角下的顫振速度變化圖。從圖中可以看出，對于不同展長的機翼，隨著迎角的增加，顫振速度逐漸減小，但減緩的速率越來越小，來流迎角從1°增加至3°時，顫振速度降低約20%，可見大展弦比柔性機翼的顫振速度對來流迎角較為敏感。

圖11 不同翼根迎角和展長機翼的顫振速度Fig.11 Flutter speed at various root angles of attack with different span-wise lengths

隨著機翼展長的增加，機翼的顫振速度也有明顯的變化，不同迎角下，顫振速度均逐漸減小，不過減緩速率也逐漸減小，機翼半展長由12 m增加至18 m時，顫振速度降低約30%。從氣動的角度來說，增加展長可以減小誘導(dǎo)阻力，進而提高升阻比，但機翼的氣動彈性特性惡化，在設(shè)計時應(yīng)綜合考慮。

4.3 非線性氣動彈性時域響應(yīng)特性

對于線性氣動彈性模型，當來流速度大于顫振速度時，機翼運動的振幅隨著時間呈指數(shù)形式增加；當考慮到非線性后，將呈現(xiàn)更復(fù)雜的響應(yīng)。本節(jié)將研究在初始小擾動條件下(0.02°來流迎角)，來流速度大于顫振速度時，大展弦比柔性機翼的時域響應(yīng)特性。

圖12為3°迎角下，來流速度大于機翼顫振速度24.2 m/s時，翼梢扭轉(zhuǎn)角位移的時域響應(yīng)。從圖12可以看到，在初始階段，響應(yīng)振幅逐漸增加，而后機翼運動收斂至穩(wěn)定的極限環(huán)振蕩。關(guān)于機翼運動沒有持續(xù)發(fā)散的原因，文獻[6-7,13-16]將其歸因于結(jié)構(gòu)幾何非線性或者氣動非線性，但并未進一步深入分析。本文接下來對極限環(huán)產(chǎn)生原因進行簡要敘述。

圖12 3°翼根迎角不同來流速度下翼梢扭轉(zhuǎn)角的時域響應(yīng)Fig.12 Time history of tip twist at various free-stream speeds for 3° root angle of attack

另一方面，當變形較大時，機翼展向長度減小，且剖面方向也發(fā)生變化，使得相同迎角下機翼的升力和力矩也減小。以翼根迎角為3°、來流速度為25.7 m/s時，75%展向位置處剖面俯仰力矩非線性為例進行分析。圖13給出了在20 s和60 s 附近一個周期內(nèi)剖面俯仰力矩隨扭轉(zhuǎn)角的變化曲線。在20 s時，剖面扭轉(zhuǎn)角度振動范圍較小，大變形效應(yīng)不明顯，力矩線性特性較好；在60 s 時，扭轉(zhuǎn)角振動范圍較大，由于大變形效應(yīng)，在0°扭轉(zhuǎn)角增加至最大扭轉(zhuǎn)角過程中，有明顯的力矩損失，該非線性效應(yīng)減小了扭轉(zhuǎn)角增加的趨勢，有利于維持機翼的等幅振蕩。

圖13 3°翼根迎角、U∞=25.7 m/s時，75%展向位置 處剖面扭轉(zhuǎn)力矩隨扭轉(zhuǎn)角的變化曲線Fig.13 Section moments of torque curves for U∞=25.7 m/s in 75% spanwise position with respect to twist angle for 3° root angle of attack

圖14給出了3°來流迎角不同來流速度下，機翼翼梢扭轉(zhuǎn)角的幅值變化。隨著來流速度的增加，機翼在平衡位置附近小幅振蕩發(fā)散的速率先增加后減?。辉诰€性氣彈系統(tǒng)中，振動發(fā)散的速率和系統(tǒng)的阻尼正相關(guān)，系統(tǒng)的阻尼越大，發(fā)散的速率越大；而機翼在平衡位置附近小幅振蕩范圍內(nèi)，機翼幾何非線性效應(yīng)較弱，因此可以認為氣動彈性系統(tǒng)在平衡位置處的阻尼也先增加后減小，這和文獻[9]采用幾何非線性梁和二元非定常氣動力模型計算結(jié)果的趨勢相同。從圖12(g)可以看出，當來流速度為27.48 m/s時，機翼在平衡位置附近維持小振幅振蕩，而非先發(fā)散后維持極限環(huán)振蕩，可見此時系統(tǒng)處于臨界穩(wěn)定狀態(tài)，系統(tǒng)由不穩(wěn)定變?yōu)榉€(wěn)定。隨著速度進一步增加，系統(tǒng)的振蕩收斂(圖12(h))。

圖14 3°翼根迎角不同來流速度下翼梢扭轉(zhuǎn)角的幅值 Fig.14 Amplitude of tip twist at various free-stream speeds for 3° root angle of attack

綜上可以看出，來流速度增加至24.2 m/s時，機翼由收斂運動變?yōu)殪o平衡位置附近的小幅振蕩運動，系統(tǒng)的阻尼由負變?yōu)?，處于臨界穩(wěn)定狀態(tài)； 隨著來流速度進一步增加，系統(tǒng)阻尼為正，機翼在靜平衡位置附近呈現(xiàn)發(fā)散振蕩，且發(fā)散率先增加后減??；直到來流速度為27.48 m/s時，系統(tǒng)發(fā)散率減小至0，機翼在靜平衡位置附近小幅振蕩，系統(tǒng)再次處于臨界穩(wěn)定狀態(tài)；隨著來流速度進一步增加，系統(tǒng)表現(xiàn)為收斂振蕩運動，直到速度增加至39.4 m/s時，系統(tǒng)第3次達到臨界穩(wěn)定狀態(tài)(見圖15)。

圖15 3°翼根迎角不同來流速度下翼梢扭轉(zhuǎn)角的 時域響應(yīng) (t=0～20 s)Fig.15 Time history of tip twist at various free-stream speeds for 3° root angle of attack (t=-20 s)

5 結(jié) 論

基于牛頓定律推導(dǎo)了幾何非線性梁的動力學方程，與文獻[5]的方法相比，推導(dǎo)過程簡單、直觀，物理意義明確，便于理解動力學方程，這也是在建立動力學方程中，牛頓力學方法相比于分析力學的優(yōu)勢。該方法可以與分析力學方法相互補充，相互驗證。

基于機翼結(jié)構(gòu)非線性動力學模型和非定常渦格法，采用松耦合的方法建立了大展弦比機翼的非線性氣動彈性模型，該模型考慮了機翼大變形、非定常氣動力的三維效應(yīng)以及氣動面隨結(jié)構(gòu)變形的自適應(yīng)更新。

大展弦比柔性機翼的顫振特性對來流迎角較為敏感，機翼迎角由1°增加至3°時，不同展長的柔性機翼的顫振速度均降低了約20%。同時，隨著機翼展長的增加，機翼的顫振速度也有明顯的變化，各個迎角下，機翼半展長由12 m增至18 m時，顫振速度降低約30%。

由于結(jié)構(gòu)幾何非線性效應(yīng)，來流速度超過顫振速度時，大展弦比柔性機翼的振動并不發(fā)散，而是呈現(xiàn)出穩(wěn)定的極限環(huán)振蕩。隨著來流速度進一步增加，系統(tǒng)再次穿過臨界平衡點(U∞=27.48 m/s)，不穩(wěn)定系統(tǒng)變?yōu)榉€(wěn)定系統(tǒng)，直到隨著速度的增加系統(tǒng)第3次達到臨界穩(wěn)定狀態(tài)(U∞=39.4 m/s)。

[1] NOLL T E, ISHMAEL S D, HENWOOD B, et al. Technical findings, lessons learned, and recommendations resulting from the Helios prototype vehicle mishap[C]∥Processes and Criteria for UAV Structural Integrity, 2007.

[2] PATIL M J, HODGES D H. On the importance of aerodynamic and structural geometrical nonlinearities in aeroelastic behavior of high-aspect-ratio wings[J]. Journal of Fluids and Structures, 2004, 19(7): 905-915.

[3] TANG D M, DOWELL E H. Effects of geometric structural nonlinearity on flutter and limit cycle oscillations of high-aspect-ratio wings[J]. Journal of Fluids and Structures, 2004, 19(3): 291-306.

[4] AFONSO F，VALE J， OLIVEIRA é, et al. A review on non-linear aeroelasticity of high aspect-ratio wings[J]. Progress in Aerospace Sciences, 2017, 89: 40-57.

[5] HODGES D H. A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams[J]. International Journal of Solids and Structures, 1990, 26(11): 1253-1273.

[6] PATIL M J, HODGES D H, CESNIK C E S. Limit-cycle oscillations in high-aspect-ratio wings[J]. Journal of Fluids and Structures, 2001, 15(1):107-132.

[7] XIANG J W, YAN Y J, LI D C. Recent advance in nonlinear aeroelastic analysis and control of the aircraft[J]. Chinese Journal of Aeronautics, 2014, 27(1): 12-22.

[8] PATIL M J, HODGES D H, CESNIK C E S. Nonlinear aeroelasticity and flight dynamics of high-altitude long-endurance aircraft[J]. Journal of Aircraft, 2001, 38(1): 88-94.

[9] 段靜波, 周洲, 王偉, 等. 大展弦比大柔性機翼載荷分布求解的一種方法[J]. 航空學報, 2016, 37(3): 799-809.

DUAN J B, ZHOU Z, WANG W, et al. A method for aeroelastic load redistribution of very flexible wing with a high-aspect-ratio[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37 (3): 799-809 (in Chinese).

[10] 王偉, 段卓毅, 耿建中, 等. 基于CR 理論的大柔性機翼幾何非線性結(jié)構(gòu)建模[J]. 航空學報,2017, 38(S1):721544.

WANG W, DUAN Z Y, GENG J Z, et al. Geometrically nonlinear structural model for very flexible wing based on CR theory[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(S1): 721544(in Chinese).

[11] 楊超, 王立波, 謝長川, 等. 大變形飛機配平與飛行載荷分析方法[J]. 中國科學: 技術(shù)科學, 2012, 42(10): 1137-1147.

YANG C, WANG L B, XIE C C, et al. Aeroelastic trim and flight loads analysis of flexible aircraft with large deformations[J]. Sci China Tech Sci, 2012, 42 (10): 1137-1147 (in Chinese).

[12] WANG W, ZHU X P, ZHOU Z, et al. A method for nonlinear aeroelasticity trim and stability analysis of very flexible aircraft based on co-rotational theory[J]. Journal of Fluids and Structures, 2016, 62: 209-229.

[13] TANG D M, DOWELL E H. Experimental and theoretical study on aeroelastic response of high-aspect-ratio wings[J]. AIAA Journal, 2001, 39(8): 1430-1441.

[14] ESKANDARY K, DARDEL M, PASHAEI M H, et al. Nonlinear aeroelastic analysis of high-aspect-ratio wings in low subsonic flow[J]. Acta Astronautica, 2012, 70(1): 6-22.

[15] ARENA A, LACARBONARA W, MARZOCCA P. Nonlinear aeroelastic formulation and postflutter analysis of flexible high-aspect-ratio wings[J]. Journal of Aircraft, 2013, 50(6): 1748-1764.

[16] SHAMS S, SADR M H, HADDADPOUR H. An efficient method for nonlinear aeroelasticy of slender wings[J]. Nonlinear Dynamics, 2012, 67(1): 659-681.

[17] 楊智春, 張惠, 谷迎松, 等. 考慮幾何非線性效應(yīng)的大展弦比機翼氣動彈性分析[J]. 振動與沖擊, 2014,33(16): 72-75.

YANG Z C, ZHANG H, GU Y S, et al. Aeroelastic analysis of the high aspect ratio wing considering the geometric nonlinearity[J]. Journal of Vibration and Shock, 2014, 33(16): 72-75 (in Chinese).

[18] SULEMAN A, AFONSO F, VALE J, et al. Non-linear aeroelastic analysis in the time domain of high-aspect-ratio wings: Effect of chord and taper-ratio variation[J]. Aeronautical Journal, 2017, 121(1235): 21-52.

[19] 謝長川, 胡銳, 王斐, 等. 大展弦比柔性機翼氣動彈性風洞模型設(shè)計與試驗驗證[J]. 工程力學, 2016, 33(11): 249-256.

XIE C C, HU R, WANG F, et al. Aeroelastic wind tunnel test model design and experiment on very flexible high-aspect-ratio wings[J]. Engineering Mechanics, 2016, 33 (11): 249-256 (in Chinese).

[20] LIU Y, XIE C C, YANG C, et al. Gust response analysis and wind tunnel test for a high-aspect ratio wing[J]. Chinese Journal of Aeronautics, 2016, 29(1): 91-103.

[21] KATZ J, PLOTKIN A. Low-speed aerodynamics[M]. Cambridge: Cambridge University Press, 2001.

[22] HODGES D H. Nonlinear composite beam theory[M]. Reston, VA: AIAA, 2006.

[23] GéRADIN M, CARDONA A. Flexible multibody dynamics: A finite element approach[M]. Chichester: Wiley, 2001.

Aeroelasticmodelingandanalysisofwingsconsideringgeometricnonlinearity

GUOTongbiao1,BAIJunqiang1, *,SUNZhiwei2,WANGChen1

1.SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China2.UAVResearchInstitute,NorthwesternPolytechnicalUniversity,Xi’an710065,China

Flexiblewingswithhighaspectratioarewidelyusedinhigh-altitudeandlong-enduranceUnmannedAerialVehicles(UAVs)becauseoflowstructuralweightandhighaerodynamiclift-to-dragratio.Thiskindofwingsexperiencelargegeometricdeformationinflight,andthelinearstructuremodelbasedonsmalldeformationhypothesisisnolongerapplicable.Therefore,itisnecessarytobuildthestructuremodelwhichcansimulategeometricnonlinearity.BasedontheNewtonianmethod,thedynamicequationsforthegeometricnon-linearstructuremodelarederived,whichcanbemutuallyvalidatedbyandcomplementedwiththemethodbasedonHamilton’sprinciplederivedbyHodges.Tosimulatetheaerodynamicsofflexiblewingsmoreprecisely,amodelforthree-dimensionalunsteadyaerodynamics,whichcanconsiderlargedeformationofthewing,isbuilt.Basedonthenonlinearstructuremodelandtheunsteadyaerodynamicmodel,thenonlinearaeroelasticmodelisbuiltthroughloosecoupling.Theprecisionoftheaeroelasticmodelisverifiedthroughtests.Theresultsshowthattheflutterspeedofflexiblewingsissensitivetothefree-streamanglesofattackandspan-wiselength.Whenthefree-streamspeedsexceedtheflutterspeed,thewing’svibrationsarestableLimitedCycleOscillations(LCO)ratherthandivergence.However,asthefree-streamspeedscontinuesincreasing,thewing’svibrationsconvergeagainandthedampingturnedtobepositive.

nonlinearaeroelasticity;geometricnonlinearity;flexiblewing;flutterspeed;limitedcycleoscillations;timemarchingresponse

2017-04-25；Revised2017-05-15；Accepted2017-07-07；Publishedonline2017-07-121459

URL：http://hkxb.buaa.edu.cn/CN/html/20171112.html

Industry,EducationandResearchFund(cxy2014XGD09)

.E-mailjunqiang@nwpu.edu.cn

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2017.121351

V211.47

A

1000-6893(2017)11-121351-11

2017-04-25；退修日期2017-05-15；錄用日期2017-07-07；< class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2017-07-121459

http://hkxb.buaa.edu.cn/CN/html/20171112.html

產(chǎn)學研基金(cxy2014XGD09)

.E-mailjunqiang@nwpu.edu.cn

郭同彪，白俊強，孫智偉，等．考慮機翼幾何非線性的氣動彈性建模與分析J． 航空學報,2017,38(11):121351.GUOTB,BAIJQ,SUNZW,etal.AeroelasticmodelingandanalysisofwingsconsideringgeometricnonlinearityJ.ActaAeronauticaetAstronauticaSinica,2017,38(11):121351.

(責任編輯：李明敏)