瞬態(tài)熱應(yīng)力問題的組合雜交有限元方法

張玲，聶玉峰

(西北工業(yè)大學(xué) 理學(xué)院，西安 710129)

瞬態(tài)熱應(yīng)力問題的組合雜交有限元方法

張玲，聶玉峰

(西北工業(yè)大學(xué) 理學(xué)院，西安 710129)

求解瞬態(tài)熱應(yīng)力問題的雜交元方法面臨Ladyshenskaya-Babuska-Brezzi(LBB)條件的檢查，其穩(wěn)定性難以保證，從而提出瞬態(tài)熱應(yīng)力問題的組合雜交有限元方法。建立瞬態(tài)熱應(yīng)力問題的基于區(qū)域分解的組合雜交變分原理及其有限元離散；通過數(shù)值算例驗(yàn)證求解瞬態(tài)熱應(yīng)力問題的組合雜交元的數(shù)值性能。結(jié)果表明：在大寬厚比網(wǎng)格下，八節(jié)點(diǎn)六面體組合雜交元(CHH(0-1))可以獲得與二十節(jié)點(diǎn)六面體元(B20)精度相當(dāng)?shù)奈灰朴?jì)算結(jié)果和更好的應(yīng)力計(jì)算結(jié)果。

瞬態(tài)；熱應(yīng)力問題；組合變分原理；穩(wěn)定性；組合雜交元；大寬厚比；應(yīng)力精度

0 引 言

在我國的很多國民經(jīng)濟(jì)部門，例如電力、航空航天、冶金以及化工等都廣泛使用動(dòng)力機(jī)械或者熱力發(fā)動(dòng)設(shè)備[1-2]，這些設(shè)備多數(shù)在高溫高壓的環(huán)境下進(jìn)行作業(yè)，設(shè)備在啟動(dòng)、停止和運(yùn)行過程中的溫度變化都會(huì)引起熱應(yīng)力[3]。特別地，空天飛行器在爬升和返回時(shí)，會(huì)經(jīng)受嚴(yán)峻的氣動(dòng)加熱，因此空天飛行器熱防護(hù)系統(tǒng)的材料選取、結(jié)構(gòu)設(shè)計(jì)和制造方法的研究涉及到力學(xué)環(huán)境以及熱學(xué)環(huán)境[4]。只有掌握溫度變化時(shí)所產(chǎn)生的熱應(yīng)力的分布規(guī)律，才能保證空天飛行器的安全，因此發(fā)展熱應(yīng)力問題的有效求解算法具有重要的意義。

近年來，模擬熱力耦合問題有多種方法，例如：有限元方法[5]、無網(wǎng)格方法[6]、多尺度分析方法[7-8]等。與有限元方法相比，無網(wǎng)格方法可以部分或徹底地?cái)[脫網(wǎng)格間的拓?fù)溥B接關(guān)系，但僅適用于處理大梯度和大變形問題，多尺度分析方法則可以以較少的計(jì)算量捕捉到材料的細(xì)觀尺度行為，但其均勻化方程的求解仍需依靠有限元方法。然而，傳統(tǒng)有限元方法存在各種閉鎖和對(duì)網(wǎng)格畸變敏感等現(xiàn)象，雜交有限元方法不能克服LBB條件這一困難[9]。

基于區(qū)域分解的Hellinger-Reissner變分原理及其對(duì)偶的變分原理，可以建立一種新的雜交元方法：組合雜交有限元方法[10]。對(duì)于彈性力學(xué)問題，它具有粗網(wǎng)格高精度、對(duì)網(wǎng)格畸變不敏感、不發(fā)生locking現(xiàn)象、避免LBB條件檢查的特性。從二維到三維，從梁問題到板問題，組合雜交有限元方法在彈性力學(xué)問題求解中顯示了優(yōu)越性[11-12]，近年來，它被發(fā)展到彈塑性力學(xué)問題以及穩(wěn)態(tài)熱應(yīng)力問題的求解方面[13-14]。然而，對(duì)于不同的物理問題，組合雜交有限元方法的變分原理不同，有限元離散后，剛度矩陣及右端項(xiàng)不同。

本文首先用傳統(tǒng)有限元方法求解瞬態(tài)熱傳導(dǎo)問題，得到初應(yīng)變之后，發(fā)展熱應(yīng)力問題的組合雜交有限元方法，以提高數(shù)值求解的精度及效率。

1 瞬態(tài)熱傳導(dǎo)問題求解

三維瞬態(tài)熱傳導(dǎo)問題如下：

(1)

式(1)的等效積分弱形式為

(2)

對(duì)上式進(jìn)行分部積分后，可得

(3)

(4)

對(duì)于單元內(nèi)部的溫度，采用單元節(jié)點(diǎn)的溫度通過三線性插值獲得，即

(5)

則

(6)

式中：Ni(i=1,2,…,8)為三線性插值基函數(shù)；Φi(i=1,2,…,8)為單元節(jié)點(diǎn)溫度值。

將式(5)和(6)代入式(4)，由于δΦi的任意性，消去該項(xiàng)后，在求解域Ω上，得到瞬態(tài)熱傳導(dǎo)問題的有限元方程。

(7)

2 熱應(yīng)力問題求解

當(dāng)物體溫度發(fā)生變化時(shí)，物體由于熱變形會(huì)產(chǎn)生線應(yīng)變，如果各部分的熱變形不受約束，將不會(huì)產(chǎn)生應(yīng)力。但如果受到的約束或者各部分溫度變化不均，即物體熱變形不能自由進(jìn)行，則產(chǎn)生應(yīng)力，把這部分應(yīng)力稱為熱應(yīng)力。物體由于熱膨脹只產(chǎn)生線應(yīng)變，而剪切應(yīng)變?yōu)?，這種由于熱變形產(chǎn)生的應(yīng)變可以看作物體的初應(yīng)變[15]。

步驟1 建立熱應(yīng)力問題的組合雜交變分原理

由于存在初應(yīng)變，與彈性力學(xué)問題相比，應(yīng)力應(yīng)變關(guān)系發(fā)生改變，熱應(yīng)力問題的模型為

(8)

(9)

ε0=β(Φ-Φ0)[1 1 1 0 0 0]T

(10)

式中：β為線膨脹系數(shù)；Φ0為初始溫度。

(11)

(12)

然后，建立熱應(yīng)力問題的組合雜交變分原理。

為了放松場(chǎng)函數(shù)的連續(xù)性，引入分片定義的Sobolev空間。

通過拉格朗日乘子法建立基于區(qū)域分解的變分原理。

(13)

(14)

令

對(duì)式(13)～式(14)求變分，得

(15)

(16)

根據(jù)v和τ的任意性，式(15)和式(16)的等價(jià)表達(dá)形式分別為

求(σ,u,uc)∈?！罸×Uc，使得

(17)

求(σ,u,uc)∈Γ×U×Uc，使得

(18)

根據(jù)文獻(xiàn)[16]，直接離散鞍點(diǎn)問題(17)或問題(18)均需要滿足LBB條件，為了避免LBB條件的檢查，引入權(quán)系數(shù)α∈(0,1)，建立組合雜交變分原理。

求(σ,u,uc)∈?！罸×Uc，使得

(19)

求(uh,σh)∈Uh×Γh，使得

(21)

同理，問題(20)～問題(21)解的適定性可以通過文獻(xiàn)[17]得到證明，并根據(jù)彈性力學(xué)問題誤差估計(jì)的證明思路，建立如下誤差估計(jì)：

(22)

步驟2 建立高性能組合雜交元

選取線性插值空間為應(yīng)力插值函數(shù)空間，Wilson bubble增強(qiáng)的三線性插值空間作為位移插值函數(shù)空間。

(23)

(24)

(25)

令b1[τ,v-Tc(v)]=0，簡化的變分方程如下：

求(uh,σh)∈Uh×Γh，使得

(26)

同時(shí)，簡化的誤差估計(jì)為

(27)

此時(shí)，誤差估計(jì)階數(shù)提高為O(ht):1lt;tlt;2，并且當(dāng)h趨于0，得到locking-free收斂性

(28)

令b1[τ,v-Tc(v)]=0，則長方體剖分情況下，應(yīng)力的顯式離散表達(dá)式見文獻(xiàn)[17]。

3 數(shù)值實(shí)驗(yàn)

通過數(shù)值算例來驗(yàn)證求解熱應(yīng)力問題的組合雜交元的數(shù)值性能。實(shí)驗(yàn)中采用的單元為：組合雜交元(CHH(0-1))、三線性元(H8)、Wilson元(B8)和二次元(B20)，三線性元、Wilson元和二次元的計(jì)算結(jié)果來自軟件Ansys。

算例1 對(duì)稱板

對(duì)稱載荷板(由于板的對(duì)稱性，取其八分之一進(jìn)行計(jì)算)如圖1所示，其長、寬、高分別是120、80、4 m。密度為2.5×10-3kg/m3，比熱為1.146 J/(kg·℃)，上下表面的對(duì)流換熱系數(shù)為1.8×10-4W/(m2·℃)，4個(gè)側(cè)面邊界為9.0×10-5W/(m2·℃)。板的彈性模量為7.9×103Pa，泊松比為0.17，熱傳導(dǎo)率kx=ky=kz=2.46×10-3W/(m·℃)，線膨脹系數(shù)β=5.04×10-7/℃，α=0.5。現(xiàn)將板加熱到600 ℃并突然放置到20 ℃的空氣中進(jìn)行淬冷，忽略熱傳導(dǎo)和熱輻射，將傳熱過程簡化為對(duì)流傳熱，求板淬冷過程中的熱應(yīng)力變化情況。

當(dāng)t=100 s時(shí)，分析大寬厚比網(wǎng)格下單元H8、B8、B20和CHH(0-1)的位移和應(yīng)力的計(jì)算精度。

B8單元在不同網(wǎng)格剖分下點(diǎn)A(60,0,0)的位移計(jì)算結(jié)果UX如表1所示，可以看出：點(diǎn)A位移的收斂值為0.597×10-2。

表1 點(diǎn)A(60,0,0)的位移UX

單元H8、B8、B20和CHH(0-1)在大寬厚比網(wǎng)格下點(diǎn)A(60,0,0)的位移計(jì)算結(jié)果UX如表2所示，可以看出：四種單元均具有較高的位移計(jì)算精度。

表2 大寬厚比網(wǎng)格下點(diǎn)A(60,0,0)的位移UX

B8單元在不同網(wǎng)格剖分下點(diǎn)F(0,40,2)的應(yīng)力計(jì)算結(jié)果σx如表3所示，可以看出：點(diǎn)F應(yīng)力σx的收斂值為1.307 3。

表3 點(diǎn)F(0,40,2)的應(yīng)力σx

單元H8、B8、B20和CHH(0-1)在大寬厚比網(wǎng)格下點(diǎn)F(0,40,2)的應(yīng)力計(jì)算結(jié)果σx如表4所示，可以看出：在10∶1大寬厚比網(wǎng)格剖分下，單元B20和CHH(0-1)的應(yīng)力精度高于單元H8和B8的應(yīng)力精度，但單元B20剛度矩陣的計(jì)算量是單元CHH(0-1)的5倍，線性方程組階數(shù)是單元CHH(0-1)的4倍[14]。

表4 大寬厚比網(wǎng)格下點(diǎn)F(0,40,2)的應(yīng)力σx

算例2 非對(duì)稱板

非對(duì)稱載荷板如圖2所示[14]，長、寬、高分別是6、1、0.1 m。密度為2.5×10-3，比熱為1.146 J/(kg·℃)，上下表面的對(duì)流換熱系數(shù)為1.8×10-4W/(m2·℃)，4個(gè)側(cè)面邊界為9.0×10-5W/(m2·℃)。板的彈性模量為1.0×107Pa，泊松比為0.3，熱傳導(dǎo)率kx=ky=kz=2.46×10-3W/(m·℃)，線膨脹系數(shù)β=5.04×10-1℃，α=0.5。板右端給定剪切載荷P=50 N?，F(xiàn)將板加熱到600 ℃并突然放置到20 ℃的空氣中進(jìn)行淬冷，忽略熱傳導(dǎo)和熱輻射，將傳熱過程簡化為對(duì)流傳熱，求板淬冷過程中的熱應(yīng)力變化情況。

當(dāng)t=100 s時(shí)，分析大寬厚比網(wǎng)格下單元H8、B8、B20和CHH(0-1)的位移和應(yīng)力的計(jì)算精度。

B8單元在不同網(wǎng)格剖分下點(diǎn)B(6,0,0.1)的位移計(jì)算結(jié)果UZ如表5所示，可以看出：點(diǎn)B位移UZ的收斂值為1.036 6。

表5 點(diǎn)B(6,0,0.1)的位移UZ

單元H8、B8、B20和CHH(0-1)在大寬厚比網(wǎng)格下點(diǎn)B(6,0,0.1)的位移計(jì)算結(jié)果UZ如表6所示，可以看出：單元B8、B20和CHH(0-1)在大寬厚比網(wǎng)格下均具有較高的位移計(jì)算精度，但單元H8位移精度較差。

表6 大寬厚比網(wǎng)格下點(diǎn)B(6,0,0.1)的位移UZ

B8單元在不同網(wǎng)格剖分下點(diǎn)A(0.25,0,0)的應(yīng)力計(jì)算結(jié)果σx如表7所示，可以看出：點(diǎn)A應(yīng)力σx的收斂值為0.308 6×108。

表7 點(diǎn)A(0.25,0,0)的應(yīng)力σx

單元H8、B8、B20和CHH(0-1)在大寬厚比網(wǎng)格下點(diǎn)A(0.25,0,0)的應(yīng)力計(jì)算結(jié)果σx如表8所示，可以看出：在10∶1大寬厚比網(wǎng)格剖分下，單元CHH(0-1)應(yīng)力計(jì)算精度高于單元H8、B8、B20的應(yīng)力計(jì)算精度。

表8 大寬厚比網(wǎng)格下點(diǎn)A(0.25,0,0)的應(yīng)力σx

綜上所述，對(duì)于瞬態(tài)熱應(yīng)力問題，組合雜交元在大寬厚比網(wǎng)格下可以獲得較高的位移和應(yīng)力計(jì)算精度。

4 結(jié) 論

本文首先建立瞬態(tài)熱應(yīng)力問題的組合雜交變分原理，進(jìn)而給出其組合雜交有限元方法。組合雜交有限元方法自動(dòng)滿足LBB條件，是穩(wěn)定的有限元方法。

對(duì)于瞬態(tài)熱應(yīng)力問題，高性能八節(jié)點(diǎn)六面體組合雜交元在大寬厚比網(wǎng)格下可以獲得較高的位移和應(yīng)力計(jì)算精度。與傳統(tǒng)的有限元方法相比，組合雜交有限元方法可以極大地節(jié)約計(jì)算量。

[1] Liu W, Qin Y P. Multi-physics coupling model of coal spontaneous combustion in longwall gob area based on moving coordinates[J]. Fuel, 2017, 188: 553-566.

[2] Su H, Rahmani R, Rahnejat H. Thermohydrodynamics of bidirectional groove dry gas seals with slip flow[J]. International Journal of Thermal Sciences, 2016, 110: 270-284.

[3] 崔健. 熱應(yīng)力分析中的應(yīng)力雜交有限元方法研究[D]. 北京: 華北電力大學(xué), 2002.

Cui Jian. The thermal stress by hybrid finite element method[D]. Beijing: North China Electric Power University, 2002.(in Chinese)

[4] 關(guān)春龍, 李壵, 赫曉東. 可重復(fù)使用熱防護(hù)系統(tǒng)防熱結(jié)構(gòu)及材料的研究現(xiàn)狀[J]. 宇航材料工藝, 2003, 33(6): 7-11, 42.

Guan Chunlong, Li Yao, He Xiaodong. Research status of structures and materials for reusable TPS[J]. Areospace Materials amp; Technology, 2003, 33(6): 7-11, 42.(in Chinese)

[5] Faruk S, Metin S. Elasto-plastic thermal stress analysis in a thermoplastic composite disc under uniform temperature using FEM[J]. Computers amp; Mathematics with Applications, 2006, 11(1): 31-39.

[6] Ren B, Qian J, Zeng X, et al. Recent developments on thermo-mechanical simulations of ductile failure by meshfree method[J]. Computer Modeling in Engineering amp; Sciences, 2011, 71(3): 253-277.

[7] Yang Z Q, Cui J Z, Sun Y, et al. Multiscale analysis method for thermo-mechanical performance of periodic porous materials with interior surface radiation[J]. International Journal for Numerical Methods in Engineering, 2016, 105(5): 323-350.

[8] Yang Z Q, Cui J Z, Zhou S. Thermo-mechanical analysis of periodic porous materials with microscale heat transfer by multiscale asymptotic expansion method[J]. International Journal of Heat and Mass Transfer, 2016, 92: 904-919.

[9] Brezzi F, Fortin M. Mixed and hybrid finite element method[M]. German: Springer-Verlag, 1991.

[10] Zhou T X. Stabilized hybrid finite element methods based on the combination of saddle point principles of elasticity problems[J]. Mathematics of Computation, 2003, 72(244): 1655-1673.

[11] Zhou T X, Xie X P. Zero energy-error mechanism of the combined hybrid method and improvement of Allman’s membrane element with drilling d.o.f.’s[J]. Communications in Numerical Methods in Engineering, 2004, 20(3): 241-250.

[12] Xie X P, Hu J C. Energy-adjustable mechanism of the combined hybrid finite element method and improvement of Zienkiewicz’s plate-element[J]. Communications in Numerical Methods in Engineering, 2005, 21(10): 531-544.

[13] 周天孝. 彈塑性力學(xué)的典則變分原理[J]. 數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用, 2011, 32(1): 1-7.

Zhou Tianxiao. Canonical variational principles in elasticity and plasticity[J]. Journal on Numerical Methods and Computer Applications, 2011, 32(1): 1-7.(in Chinese)

[14] Zhang L, Nie Y F, Yuan Z B, et al. Combined hybrid finite element method applied in elastic thermal stress problem[J]. International Journal of Computational Methods, 2017, 14(3): 1750071.

[15] 王勖成. 有限單元法[M]. 北京: 清華大學(xué)出版社, 2008.

Wang Xucheng. Finite element method[M]. Beijing: Tsinghua University Press, 2008.(in Chinese)

[16] Brezzi F, Fortin M. Mixed and hybrid finite element method[M]. German: Springer-Verlag, 1991.

[17] 聶玉峰. 高性能組合雜交有限元方法[D]. 西安: 中國航空計(jì)算技術(shù)研究所, 2000.

Nie Yufeng. Combined hybrid approach to finite element schemes of high performance[D]. Xi’an: Aeronaut Computing Technical Research Institute, 2000. (in Chinese)

張玲(1987-)，女，博士研究生。主要研究方向：熱力學(xué)問題的高效有限元方法。聶玉峰(1968-)，男，博士，教授。主要研究方向：科學(xué)工程計(jì)算的模型、理論與方法。

(編輯：趙毓梅)

CombinedHybridFiniteElementMethodAppliedinTransientThermalStressProblem

Zhang Ling, Nie Yufeng

(School of Natural and Applied Sciences, Northwestern Polytechnical University, Xi’an 710129, China)

The stability of hybrid finite element method of transient thermal stress problem is hard to be guaranteed because of the difficulty of satisfying the ladyshenskaya-babuska-brezzi(LBB) conditions, accordingly, the combined hybrid finite element method of transient thermal stress problem is proposed. Combinative variation principle and the corresponding finite element discretization of transient thermal stress problem is constructed on the basis of domain decomposition technique. The numerical performance of combined hybrid element for solving transient thermal stress problem is verified by numerical experiments. The numerical results indicate that combined hybrid element with 8 nodes(CHH(0-1)) can give almost the same computing accuracy of displacement and better computing accuracy of stress compared with cuboid element with 20 nodes(B20).

transient; thermal stress problem; combinative variation principle; stability; combined hybrid element; big width-to-thickness ratio; stress accuracy

2017-09-26；

2017-10-12

國家自然科學(xué)基金(11471262，11501450)西北工業(yè)大學(xué)博士論文創(chuàng)新基金(CX201524)

聶玉峰,yfnie@nwpu.edu.cn

1674-8190(2017)04-367-08

O242.21

A

10.16615/j.cnki.1674-8190.2017.04.001