飛機翼面氣動載荷壓力分布有限元節點插值方法研究

張建剛

航空工業第一飛機設計研究院，陜西 西安 710089

結構強度專業需要將氣動載荷準確施加在結構有限元模型上作為輸入，以進行下一步的設計，這是工程實踐中一項重要而繁瑣的工作。而目前結構強度專業使用的軟件還不具備此功能。為了提高計算精度、節約人力資源，需要對氣動載荷在有限元節點上加載的問題進行專門的研究。

氣動彈性專業計算中有類似的問題，氣彈計算每一迭代步需要將氣動載荷施加在結構節點上，通常的做法是在氣動網格上直接積分，形成氣動載荷的集中力，然后用樣條插值矩陣將該集中氣動力變換為結構網格點上的等效值[1～3]。等效值計算利用的是虛功原理，國內結構強度專業常采用“多點排”的算法[4]，以應變能最小和靜力等效作為約束條件，將氣動載荷分配到若干個有限元節點上。王專利[5]、徐建新[6]等利用該方法計算得到了翼面結構的有限元節點載荷。高尚君[7]等在該方法的基礎上，利用基于彈簧—懸臂梁模型最小變形能的分配方法，解決了多點排方法中氣動點與部分有限元點重合而導致的難以求解的問題。以上方法均是首先在氣動網格上積分得到集中力，然后再等效到有限元網格上。這些方法雖然可以有效解決工程中的一些問題，但仍有需改進之處。這是由于氣動網格和結構有限元網格的邊界通常是有差異的，直接在氣動網格上積分求得集中力的方法可能會造成氣動載荷的跨區域傳遞。另外，載荷專業提供的氣動載荷形式是一些離散氣動點上的壓力分布及翼面的總載、總矩。氣動點上的壓力值通常由試驗或理論計算得到，這些離散點的位置和有限元節點是不重合的，因此，需要由氣動點上的壓力值計算得到有限元節點上的壓力值。

本文分別采用薄板樣條插值函數和彎曲板單元插值形函數的方法進行有限元節點的壓力值插值，并將兩者的結果與試驗結果進行比較，以確定方法的有效性。

1 薄板樣條插值函數

薄板樣條插值方法是航空航天工程中廣泛使用的方法[8]。設在M維歐氏空間內域D上有n個矢量（每個矢量代表一個節點坐標）及其對應的函數值，即已知：

設矢量X的函數X=W（X）（稱為真實函數）是單值連續的，則對D上的任一矢量Xk（稱為插值點）必有且只有一個確定的函數Wk與之對應。由于在插值問題中真實函數是未知的，必須利用已知條件構造一個逼近函數，用逼近函數代替真實函數。在這里，逼近函數的形式如下：

式中：c1，c2，…，cM+1+n為待定系數；M為X的維數，由于本文中是平面問題，M=2；ε為給定的常數，稱為參量，對于一般平坦曲面ε=10-2～1，對于有奇性的曲面可取ε=10-5～10-6；xpi為第i個已知節點的第P維坐標；xpk為待求節點Xk的第P維坐標。

式（2）中的待定系數由無窮遠處的邊界條件及n個已知節點函數值來確定：

式中：xpi為第i個已知節點的第P維坐標；hj為對應于第j個節點的加權系數。

將上述方程組寫成矩陣形式如下：

記為：

式中：s為待定系數組成的矢量；f為已知函數值和0組成的矢量

求解式（5），可得到式（2）中的各待定系數，即逼近函數。對于任一待求節點，利用式（2）即可求得對應的函數值。

2 彎曲板單元插值

彎曲板單元插值方法借鑒了有限元方法，即在單元內，選擇簡單近似函數來分片逼近未知的求解函數，即分片近似，這是有限元法的創意和精華所在[10]。而整體區域上的解函數就是這些單元上的簡單近似函數的組合。

首先假設氣動點上的壓力分布類似于有限元理論的彎曲薄板的位移分布。將氣動點作為虛擬的有限元節點，將氣動網格作為虛擬的有限元網格，虛擬節點上的壓力值作為彎曲板單元的節點的位移值，真實的有限元節點是虛擬單元的內部點。根據有限元形函數插值理論，可由節點的位移值計算出單元內的任一點的位移值，此位移值即為真實有限元節點上的壓力值。本文中采用的有限元單元是彎曲薄板單元中含有9個位移參數的部分協調三角形單元[9]。

設三角形單元頂點1，2和3的直角坐標分別是（x1，y1)，（x2，y2)，和（x3，y3)，節點順序按逆時針方向排列。面積坐標是三角形的無因次坐標，其定義為：

式中：A是三角形的面積，A023、A031、A012是子三角形的面積。面積坐標只有兩個是獨立的，因為它們要滿足關系：

三角形的面積為：

式（6）可展開為：

式中：w是單元內部任意一點的位移；δe是單元節點位移；N是形函數。

形函數表達式為：

式中：b1=y2-y3；b2=y3-y1；b3=y1-y2；c1=x3-x2；c2=x1-x3；c3=x2-x1。

實際計算中，以氣動點的壓力作為虛擬單元的節點位移δe，代入式（11）中即可求得真實有限元節點上的壓力值。計算前需要先判斷有限元節點的位置在哪一個虛擬單元內。設某一個三角形單元的三個節點坐標（xi，yi），i=1，2，3，真實有限元節點坐標是（x，y），如圖1所示。假設三角形面積是S，（x，y）點和三角形三個頂點連成三個三角形，若則認為（x，y）點在三角形單元內部。

圖1 真實有限元節點與虛擬單元的位置關系Fig.1 The relation between true FEM node and virtual element

3 插值方法結果分析

3.1 壓力插值結果比較

按照以上介紹的兩種方法，以某飛機機翼主翼面上、下翼面某工況的載荷為例。將原始氣動點的壓力值和插值后有限元節點的壓力值進行比較，如圖2～圖5所示。這里的壓力指的是當地靜壓減去來流靜壓后的剩余量，即壓力系數和動壓的乘積。

圖2 上翼面薄板樣條插值結果與原始氣動點壓力比較Fig.2 Contrast between the thin plate spline result and the original aerodynamic press at up surface

圖3 上翼面彎曲薄板插值結果與原始氣動點壓力比較Fig.3 Contrast between bending element spline result and the original aerodynamic press at up surface

圖4 下翼面薄板樣條插值結果與原始氣動點壓力比較Fig.4 Contrast between the thin plate spline result and the original aerodynamic press at down surface

圖5 下翼面彎曲薄板插值結果與原始氣動點壓力比較Fig.5 Contrast between bending element spline result and the original aerodynamic press at down surface

3.2 翼面總載比較

計算得到有限元節點的壓力以后，等效節點載荷的計算比較簡單。設單元受法向分布載荷q的作用，對應的等效節點載荷列矢量為：

如果分布載荷在單元內是線性變化的，即：

式中：qi，qj和qm分別為q在節點i，j和m處的大小。

將式（13）代入式（12）得：

根據式（14）即可求得真實有限元模型上的等效氣動載荷，供結構強度專業進行計算。將該機翼此工況的上下翼面計算得到的節點力和力矩求和得到總載，并將結果與試驗實測的總載比較，結果見表1。

表1 兩種計算結果總載荷與試驗值比較Table1 Contrast between the total load of two method and test result

可以看出，兩種插值結果與原始壓力均符合的較好。相比較而言，彎曲薄板的結果更為光滑，奇點更少一些。表1中的總載對比也證明了彎曲薄板的計算結果總載符合的更好一些。此外，彎曲薄板插值方法的計算量也小一些，因為薄板樣條插值的計算過程中需要求解大型線性代數方程組。

3.3 兩種方法結果對比分析

兩種插值方法結果的差異本質上是由插值方法的不同引起的。若不考慮計算的舍入誤差，兩種方法求出的插值多項式都精確滿足插值條件，即在插值數據點沒有誤差，但利用插值多項式的目的正是要在這些數據點以外的地方求值。Runge現象[11]說明當插值的數據點的數目增大到10以后，誤差點隨著插值數據點數目的增大而增大。通常一塊翼面上的氣動點數以千計，薄板樣條插值插值結果可能會出現較大的誤差，尤其是當氣動點的某些點是壞點的時候，該壞點的影響區域也會較大。用相對較少的插值數據點作插值，可以避免大的誤差，但我們總希望將所有的數據點都用上，且所用數據點越多越好，解決這一矛盾的辦法是采用分段插值方法，即用分段多項式來代替單個多項式作插值。所謂分段插值是由一些在相互連接的區間上的不同多項式連接而成的一條連續曲線或曲面，這也是有限元方法分片近似的核心。

4 結束語

本文采用了樣條薄板插值和彎曲薄板插值兩種計算方法，完成了由氣動點壓力值求取有限元模型節點的壓力值。兩種方法的計算結果均與試驗結果符合較好，以往的工程實踐表明，二者都能較為圓滿地解決工程問題。相比較而言，彎曲薄板插值方法在計算結果準確性和計算效率方面更好。