李慶華,陳莘莘,徐青

(華東交通大學(xué) 土木建筑學(xué)院， 南昌 330013)

?

三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題的自然單元法

李慶華,陳莘莘,徐青

(華東交通大學(xué) 土木建筑學(xué)院， 南昌 330013)

摘要：為了更有效地求解三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題，對(duì)無(wú)網(wǎng)格自然單元法應(yīng)用于此類問(wèn)題進(jìn)行了研究，并發(fā)展了相應(yīng)的計(jì)算方法?；趲缀涡螤詈瓦吔鐥l件的軸對(duì)稱性，三維的軸對(duì)稱問(wèn)題可降為二維平面問(wèn)題。為了簡(jiǎn)化本質(zhì)邊界條件的施加，軸對(duì)稱面上的溫度場(chǎng)采用自然鄰近插值進(jìn)行離散。功能梯度材料特性的變化由高斯點(diǎn)的材料參數(shù)進(jìn)行模擬。時(shí)間域上，采用傳統(tǒng)的兩點(diǎn)差分法進(jìn)行離散求解，進(jìn)而得到瞬態(tài)溫度場(chǎng)的響應(yīng)。數(shù)值算例結(jié)果表明，提出的方法是行之有效的,理論及方法不僅拓展了自然單元法的應(yīng)用范圍，而且對(duì)三維軸對(duì)稱瞬態(tài)熱傳導(dǎo)分析具有普遍意義。

關(guān)鍵詞：自然單元法；軸對(duì)稱；功能梯度材料；瞬態(tài)熱傳導(dǎo)

功能梯度材料是通過(guò)特定的材料制備工藝將不同性能的兩種或兩種以上材料按一定的設(shè)計(jì)規(guī)律組合起來(lái)的新型非均勻復(fù)合材料[1]。功能梯度材料的最大特點(diǎn)是材料參數(shù)的連續(xù)性，完全避免了層合復(fù)合材料的材料參數(shù)在層層之間的間斷面處不連續(xù)的問(wèn)題，提高了材料強(qiáng)度和耐熱性。因此，功能梯度材料在航空、航天及核反應(yīng)堆等高溫環(huán)境中具有廣泛的應(yīng)用潛力，對(duì)功能梯度材料的熱力學(xué)行為進(jìn)行研究十分必要[2-5]。然而，相對(duì)于二維平面問(wèn)題，目前對(duì)于三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題的數(shù)值方法研究相對(duì)較少[6-8]。

自然單元法[9-10]是一種新興的數(shù)值分析方法，因其獨(dú)特的優(yōu)勢(shì)，得到了國(guó)內(nèi)外許多學(xué)者的極大關(guān)注[11-15]。這種方法基于離散節(jié)點(diǎn)的Voronoi圖和Delaunay三角化幾何結(jié)構(gòu)，采用自然鄰近插值構(gòu)造全域近似函數(shù)和試函數(shù)。自然鄰近插值方案構(gòu)造簡(jiǎn)單，不涉及到復(fù)雜的矩陣求逆運(yùn)算，而且不需要任何人為參數(shù)，從而避免了無(wú)單元伽遼金法中由于不確定的影響半徑造成的影響域計(jì)算的不確定性。此外，自然單元法的形函數(shù)滿足插值性質(zhì)，可以準(zhǔn)確地施加本質(zhì)邊界條件，無(wú)需其他無(wú)網(wǎng)格法類似的特殊處理過(guò)程。自然單元法已經(jīng)被成功地應(yīng)用于很多領(lǐng)域，但目前尚未見(jiàn)到三維軸對(duì)稱熱傳導(dǎo)分析的無(wú)網(wǎng)格自然單元法的研究成果。

為了進(jìn)一步拓展自然單元法的應(yīng)用范圍，本文基于加權(quán)殘值法詳細(xì)推導(dǎo)了三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)分析的自然單元法理論公式，并給出了其詳細(xì)的數(shù)值實(shí)現(xiàn)過(guò)程。在此基礎(chǔ)上，采用FORTRAN自編了相關(guān)的計(jì)算程序。最后，通過(guò)典型算例的計(jì)算和對(duì)比分析，不僅驗(yàn)證了自然單元法應(yīng)用于三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)分析的有效性和合理性，并且討論了梯度參數(shù)的變化對(duì)計(jì)算結(jié)果的影響。

1自然鄰近插值

TI={x∈R2:d(x,xI)

(1)

式中:d(x,xI)為點(diǎn)x與節(jié)點(diǎn)xI的距離。

為確定Sibson插值形函數(shù)，定義二次Voronoi結(jié)構(gòu)為

(2)

在幾何意義上，TIJ實(shí)際上是以節(jié)點(diǎn)xI為最近點(diǎn)，節(jié)點(diǎn)xJ為次近點(diǎn)的空間點(diǎn)位置的集合。圖1所示為平面7個(gè)節(jié)點(diǎn)的Voronoi結(jié)構(gòu)和待插點(diǎn)x的二次Voronoi結(jié)構(gòu)。

構(gòu)造出插值點(diǎn)x的一次和二次Voronoi結(jié)構(gòu)后，插值點(diǎn)x的形函數(shù)及其導(dǎo)數(shù)可以寫(xiě)為[9]

(3)

(4)

式中:AI(x)表示插值點(diǎn)x的二次Voronoi結(jié)構(gòu)TxI的面積;A(x)表示插值點(diǎn)x的一次Voronoi結(jié)構(gòu)Tx的面積，即TxI的總和。

圖1　點(diǎn)x的一次和二次Voronoi結(jié)構(gòu)Fig.1　First-order and second-order Voronoi cells about

定義了各節(jié)點(diǎn)的插值函數(shù)后，點(diǎn)x的溫度函數(shù)類似于有限元法可寫(xiě)為

(5)

式中:TI(I=1,2,…,n)是計(jì)算點(diǎn)x周圍自然鄰節(jié)點(diǎn)I的節(jié)點(diǎn)溫度。

2控制方程的弱形式及其離散化

在三維軸對(duì)稱瞬態(tài)熱傳導(dǎo)問(wèn)題中，空間參量只有半徑r和軸向參量z。在軸對(duì)稱面Ω上，Γ1和Γ2分別表示問(wèn)題域的本質(zhì)邊界和自然邊界?？臻g的溫度分布應(yīng)該滿足下面的偏微分方程

(6)

式中：T(x,t)為瞬態(tài)溫度場(chǎng);t為時(shí)間，Q(x,t)為單位體積的內(nèi)部熱源;k(x)為熱傳導(dǎo)系數(shù)，α(x)為熱擴(kuò)散系數(shù)。對(duì)應(yīng)的邊界條件和初始條件為

(7a)

(7b)

(7c)

式(6)的等效積分形式的伽遼金提法可表示為

(8)

式中:δT代表真實(shí)溫度的變分。對(duì)式(8)進(jìn)行分部積分并考慮式(7b)，可得

(9)

由于只對(duì)空間域進(jìn)行離散，求解域Ω內(nèi)的試函數(shù)Th(x,t)可由式(5)表示為

(10)

將空間離散后的溫度表達(dá)式(10)代入式(9)，并考慮到δTI的任意性，就可以得到最終的線性常微分方程組

(11)

式中:

(12a)

(12b)

(12c)

3時(shí)間積分方案

經(jīng)過(guò)對(duì)空間域的無(wú)網(wǎng)格自然單元法離散，已將溫度場(chǎng)的偏微分方程轉(zhuǎn)化為一組以時(shí)間t為獨(dú)立變量的一階線性常微分方程組。本文中采用傳統(tǒng)的兩點(diǎn)差分法[16]對(duì)時(shí)間域進(jìn)行離散。

圖2　空心圓筒的節(jié)點(diǎn)布置方案Fig.2　Nodal distribution for a hollow

一般情況下，兩點(diǎn)差分格式可以寫(xiě)成

(13)

式中:Δt是時(shí)間步長(zhǎng)，θ為時(shí)間加權(quán)系數(shù)。將上式代入式(11)，經(jīng)過(guò)整理后可得到

(14)

式中

(15a)

(15b)

通過(guò)從0到1更改參數(shù)θ，有關(guān)這個(gè)系列方法都可以得到。本文取θ=2/3。

4數(shù)值算例

4.1空心圓筒

有一無(wú)限長(zhǎng)的功能梯度空心圓筒，其內(nèi)半徑a=8×10-2m，外半徑b=10×10-2m，如圖2所示。為了該問(wèn)題的分析，截取長(zhǎng)L=2×10-2m的一段，且上下兩截面視為絕熱面。內(nèi)表面溫度T0=0 ℃，外表面溫度T1=1 ℃，初始溫度給定為零。計(jì)算中，均勻布置45(9×5)個(gè)節(jié)點(diǎn)，時(shí)間步長(zhǎng)取Δt=0.5s。

假設(shè)熱傳導(dǎo)系數(shù)與熱擴(kuò)散系數(shù)沿徑向指數(shù)變化，即

(16)

式中:k0=17W/m· ℃，α0=1.7×10-5m2/s。為了檢驗(yàn)本文方法的有效性，取k=k0和α=α0計(jì)算得到了均質(zhì)材料(γ=0)的情況下空心圓筒在r=9cm處溫度隨時(shí)間的變化，如圖3所示。為了進(jìn)行對(duì)比，圖3中還給出了相應(yīng)的解析解[17]。可以看出，本文方法與解析解精確吻合，證明本文方法是有效的。

圖3　均質(zhì)空心圓筒在r=9.0 cm處溫度隨時(shí)間的變化Fig.3　Time variations of the temperature at r=9.0 cm in the hollow cylinder with homogeneous material

為了研究梯度參數(shù)γ對(duì)溫度變化的影響，采用自然單元法計(jì)算得到了γ=0.0，γ=0.2 cm-1和γ=0.5 cm-1時(shí)空心圓筒在r=9 cm處溫度隨時(shí)間的變化，如圖4所示。正如我們期待的，在圖4中，隨著梯度參數(shù)的增加，材料導(dǎo)熱性增加，對(duì)應(yīng)的溫度也隨之上升。

圖4　功能梯度空心圓筒在r=9.0 cm處溫度隨時(shí)間的變化Fig.4　Time variations of the temperature at r=9.0 cm in a functionally graded hollow

4.2實(shí)心圓筒

為了進(jìn)一步驗(yàn)證本文數(shù)值方法的有效性，考慮如圖5所示的一無(wú)限長(zhǎng)的功能梯度實(shí)心圓筒，其半徑a=1 m。截取長(zhǎng)L=1 m的一段，外表面溫度T=100 ℃，其余邊界絕熱，初始溫度給定為零。計(jì)算中，沿徑向均勻布置9個(gè)節(jié)點(diǎn)，沿軸向均勻布置5個(gè)節(jié)點(diǎn)，時(shí)間步長(zhǎng)取Δt=0.02 s。

圖5　實(shí)心圓筒的節(jié)點(diǎn)布置方案Fig.5　Nodal distribution for a full cylinder

假設(shè)熱傳導(dǎo)系數(shù)與熱擴(kuò)散系數(shù)沿徑向指數(shù)變化，即

(17)

式中:k0=1.0 W/m·℃，α0=1.0 m2/s。與前一算例類似，取k=k0和α=α0對(duì)均質(zhì)材料(γ=0)情況下的實(shí)心圓筒進(jìn)行了瞬態(tài)熱傳導(dǎo)分析。圖6給出了均質(zhì)實(shí)心圓筒在對(duì)稱軸處的溫度隨時(shí)間變化的本文數(shù)值解與解析解[17]的比較，可以看出采用自然單元法的數(shù)值解與解析解吻合很好，說(shuō)明自然單元法在三維軸對(duì)稱瞬態(tài)熱傳導(dǎo)分析中可以取得很高的計(jì)算精度。

圖6　均質(zhì)實(shí)心圓筒對(duì)稱軸處溫度隨時(shí)間的變化Fig.6　Time variations of the temperature on the axis of the full cylinder with homogeneous material

為了進(jìn)一步研究梯度參數(shù)γ對(duì)溫度變化的影響，圖7給出了γ=0.0，γ=-1 m-1和γ=-2 m-1時(shí)采用自然單元法計(jì)算得到的實(shí)心圓筒在對(duì)稱軸處溫度隨時(shí)間的變化。顯然，圖7也反映了這種溫度隨著梯度參數(shù)增加而上升的現(xiàn)象。

圖7　功能梯度實(shí)心圓筒在對(duì)稱軸處溫度隨時(shí)間的變化Fig.7　Time variations of the temperature on the axis of a functionally graded full cylinder

5結(jié)論

作為介于有限元法與無(wú)網(wǎng)格法之間的一種數(shù)值方法，自然單元法的節(jié)點(diǎn)影響域是由節(jié)點(diǎn)的Voronoi結(jié)構(gòu)所規(guī)定的自然相鄰關(guān)系給出，不受人為參數(shù)的影響，具有其他無(wú)網(wǎng)格法不可比擬的優(yōu)越性。根據(jù)三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)方程及其邊界條件，利用加權(quán)殘值法，選取自然鄰近插值對(duì)軸對(duì)稱面上的溫度場(chǎng)進(jìn)行離散，首次詳細(xì)推導(dǎo)了三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題的自然單元法計(jì)算公式，并編制了相應(yīng)的FORTRAN計(jì)算程序。本文分析和算例求解結(jié)果表明，采用自然單元法求解三維軸對(duì)稱功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題是可行的，具有精度高和穩(wěn)定性好的優(yōu)點(diǎn)。本文方法還可以容易地推廣到三維軸對(duì)稱功能梯度材料熱彈性問(wèn)題的求解計(jì)算。

參考文獻(xiàn)：

[1] 張靜華，魏軍楊. DQ法求解FGM Levinson梁的靜態(tài)彎曲問(wèn)題[J]. 華東交通大學(xué)學(xué)報(bào), 2014, 31(3): 80-87.

ZHANG J H, WEI J Y. Static bending solution of functionally graded material Levinson beams based on differential quadrature method [J]. Journal of East China Jiaotong University, 2014, 31(3): 80-87. (in Chinese)

[2] WANG H, QIN Q H, KANG Y L. A meshless model for transient heat conduction in functionally graded materials [J]. Computational Mechanics, 2006, 38: 51-60.

[3] 藍(lán)林華，富明慧，程正陽(yáng). 功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題的降維精細(xì)積分法[J]. 固體力學(xué)學(xué)報(bào), 2010, 31(4): 406-410.

LAN L H, FU M H, CHENG Z Y. Decrement-dimensional precise time integration of 2D transient heat conduction equation for functionally graded materials [J]. Chinese Journal of Solid Mechanics, 2010, 31(4): 406-410. (in Chinese)

[4] 陳建橋，丁亮. 功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題的MLPG方法[J]. 華中科技大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 35(4): 119-121.

CHEN J Q, DING L. A MLPG method of transient heat transfer in FGMs [J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2007, 35(4): 119-121. (in Chinese)

[5] 李慶華，陳莘莘，薛志清，等. 基于自然元法和精細(xì)積分法的功能梯度材料瞬態(tài)熱傳導(dǎo)問(wèn)題求解[J]. 計(jì)算機(jī)輔助工程, 2011, 20(3): 25-28.

LI Q H, CHEN S S, XUE Z Q, et al. Transient heat conduction solution to functionally graded materials based on natural element method and precise integration method [J]. Computer Aided Engineering, 2011, 20(3): 25-28. (in Chinese)

[6] LI Q H, CHEN S S, ZENG J H. A meshless model for transient heat conduction analyses of 3D axisymmetric functionally graded solids [J]. Chinese Physics B, 2013, 22: 120204.

[7] SLADEK J, SLADEK V, HELLMICH C h, et al. Heat conduction analysis of 3D axisymmetric and anisotropic FGM bodies by meshless local Petrov-Galerkin method [J]. Computational Mechanics, 2007, 39: 323-333.

[8] SLADEK J, SLADEK V, KRIVACEK J, et al. Local BIEM for transient heat conduction analysis in 3D axisymmetric functionally graded solids [J]. Computational Mechanics, 2003, 32: 169-176.

[9] SUKUMAR N, MORAN T, BELYTSCHKO T. The natural element method in solid mechanics [J]. International Journal for Numerical Methods in Engineering, 1998, 43(5): 839-887.

[10] SUKUMAR N, MORAN B. C1 natural neighbour interpolation for partial differential equations [J]. Numerical Methods for Partial Differential Equations, 1999, 15(4): 417-447.

[11] 江濤，章青. 直接增強(qiáng)自然單元法計(jì)算應(yīng)力強(qiáng)度因子[J]. 計(jì)算力學(xué)學(xué)報(bào), 2010, 27(2): 264-269.

JIANG T, ZHANG Q. Computing stress intensity factors by enriched natural element method [J]. Chinese Journal of Computational Mechanics, 2010, 27(2): 264-269. (in Chinese)

[12] 張勇，易紅亮，談和平. 求解輻射導(dǎo)熱耦合換熱的自然單元法[J]. 工程熱物理學(xué)報(bào), 2013, 34(5): 918-922.

ZHANG Y, YI H L, TAN H P. Natural element method for coupled radiative and conductive heat transfer [J]. Journal of Engineering Thermophysics, 2013, 34(5): 918-922. (in Chinese)

[13] 董軼，馬永其，馮偉. 彈性力學(xué)的雜交自然單元法[J]. 力學(xué)學(xué)報(bào), 2012, 44(3): 568-575.

DONG Y, MA Y Q, FENG W. The hybrid natural element method for elasticity [J]. Chinese Journal of Theoretical and Applied Mechanics, 2012, 44(3): 568-575. (in Chinese)

[14] 曾祥勇，朱愛(ài)軍，鄧安福. Winkler地基上厚板分析的自然單元法[J]. 固體力學(xué)學(xué)報(bào), 2008, 29(2): 163-169.

ZENG X Y, ZHU A J, DENG A F. Natural element method for analysis of thick plates lying over Winkler foundation [J]. Chinese Journal of Solid Mechanics, 2008, 29(2): 163-169. (in Chinese)

[15] 周書(shū)濤，劉應(yīng)華，陳莘莘. 基于自然單元法的極限上限分析[J]. 固體力學(xué)學(xué)報(bào), 2012, 33(1): 39-47.

ZHOU S T, LIU Y H, CHEN S S. Upper-bound limit analysis based on natural element method [J]. Chinese Journal of Solid Mechanics, 2012, 33(1): 39-47. (in Chinese)

[16] CHEN L, LIEW K M. A local Petrov-Galerkin approach with moving Kriging interpolation for solving transient heat conduction problems [J]. Computational Mechanics, 2011, 47: 455-467.

[17] CARSLAW H S, JAEGER J C. Conduction of heat in solids [M]. Oxford: Clarendon Press, 1959.

(編輯胡玲)

Natural element method for transient heat conduction analyse of 3D axisymmetric functionally graded solids

Li Qinghua, Chen Shenshen, Xu Qing

(School of Civil Engineering and Architecture, East China Jiaotong University, Nanchang 330013, P.R. China)

Abstract:In order to solve the transient heat conduction problems in three-dimensional (3D) axisymmetric continuously nonhomogeneous functionally graded materials (FGMs) more effectively, a novel numerical method based on the meshless natural element method is proposed. Axial symmetry of geometry and boundary conditions helps to transform the 3D axisymmetric problem into a two-dimensional (2D) prolem. In order to simplify the imposition of the essential boundary conditions, the natural neighbour interpolation is adopted to discretize the temperature field within the cross section. The variations of functionally graded material properties are simulated by employing proper material parameters at Gauss points. The spatially discretized heat conduction equation is solved numerically with the traditional two-point difference technique in the time domain. The present method not only broadens the application scope of the natural element method, but also will be generally available to transient heat conduction analyses of 3D axisymmetric solids.

Keywords:natural element method; axisymmetric; functionally graded materials; transient heat conduction

doi:10.11835/j.issn.1674-4764.2016.02.009 10.11835/j.issn.1674-4764.2016.02.010

收稿日期：2015-11-15

基金項(xiàng)目：國(guó)家自然科學(xué)基金(21466012)；江西省教育廳項(xiàng)目(KJLD14041).

作者簡(jiǎn)介：李慶華(1976-)，女，主要從事數(shù)值傳熱學(xué)研究，(E-mail) jessylqh@126.com。

中圖分類號(hào)：TP301.6

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1674-4764(2016)02-0069-06

Received:2015-11-15

Foundation item：National Natural Science Foundation of China (No.21466012); Foundation of Jiangxi Educational Committee (No.KJLD14041)

Author brief：Li Qinghua(1976-), main research interest: numerical heat transfer, (E-mail) jessylqh@126.com.