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

金賽英，陳鐵鋒

軟件應(yīng)用

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

金賽英，陳鐵鋒

（中國(guó)航發(fā)商用航空發(fā)動(dòng)機(jī)有限責(zé)任公司，上海200241）

基于拉格朗日插值方法，結(jié)合等參單元逆變換方法，將目標(biāo)插值節(jié)點(diǎn)的全局坐標(biāo)變化成等參單元的局部坐標(biāo)，獲得拉格朗日插值所需的形函數(shù)，從而得到目標(biāo)插值節(jié)點(diǎn)所對(duì)應(yīng)的熱固耦合結(jié)果。對(duì)該方法進(jìn)行了驗(yàn)證并程序化，適用于航空發(fā)動(dòng)機(jī)二維整機(jī)有限元分析模型進(jìn)行了熱固耦合，算例表明該方法具有較高的精度，工程適用性好，能用于類似的多種航空發(fā)動(dòng)機(jī)結(jié)構(gòu)進(jìn)行熱固耦合。

拉格朗日插值方法；等參單元逆變換；航空發(fā)動(dòng)機(jī)；熱固耦合

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

1 拉格朗日插值方法及形函數(shù)

1.1 拉格朗日插值方法

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

圖1 熱固耦合模型

形函數(shù)是定義于單元內(nèi)部的、坐標(biāo)的連續(xù)函數(shù)，它應(yīng)滿足一下條件：

a）在節(jié)點(diǎn)i，形函數(shù)Ni=1；在其他節(jié)點(diǎn)，Ni=0；

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

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

d）應(yīng)滿足下列等式∑Ni.

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

1.2 單元形函數(shù)

常用的二維熱模型和二維結(jié)構(gòu)模型的離散有限元網(wǎng)格有一次三角形單元、二次三角形單元、一次四邊形單元和二次四邊形單元。如圖2所示為四種單元類型的母單元。

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

a）對(duì)于一次三角形單元（3節(jié)點(diǎn)），基于面積計(jì)算得到形函數(shù)

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

2 等參單元逆變換方法

如圖3所示為四邊形等參單元的全局坐標(biāo)系在全局坐標(biāo)系和局部坐標(biāo)系下的示意圖，在四個(gè)節(jié)點(diǎn)的局部坐標(biāo)分別為（-1，-1），（1，-1），（1，1），（-1，1），在局部坐標(biāo)系下單元內(nèi)任意位置的坐標(biāo)為滿足

圖3 四邊形等參單元坐標(biāo)模型

在全局坐標(biāo)系下（x為橫坐標(biāo)，y為縱坐標(biāo)），單元內(nèi)部任意一點(diǎn)的坐標(biāo)可以用如（8）式表達(dá)：

式中，x0和y0為單元內(nèi)任意位置在全局坐標(biāo)系下的坐標(biāo)值，ξ和η為對(duì)應(yīng)的在局部坐標(biāo)系下的坐標(biāo)值。

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

3 算法驗(yàn)證及程序

3.1 拉格朗日算法驗(yàn)證

為了驗(yàn)證拉格朗日算法的準(zhǔn)確性，建立了用于驗(yàn)證的熱模型和有限元模型。如圖4所示，熱分析模型耦合了方塊狀有限元模型的網(wǎng)格，在熱分析過(guò)程中將有限元模型對(duì)應(yīng)的網(wǎng)格節(jié)點(diǎn)溫度也獲得。輸出溫度時(shí)將有限元部分的網(wǎng)格溫度信息不輸出，只輸出如圖4所示的圓塊狀分布。根據(jù)如圖4所輸出的圓塊狀溫度，需要通過(guò)插值分析獲得如圖5所示的方塊狀有限元模型的溫度。

圖4 熱分析模型

圖5 有限元分析模型

從圖6可見(jiàn)，方塊狀的溫度分布在與熱模型疊加時(shí)的溫度情況，就近插值方法不能復(fù)原熱模型的溫度場(chǎng)分布規(guī)律，而本方法能達(dá)到與初始熱模型完全吻合的程度。

圖6 熱固耦合結(jié)果

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

表1 熱固耦合插值結(jié)果

3.2 程序算法

航空發(fā)動(dòng)機(jī)結(jié)構(gòu)復(fù)雜，插值方法需要很強(qiáng)大的通用性，因此在處理算法的時(shí)候應(yīng)該考慮到各個(gè)算法模塊化，同時(shí)熱固耦合的輸入輸出規(guī)范化。

具體的程序方法實(shí)現(xiàn)如圖7所示。

圖7 程序流程圖

4 算例

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

通過(guò)對(duì)某型二維整機(jī)有限元模型的表明，本方法插值結(jié)果精確較高，程序穩(wěn)定性好，可以用于解決整機(jī)的轉(zhuǎn)子、輪盤和軸類的熱固耦合問(wèn)題。

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

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

圖10 高低壓渦輪轉(zhuǎn)靜子有限元熱固耦合結(jié)果

5 結(jié)束語(yǔ)

采用拉格朗日形函數(shù)思想，結(jié)合等參單元逆變換方法，實(shí)現(xiàn)了在熱模型離散單元內(nèi)部插值，結(jié)論如下：

（1）本文針對(duì)熱固耦合的離散單元模型，將結(jié)構(gòu)節(jié)點(diǎn)理解為熱模型離散單元連續(xù)體的某個(gè)點(diǎn)，通過(guò)采用拉格朗日形函數(shù)方法，將所有的結(jié)構(gòu)節(jié)點(diǎn)歸一化到所對(duì)應(yīng)的熱模型的某個(gè)單元內(nèi)。同時(shí)，采用等參單元逆變換方法，將目標(biāo)節(jié)點(diǎn)從全局坐標(biāo)系變換到局部坐標(biāo)系下，最終通過(guò)單元內(nèi)插值，得到所有目標(biāo)插值節(jié)點(diǎn)的熱耦合值。

（2）針對(duì)熱模型為一次以及二次的三角形單元和四邊形單元混合的離散模型，采用本方法形成了通用的工程軟件，并成功用于某型二維整機(jī)有限元模型的溫度場(chǎng)插值，插值結(jié)果非常精確，具有很強(qiáng)的工程適用性，可以廣泛應(yīng)用于類似結(jié)構(gòu)的強(qiáng)度、壽命評(píng)估，以及需要獲得精確結(jié)構(gòu)冷、熱態(tài)結(jié)構(gòu)的冷熱態(tài)換算。

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

[1]劉振宇，傅云，譚建榮，基于異構(gòu)網(wǎng)格耦合的產(chǎn)品多物理場(chǎng)有限元數(shù)據(jù)集成與可視化仿真[J].機(jī)械工程學(xué)報(bào)，2010，46（7）：114-121.

[2]方錫武，劉振宇，譚建榮，等.異構(gòu)有限元網(wǎng)格多場(chǎng)信息的等效集成方法[J].浙江大學(xué)學(xué)報(bào)，2014，48（6）：301-302

[3]沈毅，三維曲面插值在發(fā)動(dòng)機(jī)渦輪葉片有限元溫度場(chǎng)計(jì)算中的應(yīng)用[J].航空發(fā)動(dòng)機(jī)，2002（2）：43－45.

[4]包勁青，楊強(qiáng)，官福海.線性等參單元逆變換的二分法及其修正[J].計(jì)算力學(xué)學(xué)報(bào)，2010，27（5）：770-774.

[5]李春光，鄭宏，葛修潤(rùn)，等.六面體單元等參逆變換的一種迭代解法[J].巖土力學(xué)，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].計(jì)算力學(xué)學(xué)報(bào)，1998（4）：437-441.

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

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

[11]王勖成.有限單元法[M].北京：清華大學(xué)出版社，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-），女，江西吉安人，工學(xué)碩士，工程師，研究方向：航空發(fā)動(dòng)機(jī)結(jié)構(gòu)強(qiáng)度。