考慮幾何非線性效應的大柔性太陽能無人機靜氣動彈性分析

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

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

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

1 大柔性結構幾何非線性靜力學建模

大柔性無人機機翼可以用柔性梁模型來描述[7-9],因此本文選擇空間CR梁單元建立大柔性太陽能無人機的結構模型。采用有限元技術對大柔性結構建模時,其計算精度主要取決于剛度矩陣和坐標系轉換矩陣對幾何非線性描述的精確性。CR法的平衡方程在單元坐標系內(nèi)建立,通過坐標轉換矩陣及其微分形式構造總體坐標系下的切線剛度矩陣,因而其核心是坐標轉換矩陣的建立及其微分形式的推導。

1.1 坐標系的建立

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

圖1 柔性結構變形示意圖

當結構位形確定時,單元節(jié)點的3個歐拉角及空間坐標隨之確定,按照式1旋轉總體坐標系可以得到節(jié)點坐標系。

ui(x)=Ti(x)Ixyz

(1)

Ti(x)是節(jié)點i的坐標系轉換矩陣,是3個歐拉角φx,φy,φz的函數(shù),其中:

對于2節(jié)點空間梁單元,采用(1)式可以得到節(jié)點坐標系ui和uj,假設其坐標轉換矩陣Ti和Tj對應的旋轉偽矢量為(θix,θiy,θiz)T和(θjx,θjy,θjz)T,參照文獻[11],以γ/2作為旋轉偽矢量,對節(jié)點坐標系ui進行旋轉得到平均構形坐標系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)=

那么參照文獻[11],空間梁單元的單元坐標系的定義為:

(3)

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

1.2 增量平衡方程推導

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

(4)

對(4)式進行微分得到位移增量從單元坐標系到總體坐標系的坐標轉換矩陣T:

δpl=Tδp

(5)

單元坐標系下的廣義節(jié)點力fl為:

fl=Klpl

(6)

式中:Kl是線性剛度矩陣,與線性梁單元的剛度矩陣一致。

則總體坐標系下廣義節(jié)點力向量f為:

(7)

(8)

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

KTΔp=Ft+Δt-ft

(9)

2 幾何非線性靜氣動彈性建模

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

(10)

圖2 氣動載荷插值示意圖

那么,靜氣動彈性運動方程為:

KTΔp=Qt-ft

(11)

收斂控制方程為:

(12)

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

3 代碼驗證及算例分析

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

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

圖3 懸臂梁幾何大變形

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

3.2 大柔性太陽能無人機非線性靜氣動彈性計算與結果分析

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

圖4 太陽能無人機幾何模型

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

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

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

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

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

機翼的彈性變形使氣動載荷向翼根轉移,這種幾何非線性效應改善了機翼結構的外載荷分布,對結構設計是有利的,而線彈性結構模型無法精確地模擬到這種趨勢,甚至得到相反的結論。

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

當翼尖撓度為22.5%的展長時,升力系數(shù)降低了12.4%,升阻比降低了11.8%。隨著變形增加,阻力系數(shù)基本不變,而升力系數(shù)減小比較快,因而升阻比的變化趨勢也是減小的。這主要是因為幾何大變形改變了氣動載荷的作用方向,對全機升力的影響較大。采用線彈性結構模型時,氣動力是隨彎曲變形的增加而單調增加的[7],因此采用線性結構模型與非線性結構模型研究較大彈性變形對氣動力的影響時,會得到不同的結論。

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

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

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

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

表2 橫航向氣動特性

4 結 論

本文結合CR有限元法及計算流體力學理論,對大柔性太陽能無人機的靜氣動彈性問題進行了深入的研究,得到以下結論:

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

2)幾何非線性效應改善了氣動載荷的重新分布,有利于機翼結構設計;但同時會降低無人機的氣動性能。

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

參考文獻：

[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] 昌敏,周洲,鄭志成. 太陽能飛機原理及總體參數(shù)敏度分析[J]. 西北工業(yè)大學學報, 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