李玉蓮 袁勝弢

(哈爾濱飛機(jī)工業(yè)集團(tuán)有限責(zé)任公司飛機(jī)設(shè)計(jì)研究所，哈爾濱150066)

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

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

本文以直升機(jī)艙門強(qiáng)度設(shè)計(jì)的工程實(shí)際需求為出發(fā)點(diǎn)，根據(jù)已知?dú)鈩幽Ｐ凸?jié)點(diǎn)的壓力值，通過數(shù)學(xué)分析最小二乘法，擬合出壓力分布曲面，由該曲面擬合函數(shù)計(jì)算得到結(jié)構(gòu)有限元模型每一節(jié)點(diǎn)的壓力值，并對數(shù)值回歸性進(jìn)行分析。通過該擬合壓力分布曲面函數(shù)的方法，可以準(zhǔn)確、快速地將直升機(jī)艙門的不均勻分布?xì)鈩虞d荷分配施加到有限元節(jié)點(diǎn)上，以供直升機(jī)結(jié)構(gòu)強(qiáng)度工程師參考。

1 直升機(jī)艙門氣動載荷

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

圖1 氣動節(jié)點(diǎn)與結(jié)構(gòu)有限元節(jié)點(diǎn)差異示意圖

在艙門的結(jié)構(gòu)強(qiáng)度分析有限元模型中，要實(shí)現(xiàn)施加不均勻分布的氣動載荷，需要在結(jié)構(gòu)有限元模型與氣動模型之間建立一種載荷數(shù)據(jù)傳遞關(guān)系，工程上一般采用數(shù)值插值法或數(shù)據(jù)擬合法。數(shù)據(jù)擬合法根據(jù)氣動載荷分布數(shù)據(jù)擬合出相對簡單的數(shù)學(xué)模型，不需要每一個有限元節(jié)點(diǎn)的計(jì)算均調(diào)用原始?xì)鈩臃植紨?shù)據(jù)，降低了算法時(shí)間和復(fù)雜度，擬合函數(shù)更逼近真實(shí)函數(shù)，本文中這種傳遞關(guān)系通過數(shù)學(xué)分析，采用最小二乘法擬合來實(shí)現(xiàn)[10]。

2 最小二乘法基本理論

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

假設(shè)空間區(qū)域D內(nèi)，已知數(shù)值的n個數(shù)據(jù)點(diǎn)集p i，其中i=1,2,···,n，其坐標(biāo)可以表示為(x i,y i,z i)，在n個節(jié)點(diǎn)z i=(x i,y i)有已知值

已知數(shù)值用矢量表示為

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

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

其中

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

從而求得擬合近似函數(shù)為

3 艙門氣動分布擬合

直升機(jī)艙門的氣動載荷為不均勻的分布載荷，載荷專業(yè)提供的艙門壓力為氣動模型離散點(diǎn)的分布，在某型號直升機(jī)設(shè)計(jì)中，艙門某種狀態(tài)的氣動壓力分布等高線分布圖如圖2所示。從圖中可以看出，對于直升機(jī)艙門結(jié)構(gòu)，沿著艙門縱向坐標(biāo)X以及高度方向坐標(biāo)Z的不同，氣動壓力分布呈現(xiàn)為非均布載荷。

圖2 某型機(jī)艙門氣動分布

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

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

當(dāng)m=3時(shí)，艙門的壓力分布載荷呈平面分布，基函數(shù)取二次以上時(shí)，艙門的壓力分布為曲面分布。根據(jù)已知?dú)鈩狱c(diǎn)的壓力值求出逼近函數(shù)，進(jìn)而求得有限元節(jié)點(diǎn)壓力值。有限元軟件Patran中有Field(域)加載功能，通過最小二乘法擬合，可快速實(shí)現(xiàn)不均勻分布載荷的施加處理過程。

4 艙門氣動載荷擬合結(jié)果分析

4.1 數(shù)學(xué)分析與擬合

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

表1 擬合誤差分析表

圖3 艙門氣動載荷數(shù)據(jù)與一次擬合曲面

圖4 艙門氣動載荷數(shù)據(jù)與二次擬合曲面

圖5 艙門氣動載荷數(shù)據(jù)與三次擬合曲面

圖6 艙門氣動載荷數(shù)據(jù)與四次擬合曲面

根據(jù)以上分析，二次基的擬合結(jié)果較為優(yōu)異，當(dāng)采用二次基函數(shù)擬合時(shí)，對氣動節(jié)點(diǎn)進(jìn)行擬合精度分析，計(jì)算結(jié)果見表2，擬合結(jié)果與原始數(shù)據(jù)的單點(diǎn)誤差最大值小于1%，最小二乘法的曲面擬合結(jié)果具有較高精度，能夠滿足工程上的設(shè)計(jì)要求。

表2 擬合結(jié)果與氣動節(jié)點(diǎn)壓強(qiáng)對比分析表

4.2 擬合結(jié)果分析

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

5 結(jié)論

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

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