拉格朗日插值方法和等參單元逆變換方法在熱固耦合中的應用

金賽英，陳鐵鋒

軟件應用

拉格朗日插值方法和等參單元逆變換方法在熱固耦合中的應用

金賽英，陳鐵鋒

（中國航發商用航空發動機有限責任公司，上海200241）

基于拉格朗日插值方法，結合等參單元逆變換方法，將目標插值節點的全局坐標變化成等參單元的局部坐標，獲得拉格朗日插值所需的形函數，從而得到目標插值節點所對應的熱固耦合結果。對該方法進行了驗證并程序化，適用于航空發動機二維整機有限元分析模型進行了熱固耦合，算例表明該方法具有較高的精度，工程適用性好，能用于類似的多種航空發動機結構進行熱固耦合。

拉格朗日插值方法；等參單元逆變換；航空發動機；熱固耦合

在工程設計上一般需要熱固耦合加載溫度載荷，目前常用的熱固耦合方法采用就近插值方法或者商用軟件自帶的插值方法，在解決溫度梯度大、熱模型和結構模型不一致時存在工程誤差大的特點。針對航空發動復雜精細的結構特征，常規的熱固耦合方法，無法滿足篦齒和輪筒等截面變化梯度比較大的復雜位置的溫度載荷精度要求，也無法滿足獲得精細結構輪廓的冷熱態換算計算的溫度載荷精度要求。目前國內在耦合方法方面開展了大量的研究工作，但關于熱固耦合的兩種物理結構網格上的參數轉換問題，主要針對插值方法以及單元間結構參數的有效變換等問題開展了研究[1-5]，Mutri V and Valliappan S和錢向東等人提出了基于泰勒展開的等參單元逆變換方法[6-8]，朱以文等開展了比較精細的等參單元換算方法研究[9-10]，其提出的等參有限元逆變換高效算法工程適用性非常好。本文對某型號二維整機有限元模型的溫度場進行了研究，基于拉格朗日-權函數插值方法[11]，采用等參單元逆變換思想，將目標插值節點的全局坐標變成等參單元的局部坐標，從而獲得插值形函數，最終得到目標插值節點的溫度載荷。插值結果表明基于拉格朗日-權函數插值方法和采用等參單元逆變換方法得到的轉子溫度分布準確度高，工程適用性好。航空發動機大部分結構都需要承受溫度載荷，本方法能廣泛應用于轉子、機匣、軸類等航空發動機結構，并已形成了相應的工程軟件。

1 拉格朗日插值方法及形函數

1.1 拉格朗日插值方法

采用拉格朗日形函數插值方法，需要考慮單元類型以及階次的組合以及不同單元類型和階次組合下的形函數類別。對于熱固耦合問題，基于熱模型的單元、節點布局，獲得所有結構模型節點的相應形參，從而得到相應所有目標節點的插值載荷。如圖1圓塊為熱模型，方塊為結構模型，結構模型上的一個待插值節點p處在熱模型的單元ijk內。如公式（1）所示Li、Lj、Lk為目標節點p在三角形單元△ijk內對應的形函數，LOAD（i）、LOAD（j）、LOAD（k）為相應的溫度模型單元節點的溫度。

圖1 熱固耦合模型

形函數是定義于單元內部的、坐標的連續函數，它應滿足一下條件：

a）在節點i，形函數Ni=1；在其他節點，Ni=0；

b）能保證用它定義的未知量（u，v或x，y）在相鄰單元之間的連續性；

c）它包含任意線性項，以便用它定義的單元位移可滿足常應變條件；

d）應滿足下列等式∑Ni.

因此（1）式中，Li+Lj+Lk=1.

1.2 單元形函數

常用的二維熱模型和二維結構模型的離散有限元網格有一次三角形單元、二次三角形單元、一次四邊形單元和二次四邊形單元。如圖2所示為四種單元類型的母單元。

圖2 二維母單元（從左至右分別為：一次三角形單元，二次三角形單元，一次四邊形單元，二次四邊形單元）

a）對于一次三角形單元（3節點），基于面積計算得到形函數

本文研究的熱固耦合模型基于二維結構，考慮單元類型所要用到形函數為以上四種。對于三角形單元，只要定位目標插值節點在熱模型上的精確位置，就可通過單元內的形函數積分獲得其插值結果。而對于四邊形單元，除了定位其在熱模型上的精確位置，還需獲得在等參單元對應的ξ和η值。通常的四邊形單元不是規則的母單元，需要等參變換為規則的母單元，每個母單元都具有一個獨立的局部坐標系，這里目標插值節點只有全局坐標系下的坐標值，基于拉格朗日形函數法，需要目標插值接節點所在熱模型單元對應的母單元的局部坐標系下的坐標值。因此要開展等參單元逆變換，獲得目標插值節點的局部坐標值，即ξ和η.

2 等參單元逆變換方法

如圖3所示為四邊形等參單元的全局坐標系在全局坐標系和局部坐標系下的示意圖，在四個節點的局部坐標分別為（-1，-1），（1，-1），（1，1），（-1，1），在局部坐標系下單元內任意位置的坐標為滿足

圖3 四邊形等參單元坐標模型

在全局坐標系下（x為橫坐標，y為縱坐標），單元內部任意一點的坐標可以用如（8）式表達：

式中，x0和y0為單元內任意位置在全局坐標系下的坐標值，ξ和η為對應的在局部坐標系下的坐標值。

本文所采用的等參單元逆變換，將結構模型的目標插值節點映射到熱模型相應的單元中，獲得相應的局部坐標值。如（8）式所示，已知xi和yi（i =0，1，2，3，4），求ξ和η值。理論上通過一組聯立方程，可以數值求解兩個未知量。但是之前的研究表明，這種數值求解存在較大的誤差，特別是在單元交界位置存在無解的情況。因此必須建立一種有效的等參單元逆變換的方法，在保證求解精度的同時能完成所有目標節點的變換。本文所采用的變換方法如下公式所示：

3 算法驗證及程序

3.1 拉格朗日算法驗證

為了驗證拉格朗日算法的準確性，建立了用于驗證的熱模型和有限元模型。如圖4所示，熱分析模型耦合了方塊狀有限元模型的網格，在熱分析過程中將有限元模型對應的網格節點溫度也獲得。輸出溫度時將有限元部分的網格溫度信息不輸出，只輸出如圖4所示的圓塊狀分布。根據如圖4所輸出的圓塊狀溫度，需要通過插值分析獲得如圖5所示的方塊狀有限元模型的溫度。

圖4 熱分析模型

圖5 有限元分析模型

從圖6可見，方塊狀的溫度分布在與熱模型疊加時的溫度情況，就近插值方法不能復原熱模型的溫度場分布規律，而本方法能達到與初始熱模型完全吻合的程度。

圖6 熱固耦合結果

通過對照有限元模型各個節點的溫度結果，如表1所示。在節點3、4、5、13、14、17、18、19等節點處的誤差最大，這些節點的共同特性是位于單元內部，離最近的熱模型節點較遠。采用拉格朗日插值方法能克服網格節點之間的局限性，通過形函數將目標節點與熱模型的有限元模型連接起來，作為一種同等的離散節點，從而實現目標節點的精確插值，本方法插值誤差最高為6.4%，絕大多數目標耦合點能達到零誤差。常用的就近插值方法有7個點的插值誤差超過10%，在溫度梯度大的點高達53.8%.

表1 熱固耦合插值結果

3.2 程序算法

航空發動機結構復雜，插值方法需要很強大的通用性，因此在處理算法的時候應該考慮到各個算法模塊化，同時熱固耦合的輸入輸出規范化。

具體的程序方法實現如圖7所示。

圖7 程序流程圖

4 算例

采用基于拉格朗日和等參單元逆變換方法及其完善的耦合程序應用于某型二維整機有限元模型溫度場加載，如圖8為高壓渦輪一級盤盤心位置熱固耦合情況，圖9所示為低壓渦輪一級盤輪緣及鼓筒位置熱固耦合情況，圖10為渦輪部件的整個部件熱固耦合結果。

通過對某型二維整機有限元模型的表明，本方法插值結果精確較高，程序穩定性好，可以用于解決整機的轉子、輪盤和軸類的熱固耦合問題。

圖8 高壓渦輪一級盤盤心拉格朗日插值熱固耦合

圖9 低壓渦輪一級盤輪緣及鼓筒處拉格朗日插值熱固耦合

圖10 高低壓渦輪轉靜子有限元熱固耦合結果

5 結束語

采用拉格朗日形函數思想，結合等參單元逆變換方法，實現了在熱模型離散單元內部插值，結論如下：

（1）本文針對熱固耦合的離散單元模型，將結構節點理解為熱模型離散單元連續體的某個點，通過采用拉格朗日形函數方法，將所有的結構節點歸一化到所對應的熱模型的某個單元內。同時，采用等參單元逆變換方法，將目標節點從全局坐標系變換到局部坐標系下，最終通過單元內插值，得到所有目標插值節點的熱耦合值。

（2）針對熱模型為一次以及二次的三角形單元和四邊形單元混合的離散模型，采用本方法形成了通用的工程軟件，并成功用于某型二維整機有限元模型的溫度場插值，插值結果非常精確，具有很強的工程適用性，可以廣泛應用于類似結構的強度、壽命評估，以及需要獲得精確結構冷、熱態結構的冷熱態換算。

（3）目前考慮的單元類型為二維平面內的單元，基于本文提出的方法，完善三維形函數庫和三維空間的等參單元逆變換，可用于空間殼體模型（如葉片表面單元）以及三維模型的熱固耦合問題，并擴展至空間流固耦合問題。

[1]劉振宇，傅云，譚建榮，基于異構網格耦合的產品多物理場有限元數據集成與可視化仿真[J].機械工程學報，2010，46（7）：114-121.

[2]方錫武，劉振宇，譚建榮，等.異構有限元網格多場信息的等效集成方法[J].浙江大學學報，2014，48（6）：301-302

[3]沈毅，三維曲面插值在發動機渦輪葉片有限元溫度場計算中的應用[J].航空發動機，2002（2）：43－45.

[4]包勁青，楊強，官福海.線性等參單元逆變換的二分法及其修正[J].計算力學學報，2010，27（5）：770-774.

[5]李春光，鄭宏，葛修潤，等.六面體單元等參逆變換的一種迭代解法[J].巖土力學，2004，25（7）：1050-1052.

[6]Mutri V and Valliappan S numerical inverse isoparametric mapping in remeshing and nodal quantity continuing[J].com puter&structure，1986，22：1011－1012.

[7]Mutri V，Wang Y and and Valliappan S numerical inverse isoparametric mapping in 3D FEM[J].computer&structure，1988，29：611－622.

[8]錢向東，任文青.一種高效的等參有限元逆變換算法[J].計算力學學報，1998（4）：437-441.

[9]朱以文，李偉，蔡元奇.基于解析性質的等參有限元逆變換高效算法[J].武漢大學學報，2002，35（2）：62-65.

[10]徐燕萍，項陽，劉泉聲，等.參逆變換插值法的研究及其應用[J].巖土力學，2001，22（2）：226-228.

[11]王勖成.有限單元法[M].北京：清華大學出版社，2003：468－472.

Lagrange Interpolation Method and Isoparametric Element Inverse Transformation Method Application in Thermosetting Coupling

JIN Sai－ying，CHEN Tie－feng
（Chinese Hangfa Commercial Aircraft Engine Co.，Ltd.，Shanghai 200241，China）

Based on the lagrangian interpolation method and inverse is oparametric element mapping method，to get the local coordinate of a target point in the aera of a isoparametric element by changing it’s global coordinate，then obtain the desired parameters of the lagrangian interpolation shape function，so the corresponding thermalstruture coupling result of the target mapping point be arrived.In this paper，the above method be validated and be programed，and applicable to thermal-structure coupling of the two-dimensional aero-engine rotor structure，the result show that this method has high accuracy and good adaptability in engineering utilize，it can be used in a variety of similar structure of aero-engine to thermal-structure coupling.

lagrangian interpolation method；inverse isoparametric element mapping；aero-engine thermal-structure coupling

TP301.6

B

1672-545X（2017）07-0232-05

2017-04-09

金賽英（1986-），女，江西吉安人，工學碩士，工程師，研究方向：航空發動機結構強度。