功能梯度材料熱傳導(dǎo)的統(tǒng)計(jì)多尺度邊界元分析

黃志強(qiáng)，韓要闖，聶玉峰

(1.西北工業(yè)大學(xué) 理學(xué)院，西安 710072) (2.太原科技大學(xué) 應(yīng)用科學(xué)學(xué)院，太原 030024)

0 引 言

功能梯度材料[1]是一種在材料的制備過程中通過連續(xù)地控制各組分含量的分布使材料宏觀特性在空間位置上呈現(xiàn)梯度變化的非均質(zhì)多相復(fù)合材料。由于功能梯度材料的結(jié)構(gòu)和性能呈連續(xù)變化，不存在明顯的界面，其熱力學(xué)性能和物理性能也呈梯度變化趨勢(shì)，滿足了結(jié)構(gòu)元件不同部位對(duì)材料使用性能的不同要求，進(jìn)而達(dá)到優(yōu)化結(jié)構(gòu)整體使用性能的目的。面對(duì)不斷增長(zhǎng)的工程應(yīng)用需求，如何有效表征功能梯度材料的細(xì)觀結(jié)構(gòu)并建立細(xì)觀結(jié)構(gòu)與宏觀熱傳導(dǎo)性能之間的定量關(guān)系，已成為材料科學(xué)與工程領(lǐng)域的重要課題。

研究功能梯度材料宏觀等效性能的細(xì)觀力學(xué)解析方法眾多，包括Eshelby等效夾雜理論[2]、自洽模型[3]、廣義自洽模型[4]、Mori-Tanaka模型[5]、微分法[6]、變分法[7]和n點(diǎn)界[8]等，上述模型均對(duì)材料微結(jié)構(gòu)進(jìn)行了大量簡(jiǎn)化以減少計(jì)算規(guī)模，故不能充分反映材料的真實(shí)微觀結(jié)構(gòu)特征。M.Jabbari等[9]采用直接方法分析了一維和二維穩(wěn)態(tài)功能梯度空心圓筒熱傳導(dǎo)和熱彈性問題；A.Alibeigloo[10]采用傅里葉級(jí)數(shù)法推導(dǎo)了物性參數(shù)沿厚度方向呈指數(shù)函數(shù)變化的溫度場(chǎng)解析解；陶光勇等[11]采用解析法研究了功能梯度材料板在穩(wěn)態(tài)梯度溫度場(chǎng)下的熱應(yīng)力分布狀況；Shao Z S等[12]利用級(jí)數(shù)求解法給出了功能梯度圓筒穩(wěn)態(tài)溫度場(chǎng)的解析解。采用解析法研究功能梯度材料熱傳導(dǎo)問題已取得了一些研究成果，但是由于功能梯度材料的非均勻性，利用解析法研究其熱傳導(dǎo)問題具有很大的局限性。Cao Leilei等[13]采用基本解有限梯度元方法求解了功能梯度材料的穩(wěn)態(tài)熱傳導(dǎo)問題；王魯?shù)萚14]采用有限元法研究了功能梯度熱障涂層在循環(huán)狀態(tài)下的瞬態(tài)溫度分布。由于功能梯度材料性能的不均勻性，直接使用有限元等傳統(tǒng)數(shù)值方法時(shí)需要非常精細(xì)的網(wǎng)格去捕捉功能梯度材料的局部特征，這將導(dǎo)致計(jì)算規(guī)模的大幅增加、求解困難。因此，建立一種針對(duì)功能梯度材料的高效可行的數(shù)值算法是十分必要的。

近年來(lái)，基于均勻化方法[15]的多尺度分析方法被成功應(yīng)用于預(yù)測(cè)不同復(fù)合材料結(jié)構(gòu)的物理和力學(xué)性能[16-19]。Li Youyun等[16]通過引入隨機(jī)樣本單胞模型，提出了基于有限元方法的統(tǒng)計(jì)多尺度分析方法，用于預(yù)測(cè)顆?；蚩锥措S機(jī)分布復(fù)合材料的熱傳導(dǎo)性能?？紤]到功能梯度材料微結(jié)構(gòu)的隨機(jī)性，多尺度計(jì)算結(jié)果呈現(xiàn)很大的發(fā)散性，必須進(jìn)行大量樣本的計(jì)算，最后給出統(tǒng)計(jì)意義上的結(jié)果，才能真實(shí)地反映隨機(jī)復(fù)合材料的宏觀整體導(dǎo)熱性能。此外，功能梯度材料中顆粒數(shù)量巨大，尺寸較小，且體積分?jǐn)?shù)隨著位置變化而變化，結(jié)構(gòu)十分復(fù)雜，有限元體網(wǎng)格剖分難度大；同時(shí)，為了更好地逼近材料微觀結(jié)構(gòu)，需要數(shù)量巨大的有限元網(wǎng)格，導(dǎo)致計(jì)算規(guī)模龐大。與有限元方法不同，邊界元方法只需要在功能梯度材料各相材料的界面上進(jìn)行網(wǎng)格剖分[20]，網(wǎng)格數(shù)目少，劃分難度小，邊界元方法在求解功能梯度材料導(dǎo)熱問題上具有顯著優(yōu)勢(shì)。

本文針對(duì)功能梯度材料的穩(wěn)態(tài)熱傳導(dǎo)問題，描述功能梯度材料微結(jié)構(gòu)的表征方法，建立宏觀等效材料參數(shù)的多尺度邊界元計(jì)算模型并給出統(tǒng)計(jì)意義下的多尺度算法，并研究顆粒隨機(jī)分布功能梯度材料的熱傳導(dǎo)性能。

1 統(tǒng)計(jì)多尺度分析模型

圖1 功能梯度材料的細(xì)觀結(jié)構(gòu)

(1)

根據(jù)上述細(xì)觀結(jié)構(gòu)表征，考慮隨機(jī)分布復(fù)合材料結(jié)構(gòu)的穩(wěn)態(tài)熱傳導(dǎo)方程為

(2)

Tε(x,ωs)=T0(x,y,ωx′)+εT1(x,y,ωx′)+

ε2T2(x,y,ωx′)+…

(3)

由于y=x/ε，則存在鏈?zhǔn)椒▌t：

(4)

將式(3)～式(4)帶入式(2)，并整理成ε冪級(jí)數(shù)的形式，可得：

(5)

通過比較式(5)兩端ε不同冪次的系數(shù)，根據(jù)偏微分方程理論可分別定義T0和T1，則溫度場(chǎng)的多尺度漸進(jìn)展開式可定義為

(6)

(7)

(8)

(9)

利用上述期望的均勻化系數(shù)，可以定義功能梯度材料結(jié)構(gòu)的均勻化熱傳導(dǎo)方程：

(10)

2 統(tǒng)計(jì)多尺度邊界元算法

由于功能梯度材料中顆粒數(shù)量大、滿足隨機(jī)性分布且體積分?jǐn)?shù)隨著位置而變化，其微觀單胞結(jié)構(gòu)十分復(fù)雜，為了有效求解定義在單胞Yx′上的單胞函數(shù)，采用有限元方法時(shí)需要大量的體網(wǎng)格來(lái)逼近單胞結(jié)構(gòu)，計(jì)算規(guī)模巨大。此外，由于顆粒分布的隨機(jī)性，必須進(jìn)行大量樣本的計(jì)算，最后計(jì)算出統(tǒng)計(jì)意義上的結(jié)果，才能真實(shí)地反映隨機(jī)復(fù)合材料的宏觀整體性能。因此，采用有限元算法求解期望均勻化系數(shù)時(shí)計(jì)算規(guī)模巨大。而邊界元方法在求解該問題時(shí)只需要在各組分材料的邊界上進(jìn)行網(wǎng)格離散，使得整體網(wǎng)格數(shù)量大幅減少，有效地降低了計(jì)算規(guī)模。假設(shè)功能梯度材料中基體和顆粒的熱傳導(dǎo)性能滿足各向同性，則基于邊界元方法的統(tǒng)計(jì)多尺度算法如下：

(1) 根據(jù)給定的概率分布模型P(x)(x∈Ω)，在結(jié)構(gòu)Ω上選擇M個(gè)點(diǎn)xi(i=1,2,…,M)。

(11)

(12)

式中：c(y)為與邊界有關(guān)的常量(對(duì)于光滑邊界取值1/2)；v為單位法向量。

進(jìn)一步對(duì)式(12)等號(hào)右端的第三項(xiàng)采用徑向積分，可得：

(13)

式中：r(p,y)為p點(diǎn)和y之間的距離。

(14)

式中：φi(y)為徑向基函數(shù)；Φ為徑向基函數(shù)插值矩陣。

然后采用徑向積分方法將區(qū)域積分轉(zhuǎn)換為邊界積分，即

(15)

(16)

3 算例分析

為了驗(yàn)證統(tǒng)計(jì)多尺度邊界元方法預(yù)測(cè)功能梯度材料熱傳導(dǎo)性能的有效性，取圓球顆粒Al作為增強(qiáng)相顆粒的Al/Al2O3功能梯度材料，Al和Al2O3的熱導(dǎo)率分別為211.07和38.27 W/mK，圓球顆粒在基體內(nèi)均勻隨機(jī)分布且體積分?jǐn)?shù)隨著位置從0～25%線性變化。由于顆粒分布的隨機(jī)性，數(shù)值計(jì)算結(jié)果會(huì)受此影響，即使顆粒的分布模型完全相同，計(jì)算結(jié)果也會(huì)因樣本的不同而有所差別。因此，一次樣本的計(jì)算結(jié)果無(wú)法完全反映隨機(jī)復(fù)合材料的宏觀特性，必須進(jìn)行大量樣本的計(jì)算，再最終計(jì)算出統(tǒng)計(jì)意義上的結(jié)果，才能真實(shí)反映隨機(jī)復(fù)合材料的宏觀整體性能。算例中的計(jì)算結(jié)果均由隨機(jī)抽樣50次統(tǒng)計(jì)所得。

選取不同樣本數(shù)量時(shí)的等效導(dǎo)熱系數(shù)計(jì)算結(jié)果如圖2所示，可以看出：隨著樣本數(shù)量的增加，計(jì)算結(jié)果分散性越來(lái)越小。

圖2 樣本數(shù)量對(duì)等效導(dǎo)熱系數(shù)的影響

分別使用多尺度有限元和多尺度邊界元方法計(jì)算等效導(dǎo)熱系數(shù)的數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果的比較如圖3所示，可以看出：在功能梯度材料內(nèi)部有效熱傳導(dǎo)系數(shù)隨著位置的變化呈上升趨勢(shì)，兩種多尺度方法均給出了穩(wěn)定的預(yù)測(cè)曲線，與實(shí)驗(yàn)數(shù)據(jù)[22]符合較好。

圖3 等效導(dǎo)熱系數(shù)數(shù)值結(jié)果與實(shí)驗(yàn)數(shù)據(jù)比較

多尺度有限元網(wǎng)格與多尺度邊界元網(wǎng)格數(shù)目的比較如圖4所示，可以看出：邊界元網(wǎng)格數(shù)目遠(yuǎn)少于有限元網(wǎng)格數(shù)目，采用多尺度邊界元算法可以大幅節(jié)約計(jì)算成本，減少計(jì)算量。特別地，針對(duì)隨機(jī)復(fù)合材料模擬需要大量取樣的特性，多尺度邊界元算法節(jié)約的計(jì)算量是相當(dāng)可觀的。

圖4 有限元網(wǎng)格與邊界元網(wǎng)格數(shù)比較

從圖3～圖4可以看出：統(tǒng)計(jì)多尺度邊界元方法可以用較少的計(jì)算時(shí)間得到滿意的結(jié)果，是一種預(yù)測(cè)功能梯度材料熱傳導(dǎo)性能的高效率、高精度的數(shù)值方法。

4 結(jié) 論

(1) 使用統(tǒng)計(jì)多尺度邊界元方法預(yù)測(cè)了顆粒隨機(jī)分布功能梯度材料的等效導(dǎo)熱系數(shù)，數(shù)值結(jié)果與實(shí)驗(yàn)結(jié)果的比較表明了多尺度邊界元模型的正確性和算法的有效性。

(2) 針對(duì)具有大量顆粒隨機(jī)分布且結(jié)構(gòu)和性能在空間位置上呈現(xiàn)梯度變化的功能梯度材料的熱傳導(dǎo)問題，統(tǒng)計(jì)多尺度邊界元算法可以用較少的計(jì)算時(shí)間得到滿意的結(jié)果，是一種高效率、高精度的數(shù)值算法。

[1] 孫蘭. 功能材料及應(yīng)用[M]. 成都: 四川大學(xué)出版社, 2015.

Sun Lan. Functional materials and applications[M]. Chengdu: Sichuan University Press, 2015.(in Chinese)

[2] Eshelby J D. The determination of the elastic field of an ellipsoidal inclusion and related problems[J]. Proceedings of the Royal Society, 1957, 2: 267-396.

[3] Budiansky B. On the elastic moduli of some heterogeneous materials[J]. Journal of the Mechanics and Physics of Solids, 1965, 13(4): 223-227.

[4] Kerner E H. The electrical conductivity of composite media[J]. Proceedings of the Physical Society, 1956(69): 802-808.

[5] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions[J]. Acta Metallurgica, 1973, 21(5): 571-574.

[6] Mclaughlin R. A study of the differential scheme for composite materials[J]. International Journal of Engineering Science, 1977(15): 237-244.

[7] Hashin Z, Shtrikman S. A variational approach to the theory of the elastic behaviour of multiphase materials[J]. Journal of the Mechanics and Physics of Solids,1963, 11: 86-93.

[8] Torquato S. Random heterogeneous materials: microstructure and macroscopic properties[M]. New York: Springer, 2002.

[9] Jabbari M, Sorabpour S, Eslami M R. General solution for mechanical and thermal stresses in a functionally graded hollow cylinder due to nonaxisymmetric steady-state loads[J]. Journal of Applied Mechanics, 2003, 70: 111-117.

[10] Alibeigloo A. Exact solution for thermo-elastic response of functionally graded rectangular plates[J]. Composite Structures, 2010, 92(1): 113-121.

[11] 陶光勇, 鄭子樵, 劉孫和. W/Cu梯度功能材料板穩(wěn)態(tài)熱應(yīng)力分析[J]. 中國(guó)有色金屬學(xué)報(bào), 2006, 16(4): 694-700.

Tao Guangyong, Zheng Ziqiao, Liu Sunhe. Analysis of steady thermal stress in plate of W/Cu functionally graded material[J]. The Chinese Journal of Nongerrous Metals, 2006, 16(4): 694-700.(in Chinese)

[12] Z S Shao, L F Fan, T J Wang. Analytical solutions of stresses in functionally graded circular hollow cylinder with finite length[J]. Key Engineering Materials, 2004(261/263): 651-656.

[13] Cao Leilei, Wang Hui, Qin Qinghua. Fundamental solution based graded element model for steady-state heat transfer in FGM[J]. Acta Mechanica Solida Sinica, 2012, 25(4): 377-378.

[14] 王魯, 呂廣庶, 王富恥, 等. 功能梯度熱障涂層熱負(fù)荷下的有限元分析[J]. 兵工學(xué)報(bào), 1999, 20(1): 51-54.

Wang Lu, Lu Guangshu, Wang Fuchi, et al. Finite element analysis for functionally gradient thermal barrier coatings under thermal loads[J]. Acta Armamentarii, 1999, 20(1): 51-54.(in Chinese)

[15] Bensoussan A, Lions J L, Papanicolaou G. Asymptotic analysis for periodic structures[M]. Rhode Island: American Mathematical Society, 2011.

[16] Li Youyun, Long Shuyao. A finite element model based on statistical two-scale analysis for equivalent heat transfer parameters of composite material with random grains[J]. Applied Mathematical Modelling, 2009, 33(7): 3157-3165.

[17] Han Fei, Cui Junzhi, Yu Yan. The statistical second-order two-scale method for thermomechanical properties of statistically inhomogeneous materials[J]. Computational Materials Science, 2009, 46: 654-659.

[18] Dong Hao, Cui Junzhi, Nie Yufeng. Second-order two-scale computational method for nonlinear dynamic thermo-mechanical problems of composites with cylindrical periodicity[J]. Communications in Computational Physics, 2017, 21(4): 1173-1206.

[19] Yang Zihao, Ma Liang, Ma Qiang. Multiscale nonlinear th-ermo-mechanical coupling analysis of composite structures with quasi-periodic properties[J]. CMC-Computers Materials & Continua, 2017, 53: 219-248.

[20] 高效偉, 王靜, 彭海峰, 等. 高等邊界單元法: 理論與程序[M]. 北京: 科學(xué)出版社, 2015.

Gao Xiaowei, Wang Jing, Peng Haifeng, et al. Advanced boundary element method: theory and program[M]. Beijing: Science Press, 2015.(in Chinese)

[21] Yu Yan, Cui Junzhi, Han Fei. The statistical second-order two-scale analysis method for heat conduction performances of the composite structure with inconsistent random distribution[J]. Computational Materials Science, 2009, 46(1): 151-161.

[22] Neubrand A, Becker H, Tschudi T. Spatially resolved thermal diffusivity measurements on functionally graded materials[J]. Journal of Materials Science, 2003, 38(20): 4193-4201.