太陽能無人機機翼顫振動力學建模與分析

冒 森，張 斌，肖良華，陳 斌，王 玨

(成都飛機工業(集團)有限責任公司，成都 610073)

0 引言

太陽能飛機是采用太陽能為能源的飛行器。相對比常規布局飛機，太陽能飛機超長的航時和超高的飛行高度使其具有無與倫比的優勢，而且契合綠色環保理念，因此得到國內外研究人員的廣泛關注。此類飛機的設計，以能源設計為中心，采用超大展弦比的機翼及極低的翼載設計(翼載較常規動力飛機的低2個量級左右)，例如“太陽神”太陽能無人機結構面密度為 3.2 kg/m2[1]，而現階段最長航時常規動力無人機“全球鷹” 的結構面密度高達 53 kg/m2[2]，因此太陽能飛機工程實現難度非常高。從全機重量分配來看，結構重量和蓄電池重量分別占全機重量的35%以上，且蓄電池重量隨總重的增加而增大，因此結構的輕量化是實現全機減重的有效手段。

輕量化結構和大展弦比構型使得太陽能無人機顫振、突風響應等氣動彈性問題非常顯著。其中顫振是最重要的一類結構動強度問題。顫振現象是是機體結構在受到彈性力、慣性力、氣動力共同作用下可能出現的一種自激振動，是一種飛機的固有屬性[3]。太陽能無人機因顫振而發生的事故最早被是2003 年 6 月，美國“太陽神”無人機經過多次成功飛行后，在夏威夷附近海域進行低空飛行時，機體突然發生向上彎曲，進而出現俯仰振蕩發散，最終導致墜毀[4]。距現在較近的一次太陽能飛機因顫振而發生事故的是上海奧科賽飛機公司設計的太陽能飛機2號機。此飛機的翼展為15米，在2017年7月的一次飛行過程中，在飛行高度為204.2 m、飛行速度為92 km/h時發生了顫振，最終也導致了飛機墜毀[5]。常規飛行器設計中對于結構顫振問題，一般采取試驗與仿真相結合的方式進行解決，在風洞中可直接采用1：1的結構模型來進行試驗。而太陽能無人機的翼展遠大于現有的風洞尺寸，無法直接開展顫振風洞實驗，需要設計并加工縮比模型，而且縮比模型的重量縮比量是尺寸縮比的三次方[6]，由于太陽能飛機的面密度極低，導致縮比模型加工難度很高，如此苛刻的條件讓此類飛機基本不能開展全模的顫振試驗[7]。因此在方案設計階段，通過仿真計算的手段準確得到太陽能無人機的顫振速度并且給出防顫振設計規律具有重要意義。

目前在大展弦比飛機的顫振分析方面，國內外學者已經開展了廣泛而深入的工作，太陽能無人機在結構布局上可歸屬于大展弦比常規布局飛機，現有的針對大展弦比飛機的氣動彈性分析方法也適用于太陽能飛機。文獻[8]和[9]首次使用動力學線化理論，將大變形下的機翼非線性有限元分析進行線性化處理，針對大展弦比機翼，分析了機翼的非線性靜氣彈和動氣彈問題。文獻[10]中在計算大展弦比機翼的氣動載荷時，考慮曲面效應，應用了推廣的三維升力線理論，并且在結構/氣動耦合計算中首次使用曲面樣條插值理論，從而計算出氣動力計算所必需的網格節點位移，在此基礎上對一個簡單的魚刺模型進行了相應計算，計算結果證明了該方法能很好的解決大展弦比機翼靜氣動彈性問題。文獻[11]中則是國內第一個針對太陽能無人機開展氣動彈性分析。首先建立太陽能無人機的結構有限元模型和氣動計算模型，采用MSC.Nastran軟件計算了靜氣動彈性變形和顫振速度。通過對計算結果的分析，總結了彈性變形程度對機翼上升力分布的影響，得到了不同飛行條件下的顫振速度和頻率。文獻[12]針對國內翼展為15米的太陽能飛機，進行了該太陽能飛機的動力學反向建模，并以所建模型為基礎，進行顫振計算分析，最終建立了一套完整的關于太陽能飛機顫振分析的設計方法。文獻[13]則在理論層面上更深層次地推導了太陽能飛機中機翼在受力變形后的切線剛度矩陣和質量矩陣，在不同飛行工況靜氣彈平衡位置使用準模態假設，采用P-K法計算了該飛機的在不同飛行載荷下的非線性顫振問題。文獻[14]創建了一套使用大變形插值理論計算考慮幾何非線性效應時的氣動彈性分析框架，并以某太陽能無人機機翼縮比模型為對象使用此方法進行仿真分析，之后又開展地面模態試驗與風洞試驗，結果表明使用該方法進行仿真計算不僅可以顯著提高計算效率，而且能夠有效降低常規方法中載荷選擇時產生的計算偏差。

可以看出目前針對太陽能無人機的顫振分析研究很多也非常深入，但是上述大多數的研究是根據確定性的結構參數進行正向結構建模和分析，且采用的模型都是簡化的魚刺模型。實際上在太陽能飛機的方案設計階段，是無法開展地面GVT(ground vibration test)試驗獲得機翼的準確的剛度數據和重量數據，難以建立準確的魚刺模型來進行顫振分析。因此在設計階段，現有的方法不能很好的解決太陽能無人機的氣動彈性建模和顫振分析的問題。

此外，太陽能無人機區別于傳統飛機的一大特點在于蒙皮采用薄膜結構形式，體現在結構有限元中即為將蒙皮采用剛度相對較小的板單元[15-16]進行模擬，這會導致初始有限元模型在模態分析時出現大量蒙皮單元局部模態，在計算較多階模態的情況下仍無法獲取足夠多的整體模態信息。此外，全機復雜有限元模型的自由度一般較多(通常為數十萬自由度)，且包含大量用于部件間或集中質量點連接的MPC單元，因此有必要建立一種能夠準確反映初始有限元模型動力學特性的工程化建模方法。

另一方面，由于具有明顯彈性軸的太陽能機翼通常呈現出“長且柔”的懸臂梁特點，在機翼上布置吊艙可以起到氣動力卸載的作用，但因機翼本身已具有較低的固有頻率，系掛吊艙可能會進一步降低機翼的固有頻率，使得機翼顫振速度發生大的改變[17]。所以研究太陽能機翼上吊艙不同位置的分布對于顫振特性影響也是穩定性分析中不可缺少的一環的。

針對太陽能飛機的獨特之處，本文提出一種面向方案階段的太陽能無人機機翼顫振分析辦法，以結構靜力模型為基礎，通過工程化處理辦法，開展太陽能無人機機翼氣動彈性建模和顫振研究，在確保消除局部模態的同時不改變全機模態分布特性的基礎上，計算不同機翼吊艙布置下機翼顫振速度，并以此為依據，總結提出太陽能飛機防顫振設計規律。

1 動力學建模

1.1 結構設計

該太陽能無人機機翼結構如圖1所示。采用常規布局的單梁式機翼。主梁為變截面梁，翼肋是非主承力件。除連接起落架吊艙的肋為金屬件外，其余翼肋均為一體成型復材件。能源/動力/起落架吊艙位于單邊機翼展向10 m處。吊艙用于連接推進系統、起落架和安裝部分儲能電池。

圖1 機翼結構設計圖

蒙皮采用抗撕裂薄膜。結構材料以復合材料為主，局部使用金屬材料，金屬材料選擇常規的鋁合金。機翼主梁和接頭選用中溫固化環氧樹脂單向碳纖維預浸料，超輕質肋和泡沫夾芯結構的復材選用中溫固化環氧樹脂單向碳纖維預浸料。

1.2 有限元建模

對于該機翼使用商業軟件MSC.Patran開展有限元建模，有限元模型中按照實際情況建出，如圖2所示。

圖2 太陽能無人機機翼有限元模型

機翼肋緣條和腹板均按照實際數模建出，如圖3所示。

圖3 梁和肋有限元建模

太陽能電池陣鋪設在機翼上蒙皮，為模擬真實情況，有限元模型中在機翼上蒙皮單元每一個節點建立質量元，平均分配太陽能電池陣的重量，如圖4所示。

圖4 太陽能薄膜集中質量單元布置

1.3 建模分析

有限元建模完成后，使用MSC.Nastran的101模塊進行靜力學分析，檢驗機翼強度。平飛狀態的機翼氣動力采用面壓的方式加載在機翼蒙皮上，氣動力為CFD由計算得到，均以外部載荷的形式加載到結構。機翼結構變形如圖5所示。從靜力學角度可以看出，蒙皮出現了失穩。

圖5 機翼受氣動載荷變形圖

運用MSC.Nastran 103模塊開展有限元模態計算，檢驗建模準確性。計算結果發現在從4階模態之后，頻率2.0 Hz到2.5 Hz之間出現大量的局部模態，如圖6所示。這些局部模態，屬于蒙皮失穩模態。局部模態極大的影響顫振計算結果，造成顫振計算結果不準確。

圖6 蒙皮局部模態

需要注意的是，在開始進行模態計算中發現扭轉頻率偏低，排查原因是扭轉剛度存在問題，在有限元模型中對于主梁采用桿單元進行模擬，肋采用板單元進行模擬，在梁與肋之間采用單點連接，但由于主梁幾何尺寸比較大，僅簡化成桿單元和板單元之間采用單點連接的形式會造成模型扭轉剛度低于實際結構。

2 氣動彈性建模修正

2.1 工程化處理方法

為了體現梁肋連接的抗扭能力，在肋腹板處挖出圓孔，將腹板圓孔周圍節點與機翼主梁節點采用MPC固聯，如圖7所示。

圖7 梁和肋有限元建模優化

太陽能機翼的薄膜蒙皮只承受拉力，不能承受壓力，而且在安裝時是進行預緊安裝，蒙皮失穩只存在于仿真計算中，實際情況中不存在蒙皮失穩，因此在顫振計算中需要排除蒙皮局部模態的影響。現階段國內外用于顫振分析的商業軟件有ZAERO和MSC.Nastran，采用ZAERO計算顫振需要MSC.Nastran計算的模態信息，而MSC.Nastran無法針對這種情況直接進行顫振計算。另一種處理局部模態的辦法是通過MPC約束蒙皮自由度，但是也會增加機翼剛度，造成顫振計算不準確。

本文提出從工程處理角度，對失穩處的蒙皮，采用增加殼單元面外剛度的方法消除局部模態。在薄殼單元有限元計算中，單元面內剛度矩陣Lm和面外剛度矩陣Lb如下[18]：

(1)

(2)

薄板類殼單元的局部模態通常表現為單個單元的彎曲或鼓包，如果增加殼單元的面外剛度，相當于提高了殼單元抵抗彎曲或鼓包變形的能力，這樣即可達到消除局部模態的目的。同時，對于全機而言，帶有蒙皮單元的機翼或尾翼部件整體模態(彎曲、扭轉等)主要受蒙皮單元面內剛度影響，所以上述增大面外剛度的修改方式不會對全機整體模態產生太大影響。Nastran中可通過修改殼單元的PSHELL屬性，增大彎曲剛度的比例，消除局部模態[19-20]。對于常規飛機的復材蒙皮或者金屬蒙皮，失穩的處理辦法一般情況是也是提高蒙皮剛度進行處理，因此具有一定可信度。

2.2 模態計算

工程方法處理之后，機翼模態結果如圖8所示。表 1列出了修改剛度前后全機模型模態分析所得前8階模態頻率和振型。可以看到，修改剛度后模態分析結果與原始模型保持高度一致性，各階模態頻率誤差均在1%以內。

表1 原始模型與修正后模型的模態結果對比

圖8 機翼模態圖

3 顫振計算

頗振分析是用來確定氣動彈性系統的動穩定性問題，一般情況下氣動彈性顫振方程可以寫為：

(3)

(4)

其中：p=ω(γ±i)=2πf(γ±i)為方程復特征值，γ為衰減率，g=2γ為結構阻尼比，QIQR為廣義氣動力阻尼矩陣和廣義氣動力剛度矩陣。選擇MSC.Nastran的SOL145計算模塊開展顫振分析，由于高階模態局部模態偏多，因此值截取中前八階模態參與顫振計算。

P-K法理論上可進行任意運動形式下的顫振計算分析，但是由于在計算的過程中使用的Theodorson 氣動力仍然是結構在諧振蕩下產生的氣動力，因此在顫振點出的分析結果是較為準確的，而非臨界點處得到的計算結果僅能夠作為參考數據，并不一定準確[21]。使用P-K法進行顫振分析雖然計算量較大，但是求得的結果與飛行顫振試驗和風洞試驗結果比較接近，是一個近似真實阻尼解，而且可以得到不同動壓下的顫振速度。P-K法具體的計算流程是在給定動壓條件下，以來流速度v作為自變量，在顫振計算時對指定主要結構模態根據提前設定好的一系列減縮頻率點反復進行迭代計算，最終得到相應的氣動力，特征根和特征值。本文中根據MSC.Nastran的SOL145模塊得到的計算結果，繪制對應的V-g和V-f圖。

在顫振分析中，還需要指定氣動網格的劃分情況。在本文計算中氣動網格劃分情況見圖9，將機翼劃分4個氣動區域，內外翼段分開處理，其中內翼面劃分為3個區，副翼單獨劃分為1個區，外翼段1個區，氣動網格總數為640個，由1.3節的有限元分析可得，蒙皮承受氣動載荷會發生失穩，所以優化插值點選擇，只選擇在主梁和翼肋上的結構節點作為插值點。圖10所示為顫振分析中前5階結構模態插值至氣動面上的結果，表明插值信息無異常。

圖9 氣動網格

圖10 氣動面插值情況

3.1 顫振計算結果

計算得到顫振計算結果如圖11、圖12所示。

圖11 V-G

圖12 V-F

從V-G圖看出兩個模態發生穿越，分別是二階彎曲和三階彎曲，二階彎曲穿越速度為26 m/s，三階彎曲穿越速度為61 m/s，所以顫振速度為26 m/s。水平一彎模態也發了穿越，但是阻尼小于0.001，屬于小阻尼模態不是顫振。

從V-F圖可以看出二階彎曲和三階彎曲對應的顫振是傳統的爆發性彎扭顫振。原理是扭轉模態正好處于二階彎曲模態與三階彎曲模態之間，隨著速度的增加這些模態頻率接近發生相互耦合，出現顫振。一階彎曲模態在54 m/s處頻率降為0，發生扭轉發散現象，扭轉發散速度為54 m/s。一階彎曲模態與扭轉模態頻率差距比較大，兩者之間沒發生顫振耦合。

綜上可得機翼的顫振速度為26 m/s，顫振速度計算值偏大，原因在于增加了機翼剛度。但是修改后模態頻率相差小于1%，真實的顫振計算結果不會出現顯著變化。機翼的顫振速度小于扭轉發散速度，所以應該優先考慮顫振問題。該無人機設計最大飛行表速14 m/s，考慮顫振安全余量15%，所以該機翼的顫振計算結果滿足顫振設計要求。

3.2 吊艙布置對于顫振速度的影響

現有研究已表明，在機翼吊裝外掛在某些程度上會改變機翼的顫振特性，分析機翼外掛對顫振的影響成為飛機尤其是具有大展弦比柔性機翼飛機氣動彈性分析的一個必要過程。常見帶外掛的戰斗機通常都采用小展弦比三角翼，翼展較小且剛度高，因此外掛對于顫振影響不是很突出，但是類似于“捕食者”這種大展弦比無人機在帶外掛后對顫振速度影響較大，因為外掛會現在影響大展弦比機翼剛度分布。此外對于小展弦比機翼，面內的彎曲變形較小，對氣動彈性的影響可以忽略，但是在大展弦比機翼分析中機翼面內的彎曲變形是不可忽略的，因為展弦比的增加會使機翼的平面內特性更加明顯，面內變形的程度也極大增加。且對于太陽能飛機這種結構來說，由于柔性增加且翼載減少，在加掛吊艙后，吊艙的重量和相對于機翼的位置對于機翼的模態會有很大的改變，因此會顯著影響到機翼的顫振速度。

本文中模型中吊艙總質量約為40 kg，與機翼結構質量相當，如此大的集中質量，吊掛在機翼的某一位置，會對機翼的顫振速度產生很大的影響。吊艙布置對于顫振的影響可分為兩部分：一部分是外伸的發動機對于顫振速度的影響；另一部分為吊艙的展向占位于顫振速度的影響。

首先保持吊艙的展向位置和發動機質量不動，更改發動機連桿長度，得到其與機翼顫振速度的關系如圖13。

圖13 發動機連桿長度與機翼顫振速度關系

保持吊艙的展向位置不動，發動機連桿長度不變，修改發動機質量得到其與機翼顫振速度的關系見圖14。

圖14 發動機重量與機翼顫振速度關系

從以上算例可以看出，發動機連桿越長，發動機越重，機翼的顫振速度越高。從模態角度看，這是因為發動機越重距離主梁越遠，機翼的扭轉模態越低，機翼越不容易發生顫振。從另一方面說明，該機翼重心在主梁后方，可以通過增加前緣配重，使得重心前移，減小扭轉重心與重心的距離，從而提高顫振速度。值得注意的是發動機是由總體需求確定的，重量不可以隨意更改，可以通過在相同位置處增加配重，提高機翼顫振速度。

最后保持吊發動機連桿長度和發動機重量不變，修改吊艙站位，得到其與機翼顫振速度的關系見圖15。

圖15 吊艙位置與機翼顫振速度關系

可以看出吊艙越靠近機身，顫振速度越高。這是因為吊艙內移，會使得機翼扭轉模態頻率顯著提高，距離相鄰的彎曲模態頻率差距增加，不容易發生耦合，因此更不容易發生顫振。雖然吊艙內移可以顯著提高機翼顫振速度，但是不是越內移越好。吊艙離機身太近會使得其對機翼的氣動力卸載作用減弱，導致機翼翼尖變形過大，對于飛行控制和安全帶來不利影響。

從以上分析中可以得出結論，二階水平彎曲是太陽能機翼需要特別注意的一種模態，這種模態與機翼的一階扭轉模態的耦合可能導致顫振速度降低。這種類型的顫振速度隨著機翼彎曲變形的增加而減小。在設計太陽能無人機時，必須特別注意扭轉剛度的影響，避免二階彎曲、一階扭轉和三階彎曲的頻率接近，以確保二階彎曲模態中涉及的顫振模態具有足夠高的顫振速度。

4 結束語

本文以國內某40米翼展太陽能無人機為研究對象，通過工程化處理方法，運用商業軟件MSC.Nastran開展氣動彈性建模修正和顫振分析，得出以下結論：

1)本文通過修改梁與肋的傳力結構、修改蒙皮屬性消除局部模態和優化插值點選擇等工程處理方法修正模型，并通過仿真計算得到太陽能無人機機翼顫振速度為26 m/s。雖然采用MSC.Nastran計算得到的顫振速度較實際結果偏大，但是計算結果遠高于設計安全值，可以認為該機翼滿足顫振設計要求。

2)對于本文中的機翼，可以通過增加發動機連桿長度，給發動機增加配重，和吊艙內移的方式，改變模型的扭轉剛度，從而提高機翼的顫振速度。

3)后續研究中可以開展地面模態試驗，進行有限元模型修正，從而獲得更準確的顫振速度。未來則可考慮開展太陽能飛機的地面車載滑跑試驗，保證試飛安全。