王偉, 周洲, 祝小平, 王睿

太陽能無人機(jī)能夠完成大范圍、長(zhǎng)時(shí)間的大氣環(huán)境探測(cè)、軍事偵查和通信中繼任務(wù),并在各國(guó)頻頻傳出新的研究進(jìn)展[1-5]。這類無人機(jī)一般具有較大的展弦比、較低的結(jié)構(gòu)面密度等特點(diǎn),比如“太陽神”太陽能無人機(jī)的展弦比達(dá)到31,而常規(guī)動(dòng)力的高空長(zhǎng)航時(shí)無人機(jī)“全球鷹”的展弦比為25,并且兩者的機(jī)翼結(jié)構(gòu)面密度分別為3.2 kg/m2和53kg/m2。太陽能無人機(jī)的大展弦比、輕質(zhì)結(jié)構(gòu)設(shè)計(jì)特點(diǎn)使其在飛行中產(chǎn)生明顯的機(jī)翼結(jié)構(gòu)靜變形,如“探路者”太陽能無人機(jī)平飛時(shí)機(jī)翼上翹量超過半展長(zhǎng)的12%,“太陽神”在極限飛行狀態(tài)下的機(jī)翼上反角可達(dá)到50°[3-4]。2003年6月,加裝實(shí)驗(yàn)設(shè)備的Helios無人機(jī)飛行時(shí),在低空出現(xiàn)振蕩導(dǎo)致全機(jī)解體墜毀。2004年1月,NASA發(fā)布的事故調(diào)查報(bào)告中提出繼續(xù)研制這種無人機(jī),并著重強(qiáng)調(diào)了對(duì)這種無人機(jī)非線性氣動(dòng)彈性特性分析方法的研究[6]。扭轉(zhuǎn)變形較小時(shí),采用線彈性結(jié)構(gòu)模型模型研究靜氣動(dòng)彈性問題,氣動(dòng)力隨撓度的增大而單調(diào)增加,變形略大時(shí)與實(shí)際嚴(yán)重不符[7]。因此,這類大變形無人機(jī)的結(jié)構(gòu)建模必須考慮大位移、大轉(zhuǎn)角等幾何非線性特征。

目前,大柔性機(jī)翼的幾何非線性結(jié)構(gòu)模型主要有Hodges-Dowell理論梁模型和基于Total Lagrange(TL法)、Updated Lagrange(UL法)和共旋坐標(biāo)法(CR)的有限元模型。其中,在NASA的資助下,基于Hodges-Dowell非線性梁模型的高空長(zhǎng)航時(shí)(HALE)飛機(jī)非線性氣動(dòng)彈性配平及穩(wěn)定性計(jì)算代碼(NATASHA)的開發(fā)得到了一定的發(fā)展,可以預(yù)測(cè)飛翼布局HALE飛機(jī)的氣動(dòng)彈性響應(yīng)特性,但變形的計(jì)算需要預(yù)先求解結(jié)構(gòu)內(nèi)力[7-10];傳統(tǒng)的TL法變形前后單元坐標(biāo)系保持不變,無法考慮變形后結(jié)構(gòu)狀態(tài)的影響;UL法中單元兩端坐標(biāo)轉(zhuǎn)換矩陣不一致,也無法考慮彎曲變形的影響;因此,這2種方法通常需要增加載荷步以提高求解精度,從而顯著增加計(jì)算時(shí)間。CR有限元法是近年來發(fā)展起來的新的幾何非線性計(jì)算方法,主要思想是:把結(jié)構(gòu)的運(yùn)動(dòng)分解為純粹的彈性變形和剛體的平移及轉(zhuǎn)動(dòng),位移-應(yīng)變關(guān)系建立在局部坐標(biāo)系內(nèi),不像傳統(tǒng)的TL法那樣從非線性位移-應(yīng)變關(guān)系出發(fā)來推導(dǎo)切線剛度矩陣,并且具有更高的求解精度和計(jì)算效率[11-12]。因此,本文針對(duì)大柔性太陽能無人機(jī)機(jī)翼結(jié)構(gòu)的大位移、大轉(zhuǎn)動(dòng)、小應(yīng)變特征,通過建立節(jié)點(diǎn)隨動(dòng)坐標(biāo)系和單元隨動(dòng)坐標(biāo)系,推導(dǎo)了基于CR有限元理論的切線剛度矩陣、內(nèi)力求解格式及增量平衡方程,編寫了基于CR有限元理論的幾何非線性結(jié)構(gòu)靜力學(xué)求解代碼;氣動(dòng)載荷在3個(gè)方向上進(jìn)行積分,以描述結(jié)構(gòu)大變形引起的氣動(dòng)載荷作用方向的改變。

1 大柔性結(jié)構(gòu)幾何非線性靜力學(xué)建模

大柔性無人機(jī)機(jī)翼可以用柔性梁模型來描述[7-9],因此本文選擇空間CR梁?jiǎn)卧⒋笕嵝蕴柲軣o人機(jī)的結(jié)構(gòu)模型。采用有限元技術(shù)對(duì)大柔性結(jié)構(gòu)建模時(shí),其計(jì)算精度主要取決于剛度矩陣和坐標(biāo)系轉(zhuǎn)換矩陣對(duì)幾何非線性描述的精確性。CR法的平衡方程在單元坐標(biāo)系內(nèi)建立,通過坐標(biāo)轉(zhuǎn)換矩陣及其微分形式構(gòu)造總體坐標(biāo)系下的切線剛度矩陣,因而其核心是坐標(biāo)轉(zhuǎn)換矩陣的建立及其微分形式的推導(dǎo)。

1.1 坐標(biāo)系的建立

坐標(biāo)轉(zhuǎn)換矩陣是單元坐標(biāo)系與總體坐標(biāo)系之間的聯(lián)系紐帶,為了較精確地描述坐標(biāo)系轉(zhuǎn)換矩陣,本文引入總體坐標(biāo)系Ixyz、節(jié)點(diǎn)坐標(biāo)系ui、單元平均構(gòu)形坐標(biāo)系um和單元坐標(biāo)系ue。

圖1 柔性結(jié)構(gòu)變形示意圖

當(dāng)結(jié)構(gòu)位形確定時(shí),單元節(jié)點(diǎn)的3個(gè)歐拉角及空間坐標(biāo)隨之確定,按照式1旋轉(zhuǎn)總體坐標(biāo)系可以得到節(jié)點(diǎn)坐標(biāo)系。

ui(x)=Ti(x)Ixyz

(1)

Ti(x)是節(jié)點(diǎn)i的坐標(biāo)系轉(zhuǎn)換矩陣,是3個(gè)歐拉角φx,φy,φz的函數(shù),其中:

對(duì)于2節(jié)點(diǎn)空間梁?jiǎn)卧?采用(1)式可以得到節(jié)點(diǎn)坐標(biāo)系ui和uj,假設(shè)其坐標(biāo)轉(zhuǎn)換矩陣Ti和Tj對(duì)應(yīng)的旋轉(zhuǎn)偽矢量為(θix,θiy,θiz)T和(θjx,θjy,θjz)T,參照文獻(xiàn)[11],以γ/2作為旋轉(zhuǎn)偽矢量,對(duì)節(jié)點(diǎn)坐標(biāo)系ui進(jìn)行旋轉(zhuǎn)得到平均構(gòu)形坐標(biāo)系um。

um=R(γ/2)ui

(2)

式中：γ=(θjx-θix,θjy-θiy,θjz-θiz)T

S(γ/2)S(γ/2)λ=((γ/2)(γ/2)T)1/2

I為單位陣。

S(γ/2)=

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

(3)

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

1.2 增量平衡方程推導(dǎo)

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

(4)

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

δpl=Tδp

(5)

單元坐標(biāo)系下的廣義節(jié)點(diǎn)力fl為:

fl=Klpl

(6)

式中:Kl是線性剛度矩陣,與線性梁?jiǎn)卧膭偠染仃囈恢隆?/p>

則總體坐標(biāo)系下廣義節(jié)點(diǎn)力向量f為:

(7)

(8)

TTKlT為材料剛度矩陣,KTσ(fl)為幾何剛度矩陣,是單元坐標(biāo)系下單元節(jié)點(diǎn)力fl的函數(shù)。大柔性結(jié)構(gòu)幾何非線性問題求解時(shí),0,Δt,2Δt,…,t時(shí)刻的結(jié)構(gòu)狀態(tài)是已知的,t+Δt時(shí)刻的狀態(tài)未知量是待求解的。那么,結(jié)構(gòu)靜力學(xué)增量平衡方程為:

KTΔp=Ft+Δt-ft

(9)

2 幾何非線性靜氣動(dòng)彈性建模

假設(shè)機(jī)翼剖面為剛性,變形后的機(jī)翼外形通過一系列控制剖面的平移和轉(zhuǎn)動(dòng)經(jīng)過擬合得到;如圖2所示,剖面i、j分別由結(jié)構(gòu)節(jié)點(diǎn)i、j的線位移和角位移控制。流場(chǎng)控制方程采用三維可壓N-S定常方程,湍流模型采用k-kl-ω轉(zhuǎn)捩模型,流場(chǎng)計(jì)算采用Fluent軟件。氣動(dòng)力在總體坐標(biāo)系的3個(gè)方向分別進(jìn)行積分,然后根據(jù)功互等原理,將分布的氣動(dòng)力載荷轉(zhuǎn)換到節(jié)點(diǎn)上,得到等效節(jié)點(diǎn)力,計(jì)算公式為:

(10)

圖2 氣動(dòng)載荷插值示意圖

那么,靜氣動(dòng)彈性運(yùn)動(dòng)方程為:

KTΔp=Qt-ft

(11)

收斂控制方程為:

(12)

一般,前2～3次主迭代時(shí)Qt-ft較大,需要對(duì)載荷進(jìn)行分步施加,以保證求解精度及收斂效率。

3 代碼驗(yàn)證及算例分析

3.1 懸臂梁在集中彎矩下的位移求解

根據(jù)上述計(jì)算方法編制了非線性有限元求解代碼,為了驗(yàn)證其計(jì)算精度和效率,選擇具有解析解的大柔性懸臂梁作為算例。懸臂梁自由端受集中彎矩M作用,取l=20 m,EIz=10 Nm,EA=10 N,M分別取0.2π、0.4π、0.6π、0.8π和π。其解析解為弧長(zhǎng)20 m,圓心角θ=0.4π、0.8π、1.2π、1.6π和2π的圓弧。將懸臂梁劃分為20個(gè)單元,采用所述載荷增量法進(jìn)行求解,收斂控制方程采用(12)式,子迭代采用Newton-Raphson迭代格式,位移增量允許誤差<10-3,變形如圖3所示。

圖3 懸臂梁幾何大變形

上述算例求解中,子迭代經(jīng)過3～4次迭代收斂,且最大誤差不超過3%,可見CR有限元法具有較好的求解精度及求解效率,并且能夠非常精確地模擬到結(jié)構(gòu)的幾何非線性變形特征。

3.2 大柔性太陽能無人機(jī)非線性靜氣動(dòng)彈性計(jì)算與結(jié)果分析

大柔性太陽能無人機(jī)的幾何參數(shù)、結(jié)構(gòu)參數(shù)及飛行條件如圖4及表1所示:

圖4 太陽能無人機(jī)幾何模型

表1 太陽能無人機(jī)模型參數(shù)

以來流速度作為變化參數(shù),氣動(dòng)載荷隨來流速度的增加而增加,以描述結(jié)構(gòu)在不同氣動(dòng)載荷作用下的彈性變形,隨后凍結(jié)氣動(dòng)外形,研究該幾何大變形對(duì)氣動(dòng)特性及橫航向穩(wěn)定性的影響。不同氣動(dòng)載荷作用下的機(jī)翼變形歷程如圖5所示。

圖5 不同來流速度下的機(jī)翼變形歷程

來流速度V=0 m/s表示無人機(jī)初始剛體外形,V=15 m/s 、20 m/s、25 m/s、30 m/s、35 m/s時(shí),所引起的翼尖撓度分別為展長(zhǎng)的5%、8.9%、13.3%、18%、22.5%,引起的翼尖展向位移分別為展長(zhǎng)的0.5%、1.3%、2.7%、4.9%,7.5%,相當(dāng)于展長(zhǎng)縮短了1%、2.6%、5.4%、9.8%,15%;另外,可以看出氣動(dòng)載荷作用的方向也發(fā)生了明顯的變化。

機(jī)翼的彈性變形會(huì)引起載荷的重新分布。無人機(jī)在不同來流速度作用下,結(jié)構(gòu)發(fā)生幾何大變形時(shí)的升力重新分布情況如圖6所示。當(dāng)?shù)厣ο禂?shù)計(jì)算時(shí)選擇局部機(jī)翼面積作為參考面積,默認(rèn)無人機(jī)迎角為0°。

機(jī)翼的彈性變形使氣動(dòng)載荷向翼根轉(zhuǎn)移,這種幾何非線性效應(yīng)改善了機(jī)翼結(jié)構(gòu)的外載荷分布,對(duì)結(jié)構(gòu)設(shè)計(jì)是有利的,而線彈性結(jié)構(gòu)模型無法精確地模擬到這種趨勢(shì),甚至得到相反的結(jié)論。

圖6 升力沿展向分布 圖7 無人機(jī)極曲線與升組比曲線

當(dāng)翼尖撓度為22.5%的展長(zhǎng)時(shí),升力系數(shù)降低了12.4%,升阻比降低了11.8%。隨著變形增加,阻力系數(shù)基本不變,而升力系數(shù)減小比較快,因而升阻比的變化趨勢(shì)也是減小的。這主要是因?yàn)閹缀未笞冃胃淖兞藲鈩?dòng)載荷的作用方向,對(duì)全機(jī)升力的影響較大。采用線彈性結(jié)構(gòu)模型時(shí),氣動(dòng)力是隨彎曲變形的增加而單調(diào)增加的[7],因此采用線性結(jié)構(gòu)模型與非線性結(jié)構(gòu)模型研究較大彈性變形對(duì)氣動(dòng)力的影響時(shí),會(huì)得到不同的結(jié)論。

太陽能無人機(jī)機(jī)翼的幾何大變形相當(dāng)于增加了機(jī)翼的上反角,嚴(yán)重影響了無人機(jī)的橫航向力矩特性,如圖8所示,Cl為滾轉(zhuǎn)力矩系數(shù)、Cn為偏航力矩系數(shù),Clβ和Cnβ分別是滾轉(zhuǎn)力矩導(dǎo)數(shù)和偏航力矩導(dǎo)數(shù)。

圖8 橫航向力矩系數(shù)

隨著彈性變形的增加,Clβ和Cnβ的絕對(duì)值都呈增加趨勢(shì);從圖8中還可以看出,幾何非線性效應(yīng)對(duì)無人機(jī)橫航向穩(wěn)定性的影響較大。

當(dāng)彎曲變形為22.5%的展長(zhǎng)時(shí),滾轉(zhuǎn)力矩導(dǎo)數(shù)為剛性飛機(jī)的6.65倍,偏航力矩導(dǎo)數(shù)為剛性飛機(jī)的10.3倍,Cnβ/Clβ為剛性飛機(jī)1.66倍。因此,這類飛機(jī)的大彈性變形使?jié)L轉(zhuǎn)力矩和偏航力矩得到很大提高,且顯著改善了橫航向靜穩(wěn)定性。

表2 橫航向氣動(dòng)特性

4 結(jié) 論

本文結(jié)合CR有限元法及計(jì)算流體力學(xué)理論,對(duì)大柔性太陽能無人機(jī)的靜氣動(dòng)彈性問題進(jìn)行了深入的研究,得到以下結(jié)論:

1)大柔性太陽能無人機(jī)的結(jié)構(gòu)建模應(yīng)采用能夠考慮幾何非線性效應(yīng)的結(jié)構(gòu)模型。CR有限元法具有較好的計(jì)算精度和求解效率,能夠較好地模擬大變形無人機(jī)的幾何非線性特征,可以用來建立大柔性太陽能無人機(jī)較高精度的非線性結(jié)構(gòu)模型。

2)幾何非線性效應(yīng)改善了氣動(dòng)載荷的重新分布,有利于機(jī)翼結(jié)構(gòu)設(shè)計(jì);但同時(shí)會(huì)降低無人機(jī)的氣動(dòng)性能。

3)對(duì)具有正常橫航向穩(wěn)定性的飛機(jī),當(dāng)Cnβ/Clβ過大時(shí),飛機(jī)易產(chǎn)生螺旋不穩(wěn)定;當(dāng)Cnβ/Clβ過小時(shí),則飛機(jī)易產(chǎn)生荷蘭滾或飄擺不定,而彈性變形對(duì)Cnβ及Clβ的影響比較明顯,大柔性太陽能無人機(jī)在總體設(shè)計(jì)時(shí),應(yīng)考慮幾何大變形對(duì)氣動(dòng)特性及橫航向穩(wěn)定性的影響,特別是確定上反角等設(shè)計(jì)參數(shù)時(shí)。

參考文獻(xiàn)：

[1] Cestino E. Design of Very Long-Endurance Solar Powered UAV[D]. Politecnico di Torino, 2006

[2] Romeo G, Frulla G, Cestino E, Corsino G. HELIPLAT: Design, Aerodynamic, Structural Analysis of Long-Endurance Solar-Powered Stratospheric Platform [J]. Aircr, 2004, 41(6): 1505-1520

[3] Kirk Flittie, Bob Curtin. Pathfinder Solar-Powered Aircraft Flight Performance [R]. AIAA-1998-4446

[4] Youngblood J W, Talay T A. Solar-Powered Airplane Design for Long-Endurance, High-Altitude Flight [R]. AIAA-1982-0811

[5] 昌敏,周洲,鄭志成. 太陽能飛機(jī)原理及總體參數(shù)敏度分析[J]. 西北工業(yè)大學(xué)學(xué)報(bào), 2010, 28(5):792-796

Chang Min, Zhou Zhou, Zheng Zhicheng. Flight Principles of Solar-Powered Airplane and Sensitivity Analysis of Its Conceptual Parameters[J].Journal of Northwestern Polytechnical University, 2010, 28(5):792-796 (in Chinese)

[6] Noll T E, Brown J M, Perez-Davis M E, Ishmael S D. Investigation of the Helios Prototype Aircraft Mishap[R]. Tech Rep, NASA, 2004

[7] Patil M J, Hodges D H, Cesnik C E S. Nonlinear Aeroelasticity and Flight Dynamics of High-Altitude Long-Endurance Aircraft[R]. AIAA-1999-1470

[8] Patil M J, Hodges D H. Flight Dynamics of Highly Flexible Flying Wings [R]. Journal of Aircraft, 2006, 43(6): 1790-1798

[9] Jaworski J W. Nonlinear Aeroelastic Analysis of Flexible High Aspect Ratio Wing Including Correction with Experiment[D]. Duke University, 2009

[10] Zhang Jian, Xiang Jinwu. Nonlinear Aeroelastic Response of High-Aspect-Ratio Flexible Wings[J]. Chinese Journal of Aeronautics, 2009, 22(4): 355-363

[11] Crisfield M A, Galvanetto U, Jeleni? G. Dynamics of 3-D Co-Rotational Beams[J]. Computational Mechanics,1997,20:507-519

[12] Rafael Palacios, Joseba Murua, Robert Cook. Structural and Aerodynamic Models in Nonlinear Flight Dynamics of Very Flexible Aircraft[J]. AIAA Journal, 2010, 48(11): 2648-2659