直升機艙門不均勻分布氣動載荷加載方法研究

李玉蓮 袁勝弢

(哈爾濱飛機工業集團有限責任公司飛機設計研究所，哈爾濱150066)

應力分析是直升機結構強度設計最重要的工作內容之一，隨著商業分析軟件的快速發展，有限元已經成為復雜結構應力分析中最重要的工具[1]。結構強度有限元分析過程包括有限元前處理、求解計算和后處理。對于直升機艙門結構強度有限元分析前處理工作，氣動載荷是分析設計工作中非常重要的載荷設計輸入[2]。在國內航空界大部分科研院所，氣動載荷分布由載荷專業計算，結構強度分析專業在有限元的前處理工作中將該載荷模擬并施加在有限元模型節點(單元)上，以進行有限元的求解計算工作。在工程設計中，將不均勻分布氣動載荷施加到有限元節點上是一項非常繁瑣但又很重要的工作。截至目前，航空航天行業廣泛應用的有限元分析軟件如Patran/Nastran，Abaqus等，盡管提供了強大的有限元分析計算功能，但對于不均勻分布面載荷還不能自動實現加載。人工手動對有限元節點或者單元進行逐個加載，滿足不了實際工程設計的需求，因為直升機結構強度設計中的載荷工況較多，而且有限元細節模型的節點數目非常龐大，逐個單元或者節點的加載方式嚴重影響載荷施加的效率以及質量。如何提高有限元分析前處理工作中加載的效率與質量成為國內外專家學者重視的問題。

為了解決有限元分析前處理中施加不均勻分布氣動載荷這一難題，國內外出現了很多算法。目前航空領域針對翼面結構的氣動載荷分布，常用的算法是“多點排”方法[3]，該方法以靜力等效以及應變能最小為約束條件，將每一種氣動載荷分配到結構有限元分析模型的節點上。林小夏等[4]提出基于特征函數分布對不均勻載荷進行加載。王專利[5]、徐建新等[6]使用了“多點排”的計算方法。高尚君等[7]提出基于彈簧懸臂梁模型最小變形能的氣動載荷分配方法。張建剛等[8]針對翼面結構提出薄板樣條插值函數來計算節點壓力值，這些算法能解決實際中的問題，但需要在結構有限元模型與氣動模型之間建立一種載荷數據傳遞關系，這種傳遞關系通過樣條矩陣來實現。這些方法都是基于數值插值分析，即每一次計算均要調用原始數據，增加了算法的時間以及復雜度，而針對直升機結構氣動載荷如何實現快速加載、等效分配的方法比較少。

本文以直升機艙門強度設計的工程實際需求為出發點，根據已知氣動模型節點的壓力值，通過數學分析最小二乘法，擬合出壓力分布曲面，由該曲面擬合函數計算得到結構有限元模型每一節點的壓力值，并對數值回歸性進行分析。通過該擬合壓力分布曲面函數的方法，可以準確、快速地將直升機艙門的不均勻分布氣動載荷分配施加到有限元節點上，以供直升機結構強度工程師參考。

1 直升機艙門氣動載荷

現代飛行器結構建模和氣動模型普遍采用有限單元法[9]，兩種模型劃分是為了各自的分析任務而獨立進行的。由于兩種模型的分析目的不同，建模方式也存在很大的不同，以某型號直升機艙門為例，艙門的氣動模型與結構分析有限元模型節點存在很大的差異。艙門氣動模型的節點數為1892個，氣動壓力分布計算結果為1892組離散數據；而有限元模型為了保留結構邊界，區分玻璃區域與結構區域，模型劃分方法以及節點數目與氣動模型均存在不同，兩種模型的有限元節點差異對比如圖1所示。

圖1 氣動節點與結構有限元節點差異示意圖

在艙門的結構強度分析有限元模型中，要實現施加不均勻分布的氣動載荷，需要在結構有限元模型與氣動模型之間建立一種載荷數據傳遞關系，工程上一般采用數值插值法或數據擬合法。數據擬合法根據氣動載荷分布數據擬合出相對簡單的數學模型，不需要每一個有限元節點的計算均調用原始氣動分布數據，降低了算法時間和復雜度，擬合函數更逼近真實函數，本文中這種傳遞關系通過數學分析，采用最小二乘法擬合來實現[10]。

2 最小二乘法基本理論

最小二乘法是一種基于已知離散點而擬合出近似函數的方法，使用離散點數據擬合得到曲面。最小二乘法具有擬合精度高、通用性強的特點。已知離散點函數(稱為真實氣動載荷)是單值連續的，則對空間區域D上的任何一個向量坐標點來說，理論上必有且只有一個確定的函數與之對應。由于在計算過程中真實函數是不知道的，任一向量所對應的函數值無法直接求出。為此，必須利用已知條件(有限元的節點以及氣動載荷信息)構成一個逼近函數，用它來代替真實函數，借以計算任一向量所對應的函數的近似值。

假設空間區域D內，已知數值的n個數據點集p i，其中i=1,2,···,n，其坐標可以表示為(x i,y i,z i)，在n個節點z i=(x i,y i)有已知值

已知數值用矢量表示為

設擬合函數φ(z)=[b1(z)b2(z)···b m(z)]。已知數值用矢量表示為

其中b m(z)為m個線形無關的基函數，λ(z)為m維系數矢量，是基函數的待定系數。最小二乘法即按照偏差平方和最小原則選取擬合曲線，在計算區域內，由式(3),利用加權最小二乘法構成二次形式

其中

基函數b m(z)是已知的，為求得待定系數，式(4)中對λ(z)一次求導等于0，求得

從而求得擬合近似函數為

3 艙門氣動分布擬合

直升機艙門的氣動載荷為不均勻的分布載荷，載荷專業提供的艙門壓力為氣動模型離散點的分布，在某型號直升機設計中，艙門某種狀態的氣動壓力分布等高線分布圖如圖2所示。從圖中可以看出，對于直升機艙門結構，沿著艙門縱向坐標X以及高度方向坐標Z的不同，氣動壓力分布呈現為非均布載荷。

圖2 某型機艙門氣動分布

艙門氣動載荷擬合函數具有最小二乘法的性質，得到直升機在空中各種不同飛行狀態時壓力分布的逼近曲面函數。在計算過程中，待定系數的快速求解采用數學分析軟件Mathematica中的程序進行[11]。

分析計算中，有限元節點的坐標值參考坐標系為機身全局坐標系。艙門機身軸向向后為X坐標正向，機身高度向上為Z坐標正向。最小二乘法擬合函數可以是一次、二次函數等所有初等函數，在艙門的壓力分布中，分別對平面分布和曲面分布進行多次擬合以及回歸性分析，以計算校核逼近函數的精度。其中線性基擬合時，式(3)中m取3，基函數表達式為式(9)；二次基擬合時，式(3)中m取6，基函數表達式為式(10)；立方基擬合時，式(3)中m取10，基函數表達式為式(11)。

當m=3時，艙門的壓力分布載荷呈平面分布，基函數取二次以上時，艙門的壓力分布為曲面分布。根據已知氣動點的壓力值求出逼近函數，進而求得有限元節點壓力值。有限元軟件Patran中有Field(域)加載功能，通過最小二乘法擬合，可快速實現不均勻分布載荷的施加處理過程。

4 艙門氣動載荷擬合結果分析

4.1 數學分析與擬合

某型號直升機艙門的非均布氣動分布載荷，當采用不同的基函數進行最小二乘法擬合時，確定系數如表1所示。圖3～圖6是不同基函數所計算得出的擬合函數與艙門原氣動點壓力值對比分析，從圖4與圖3的對比結果可以看出，隨著擬合函數系數的增高，擬合后的壓力曲面和擬合前的節點壓力分布趨勢一致，吻合程度逐漸提高，結合表1數據分析，從線性基到二次基的確定系數提高量最大，二次基后隨著擬合函數次數提高，確定系數的數值區域增加得較為平緩。

表1 擬合誤差分析表

圖3 艙門氣動載荷數據與一次擬合曲面

圖4 艙門氣動載荷數據與二次擬合曲面

圖5 艙門氣動載荷數據與三次擬合曲面

圖6 艙門氣動載荷數據與四次擬合曲面

根據以上分析，二次基的擬合結果較為優異，當采用二次基函數擬合時，對氣動節點進行擬合精度分析，計算結果見表2，擬合結果與原始數據的單點誤差最大值小于1%，最小二乘法的曲面擬合結果具有較高精度，能夠滿足工程上的設計要求。

表2 擬合結果與氣動節點壓強對比分析表

4.2 擬合結果分析

從圖6中可以看出，插值擬合的壓力曲面和插值前的氣動點壓力分布吻合程度很好，通過這種數據傳遞方式，可以快速實現由氣動模型的散點壓力分布求出有限元節點的壓力分布，進而通過逼近函數實現對有限元模型的快速加載。

5 結論

本文根據工程實際需求，基于氣動離散點壓力分布，通過最小二乘法擬合函數，可以實現直升機艙門壓力載荷擬合加載，提高有限元分析中加載的效率，擬合加載的精度較高。

Mathematica程序中的數學理論清晰，編程擬合插值簡單，可以有效解決工程中的實際問題。根據不同的工況，擬合逼近函數隨著氣動離散點壓力分布的改變而改變，尤其針對直升機的飛行狀態多變、壓力分布多變的情況，插值擬合函數具有較好的適應性，對多工況可以實現快速方便地加載，對于主要承受不均勻分布的氣動載荷的結構，提高了強度有限元分析的效率。