張建剛

航空工業(yè)第一飛機(jī)設(shè)計(jì)研究院，陜西 西安 710089

結(jié)構(gòu)強(qiáng)度專業(yè)需要將氣動載荷準(zhǔn)確施加在結(jié)構(gòu)有限元模型上作為輸入，以進(jìn)行下一步的設(shè)計(jì)，這是工程實(shí)踐中一項(xiàng)重要而繁瑣的工作。而目前結(jié)構(gòu)強(qiáng)度專業(yè)使用的軟件還不具備此功能。為了提高計(jì)算精度、節(jié)約人力資源，需要對氣動載荷在有限元節(jié)點(diǎn)上加載的問題進(jìn)行專門的研究。

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

本文分別采用薄板樣條插值函數(shù)和彎曲板單元插值形函數(shù)的方法進(jìn)行有限元節(jié)點(diǎn)的壓力值插值，并將兩者的結(jié)果與試驗(yàn)結(jié)果進(jìn)行比較，以確定方法的有效性。

1 薄板樣條插值函數(shù)

薄板樣條插值方法是航空航天工程中廣泛使用的方法[8]。設(shè)在M維歐氏空間內(nèi)域D上有n個(gè)矢量（每個(gè)矢量代表一個(gè)節(jié)點(diǎn)坐標(biāo)）及其對應(yīng)的函數(shù)值，即已知：

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

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

式（2）中的待定系數(shù)由無窮遠(yuǎn)處的邊界條件及n個(gè)已知節(jié)點(diǎn)函數(shù)值來確定：

式中：xpi為第i個(gè)已知節(jié)點(diǎn)的第P維坐標(biāo)；hj為對應(yīng)于第j個(gè)節(jié)點(diǎn)的加權(quán)系數(shù)。

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

記為：

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

求解式（5），可得到式（2）中的各待定系數(shù)，即逼近函數(shù)。對于任一待求節(jié)點(diǎn)，利用式（2）即可求得對應(yīng)的函數(shù)值。

2 彎曲板單元插值

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

首先假設(shè)氣動點(diǎn)上的壓力分布類似于有限元理論的彎曲薄板的位移分布。將氣動點(diǎn)作為虛擬的有限元節(jié)點(diǎn)，將氣動網(wǎng)格作為虛擬的有限元網(wǎng)格，虛擬節(jié)點(diǎn)上的壓力值作為彎曲板單元的節(jié)點(diǎn)的位移值，真實(shí)的有限元節(jié)點(diǎn)是虛擬單元的內(nèi)部點(diǎn)。根據(jù)有限元形函數(shù)插值理論，可由節(jié)點(diǎn)的位移值計(jì)算出單元內(nèi)的任一點(diǎn)的位移值，此位移值即為真實(shí)有限元節(jié)點(diǎn)上的壓力值。本文中采用的有限元單元是彎曲薄板單元中含有9個(gè)位移參數(shù)的部分協(xié)調(diào)三角形單元[9]。

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

式中：A是三角形的面積，A023、A031、A012是子三角形的面積。面積坐標(biāo)只有兩個(gè)是獨(dú)立的，因?yàn)樗鼈円獫M足關(guān)系：

三角形的面積為：

式（6）可展開為：

式中：w是單元內(nèi)部任意一點(diǎn)的位移；δe是單元節(jié)點(diǎn)位移；N是形函數(shù)。

形函數(shù)表達(dá)式為：

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

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

圖1 真實(shí)有限元節(jié)點(diǎn)與虛擬單元的位置關(guān)系Fig.1 The relation between true FEM node and virtual element

3 插值方法結(jié)果分析

3.1 壓力插值結(jié)果比較

按照以上介紹的兩種方法，以某飛機(jī)機(jī)翼主翼面上、下翼面某工況的載荷為例。將原始?xì)鈩狱c(diǎn)的壓力值和插值后有限元節(jié)點(diǎn)的壓力值進(jìn)行比較，如圖2～圖5所示。這里的壓力指的是當(dāng)?shù)仂o壓減去來流靜壓后的剩余量，即壓力系數(shù)和動壓的乘積。

圖2 上翼面薄板樣條插值結(jié)果與原始?xì)鈩狱c(diǎn)壓力比較Fig.2 Contrast between the thin plate spline result and the original aerodynamic press at up surface

圖3 上翼面彎曲薄板插值結(jié)果與原始?xì)鈩狱c(diǎn)壓力比較Fig.3 Contrast between bending element spline result and the original aerodynamic press at up surface

圖4 下翼面薄板樣條插值結(jié)果與原始?xì)鈩狱c(diǎn)壓力比較Fig.4 Contrast between the thin plate spline result and the original aerodynamic press at down surface

圖5 下翼面彎曲薄板插值結(jié)果與原始?xì)鈩狱c(diǎn)壓力比較Fig.5 Contrast between bending element spline result and the original aerodynamic press at down surface

3.2 翼面總載比較

計(jì)算得到有限元節(jié)點(diǎn)的壓力以后，等效節(jié)點(diǎn)載荷的計(jì)算比較簡單。設(shè)單元受法向分布載荷q的作用，對應(yīng)的等效節(jié)點(diǎn)載荷列矢量為：

如果分布載荷在單元內(nèi)是線性變化的，即：

式中：qi，qj和qm分別為q在節(jié)點(diǎn)i，j和m處的大小。

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

根據(jù)式（14）即可求得真實(shí)有限元模型上的等效氣動載荷，供結(jié)構(gòu)強(qiáng)度專業(yè)進(jìn)行計(jì)算。將該機(jī)翼此工況的上下翼面計(jì)算得到的節(jié)點(diǎn)力和力矩求和得到總載，并將結(jié)果與試驗(yàn)實(shí)測的總載比較，結(jié)果見表1。

表1 兩種計(jì)算結(jié)果總載荷與試驗(yàn)值比較Table1 Contrast between the total load of two method and test result

可以看出，兩種插值結(jié)果與原始壓力均符合的較好。相比較而言，彎曲薄板的結(jié)果更為光滑，奇點(diǎn)更少一些。表1中的總載對比也證明了彎曲薄板的計(jì)算結(jié)果總載符合的更好一些。此外，彎曲薄板插值方法的計(jì)算量也小一些，因?yàn)楸“鍢訔l插值的計(jì)算過程中需要求解大型線性代數(shù)方程組。

3.3 兩種方法結(jié)果對比分析

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

4 結(jié)束語

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