瞬態熱應力問題的組合雜交有限元方法

張玲，聶玉峰

(西北工業大學 理學院，西安 710129)

瞬態熱應力問題的組合雜交有限元方法

張玲，聶玉峰

(西北工業大學 理學院，西安 710129)

求解瞬態熱應力問題的雜交元方法面臨Ladyshenskaya-Babuska-Brezzi(LBB)條件的檢查，其穩定性難以保證，從而提出瞬態熱應力問題的組合雜交有限元方法。建立瞬態熱應力問題的基于區域分解的組合雜交變分原理及其有限元離散；通過數值算例驗證求解瞬態熱應力問題的組合雜交元的數值性能。結果表明：在大寬厚比網格下，八節點六面體組合雜交元(CHH(0-1))可以獲得與二十節點六面體元(B20)精度相當的位移計算結果和更好的應力計算結果。

瞬態；熱應力問題；組合變分原理；穩定性；組合雜交元；大寬厚比；應力精度

0 引 言

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

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

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

本文首先用傳統有限元方法求解瞬態熱傳導問題，得到初應變之后，發展熱應力問題的組合雜交有限元方法，以提高數值求解的精度及效率。

1 瞬態熱傳導問題求解

三維瞬態熱傳導問題如下：

(1)

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

(2)

對上式進行分部積分后，可得

(3)

(4)

對于單元內部的溫度，采用單元節點的溫度通過三線性插值獲得，即

(5)

則

(6)

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

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

(7)

2 熱應力問題求解

當物體溫度發生變化時，物體由于熱變形會產生線應變，如果各部分的熱變形不受約束，將不會產生應力。但如果受到的約束或者各部分溫度變化不均，即物體熱變形不能自由進行，則產生應力，把這部分應力稱為熱應力。物體由于熱膨脹只產生線應變，而剪切應變為0，這種由于熱變形產生的應變可以看作物體的初應變[15]。

步驟1 建立熱應力問題的組合雜交變分原理

由于存在初應變，與彈性力學問題相比，應力應變關系發生改變，熱應力問題的模型為

(8)

(9)

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

(10)

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

(11)

(12)

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

為了放松場函數的連續性，引入分片定義的Sobolev空間。

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

(13)

(14)

令

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

(15)

(16)

根據v和τ的任意性，式(15)和式(16)的等價表達形式分別為

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

(17)

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

(18)

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

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

(19)

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

(21)

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

(22)

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

選取線性插值空間為應力插值函數空間，Wilson bubble增強的三線性插值空間作為位移插值函數空間。

(23)

(24)

(25)

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

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

(26)

同時，簡化的誤差估計為

(27)

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

(28)

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

3 數值實驗

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

算例1 對稱板

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

當t=100 s時，分析大寬厚比網格下單元H8、B8、B20和CHH(0-1)的位移和應力的計算精度。

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

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

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

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

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

表3 點F(0,40,2)的應力σx

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

表4 大寬厚比網格下點F(0,40,2)的應力σx

算例2 非對稱板

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

當t=100 s時，分析大寬厚比網格下單元H8、B8、B20和CHH(0-1)的位移和應力的計算精度。

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

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

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

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

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

表7 點A(0.25,0,0)的應力σx

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

表8 大寬厚比網格下點A(0.25,0,0)的應力σx

綜上所述，對于瞬態熱應力問題，組合雜交元在大寬厚比網格下可以獲得較高的位移和應力計算精度。

4 結 論

本文首先建立瞬態熱應力問題的組合雜交變分原理，進而給出其組合雜交有限元方法。組合雜交有限元方法自動滿足LBB條件，是穩定的有限元方法。

對于瞬態熱應力問題，高性能八節點六面體組合雜交元在大寬厚比網格下可以獲得較高的位移和應力計算精度。與傳統的有限元方法相比，組合雜交有限元方法可以極大地節約計算量。

[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] 崔健. 熱應力分析中的應力雜交有限元方法研究[D]. 北京: 華北電力大學, 2002.

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

[4] 關春龍, 李壵, 赫曉東. 可重復使用熱防護系統防熱結構及材料的研究現狀[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] 周天孝. 彈塑性力學的典則變分原理[J]. 數值計算與計算機應用, 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]. 北京: 清華大學出版社, 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]. 西安: 中國航空計算技術研究所, 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-)，女，博士研究生。主要研究方向：熱力學問題的高效有限元方法。聶玉峰(1968-)，男，博士，教授。主要研究方向：科學工程計算的模型、理論與方法。

(編輯：趙毓梅)

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

國家自然科學基金(11471262，11501450)西北工業大學博士論文創新基金(CX201524)

聶玉峰,yfnie@nwpu.edu.cn

1674-8190(2017)04-367-08

O242.21

A

10.16615/j.cnki.1674-8190.2017.04.001