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

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

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

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

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

中圖分類號(hào):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.

基金項(xiàng)目：國(guó)家自然科學(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

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

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

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

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

大柔性機(jī)翼在氣動(dòng)載荷的作用下，變形前后的幾何位形具有明顯的差異，常規(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)建立，通過(guò)坐標(biāo)轉(zhuǎn)換矩陣及其微分形式構(gòu)造總體坐標(biāo)系下的切線剛度矩陣，因而其核心是坐標(biāo)轉(zhuǎn)換矩陣的建立及其微分形式的推導(dǎo)，其計(jì)算精度主要取決于剛度矩陣和坐標(biāo)系轉(zhuǎn)換矩陣對(duì)幾何非線性描述的精確性。

1.1CR梁?jiǎn)卧o力學(xué)平衡方程

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

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

ui(x)=Ti(x)Ixyz

(1)

那么參照文獻(xiàn)[9,14]，CR空間梁?jiǎn)卧膯卧鴺?biāo)系定義為：

ue1=(d2-d1)/ln

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

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

(2)

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

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

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

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

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)

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

δpl=Tδp

(4)

總體坐標(biāo)系下廣義節(jié)點(diǎn)力向量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é)點(diǎn)力fl的函數(shù)。

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

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

(8)

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

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

(9)

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

1.2CR梁?jiǎn)卧獎(jiǎng)恿W(xué)平衡方程

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

(10)

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

(11)

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

(12)

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

(13)

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

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

(14)

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

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

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

Fe-fmass-fi=0

(16)

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

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

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

(17)

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

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

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

(18)

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

(19)

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

(20)

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

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

(21)

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

2非線性非定常氣動(dòng)力模型

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

(22)

(23)

(24)

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

(25)

(26)

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

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

(27)

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

第二節(jié)中的動(dòng)力學(xué)方程求解是以總體坐標(biāo)系下的節(jié)點(diǎn)線位移、角位移、線速度、線加速度和單元坐標(biāo)系下角速度、角加速度為未知量的；非線性非定常氣動(dòng)力的求解輸入是局部氣流坐標(biāo)系下的片條沉浮線位移、線速度、線加速度、角位移、角速度、角加速度等，因此需要對(duì)動(dòng)力學(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)過(guò)插值得到每個(gè)片條的振動(dòng)及沉浮、扭轉(zhuǎn)運(yùn)動(dòng)：

(30)

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

(31)

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

(32)

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

Qe-fmass-fi=0

(33)

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

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

(34)

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

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

4算例及驗(yàn)證

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

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

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

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

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

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

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

4.2非線性非定常氣動(dòng)力模型驗(yàn)證

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.5積分步長(zhǎng)的選擇及收斂性

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

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

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

5結(jié)論

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

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

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

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

(4)前移彈性軸可以有效地提高臨界極限環(huán)顫振速度，增加周期性極限環(huán)顫振的速度區(qū)間，顯著地改善大柔性機(jī)翼的非線性氣動(dòng)彈性特性，可以用于指導(dǎo)大柔性機(jī)翼彈性軸在弦向的站位設(shè)計(jì)。進(jìn)一步的研究工作中可以引入高效的預(yù)估求解算子，減少動(dòng)力學(xué)響應(yīng)求解時(shí)的迭代次數(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列式三維梁?jiǎn)卧?jì)算方法[J]. 西南交通大學(xué)學(xué)報(bào)，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)的大柔性太陽(yáng)能無(wú)人機(jī)靜氣動(dòng)彈性分析[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 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]吳國(guó)榮,顏桂云，陳福全.柔性梁大柔度動(dòng)力響應(yīng)分析的多體系統(tǒng)方法[J].振動(dòng)與沖擊,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è)向隨動(dòng)力作用下大展弦比柔性機(jī)翼的穩(wěn)定性[J].航空學(xué)報(bào),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.

第一作者何洪陽(yáng)男，碩士，1988年生

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