李義強，楊自豪，楊志強

(1.西北工業(yè)大學(xué) 理學(xué)院，西安 710072)(2.哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001)

隨機復(fù)合材料統(tǒng)計多尺度邊界元算法研究

李義強1，楊自豪1，楊志強2

(1.西北工業(yè)大學(xué) 理學(xué)院，西安 710072)(2.哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001)

熱傳輸行為是空天飛行器熱防護系統(tǒng)服役過程中的重要問題，研究隨機復(fù)合材料的熱傳輸機理，可為熱防護系統(tǒng)材料與結(jié)構(gòu)的一體化設(shè)計提供理論依據(jù)。針對隨機復(fù)合材料的熱傳導(dǎo)問題，發(fā)展一種統(tǒng)計多尺度邊界元算法，首先建立等效材料參數(shù)的統(tǒng)計多尺度分析模型，然后給出統(tǒng)計意義下等效參數(shù)的多尺度邊界元預(yù)測算法，并通過與理論結(jié)果的對比來驗證算法的有效性，最后研究微結(jié)構(gòu)分布狀態(tài)對陶瓷多孔材料等效導(dǎo)熱系數(shù)的影響。結(jié)果表明：采用統(tǒng)計多尺度模型及多尺度邊界元算法預(yù)測隨機復(fù)合材料的熱傳導(dǎo)性能是有效的。

隨機復(fù)合材料；熱傳導(dǎo)性能；統(tǒng)計多尺度邊界元算法；等效導(dǎo)熱系數(shù)

0 引 言

實際材料的微細(xì)觀構(gòu)造通常是模糊不清的，隨機不確定性是很多天然/工程材料的固有屬性。復(fù)合材料的增強相離散地分布于基體中，增強相的幾何尺寸、空間方位等多具有隨機分布的特征。隨機復(fù)合材料是工程中常用的一類復(fù)合材料。面對不斷增長的工程應(yīng)用需求，如何有效表征隨機復(fù)合材料的細(xì)觀結(jié)構(gòu)，如何定量描述復(fù)合材料的局部細(xì)觀構(gòu)造特征對材料宏觀性能的影響，已成為材料科學(xué)與工程領(lǐng)域的重要課題。

針對隨機復(fù)合材料的熱傳導(dǎo)性能預(yù)測問題，國內(nèi)外已開展了諸多研究，例如，Eshelby等效夾雜方法[1]、Hashin-Shtrikman(HS)上下界法[2]、自洽場方法[3]和Mori-Tanaka理論[4]等，但上述平均意義下的方法都對真實的材料結(jié)構(gòu)進行了較大程度的簡化以減少計算量，并不能充分反映材料的真實微結(jié)構(gòu)特征。隨著計算機技術(shù)的發(fā)展，針對顆粒填充或多孔復(fù)合材料的熱學(xué)性能模擬預(yù)測問題，以有限元分析為基礎(chǔ)的數(shù)值方法采用了規(guī)則排布的單胞模型[5]，但仍然存在較大程度的簡化。然而，顆粒/多孔隨機復(fù)合材料需要更加真實、復(fù)雜的微結(jié)構(gòu)單胞模型以及快速有效的熱學(xué)性能預(yù)測方法，以便能夠更精確地描述材料的實際熱學(xué)行為。

近年來，針對不同類型的復(fù)合材料及其結(jié)構(gòu)，崔俊芝等[6-8]基于均勻化方法[9]，發(fā)展了一系列高階雙尺度計算方法，成功地預(yù)測了復(fù)合材料結(jié)構(gòu)的物理和力學(xué)性能?；诖?，Guan Xiaofei等[10]、Cui Junzhi等[11-12]通過引入均勻隨機樣本單胞模型，提出了統(tǒng)計二階雙尺度分析方法及其有限元方法，用于預(yù)測顆?；蚩锥措S機分布復(fù)合材料的物理和力學(xué)性能。由于隨機樣本單胞的分散性，為了獲得更為準(zhǔn)確的等效熱學(xué)性能，需要大量取樣計算等效導(dǎo)熱系數(shù)并求均值。此外，在使用傳統(tǒng)的有限元方法求解多尺度單胞模型時，由于顆粒或孔洞數(shù)量巨大、尺寸較小、結(jié)構(gòu)復(fù)雜，極大地增加了有限元網(wǎng)格剖分的難度，并且為了更好地逼近單胞結(jié)構(gòu)，需要數(shù)量巨大的有限元網(wǎng)格，導(dǎo)致極大的計算量。

與有限元方法不同，邊界元算法只需要在邊界和界面上進行網(wǎng)格剖分[13-14]，單胞結(jié)構(gòu)網(wǎng)格劃分難度小、網(wǎng)格數(shù)目少，因此，邊界元算法在求解隨機復(fù)合材料單胞模型時具有明顯優(yōu)勢。但是，邊界元算法在直接求解復(fù)合材料物理力學(xué)問題時，存在諸多困難，特別是在求解三維問題時尚未有成熟的方法。本文針對多孔隨機分布復(fù)合材料的熱傳導(dǎo)問題，發(fā)展一種統(tǒng)計多尺度邊界元算法，并應(yīng)用該方法計算復(fù)合材料結(jié)構(gòu)的熱學(xué)參數(shù)，分析其熱學(xué)性能，以期為研究隨機復(fù)合材料的熱傳輸機理提供有益參考。

1 隨機復(fù)合材料微結(jié)構(gòu)表征

設(shè)隨機復(fù)合材料結(jié)構(gòu)Ω由顆粒和基體組成，將顆?；蚩锥吹男螤罱y(tǒng)一設(shè)為橢球或者嵌入橢球的多面體，并假設(shè)材料內(nèi)部顆粒或孔洞的隨機分布模型處處相同。在三維空間中，一個橢球可以由9個參數(shù)(橢球的中心點，長、中、短軸的長度和方向角)唯一確定?；诠こ虦y量和統(tǒng)計方法，整個復(fù)合材料結(jié)構(gòu)Ω可以邏輯地分為一系列具有相同尺寸ε的單胞體[15]，如圖1所示。

因此，只要確定一個單胞內(nèi)橢球的概率分布模型，也就確定了整個復(fù)合材料結(jié)構(gòu)的細(xì)觀特征。按照上述表征方法，隨機分布復(fù)合材料結(jié)構(gòu)Ω可以看作是由具有相同顆?；蚩锥捶植寄Ｐ偷膯伟M合而成，即

(1)

式中：ε(Qs+zs)為區(qū)域Ω內(nèi)的第s個單胞；P為復(fù)合材料的概率分布模型。

(2)

式中：ei為單胞εQs內(nèi)的第i個橢球顆粒；A1和A2分別為顆粒和基體的材料參數(shù)。

2 統(tǒng)計多尺度分析模型

對于隨機分布復(fù)合材料，考慮如下具有混合邊界條件的穩(wěn)態(tài)熱傳導(dǎo)問題：

(3)

θε(x,ωs)= θ0(x,ξ,ωs)+εθ1(x,ξ,ωs)+

ε2θ2(x,ξ,ωs)+…

(4)

由于ξ=x/ε，存在如下鏈?zhǔn)椒▌t：

(5)

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

O(ε)=f

(6)

比較式(6)兩端不同ε冪次的系數(shù)，通過偏微分方程理論可分別定義θ0，θ1和θ2，則溫度場的多尺度漸進展開式可定義為

(7)

式中：θ0(x)為定義在Ω上的均勻化解；Hα1(ξ,ωs)和Hα1α2(ξ,ωs)為定義在單位化單胞Qs上的依賴于樣本ωs的局部單胞解。

Hα1(ξ,ωs)是控制方程式(8)的解。

(8)

通過式(9)可得與樣本ωs有關(guān)的均勻化熱傳導(dǎo)系數(shù)：

(9)

根據(jù)柯爾莫哥洛夫強大數(shù)定理，期望均勻化系數(shù)為

(10)

(11)

Hα1α2(ξ,ωs)是式(12)的解。

(12)

基于溫度場的多尺度計算模型(式(7))，可以推導(dǎo)出熱流密度的多尺度計算公式：

(13)

3 統(tǒng)計多尺度邊界元算法

隨機單胞結(jié)構(gòu)十分復(fù)雜，為了有效求解單胞問題，采用有限元方法需要大量的體網(wǎng)格來逼近單胞結(jié)構(gòu)，使得計算規(guī)模巨大。在統(tǒng)計多尺度分析模型中，計算期望均勻化系數(shù)時需要求解大量單胞問題，導(dǎo)致采用有限元方法時會產(chǎn)生巨大的計算量。與有限元方法不同，邊界元算法只需要在單胞邊界上進行離散，不需要體網(wǎng)格，有效地減小了計算規(guī)模，提高了計算效率。假設(shè)復(fù)合材料的基體和顆粒/孔洞的導(dǎo)熱性質(zhì)滿足各向同性，本節(jié)將給出統(tǒng)計多尺度邊界元算法及其計算流程。

3.1 單胞問題邊界元算法

推導(dǎo)單胞函數(shù)Hα1(ξ,ωs)滿足的式(8)的邊界積分方程形式，該方程類似于Laplace方程，取其基本解為[13]

(14)

取kij(ξ,ωs)=k(ξ,ωs)，并以基本解G(ξ,τ)為權(quán)函數(shù)，對式(8)進行加權(quán)積分，得：

(15)

利用高斯散度定理和基本解的性質(zhì)[14]，對式(15)的第一個積分項進行分部積分：

=-c(τ)k(τ,ωs)Hα1(τ,ωs)-

∫?QsG(ξ,τ)q(ξ)dS-

(τ∈?Qs)

(16)

因此，單胞函數(shù)Hα1(ξ,ωs)滿足的積分方程為

-c(τ)k(τ,ωs)Hα1(τ,ωs)

∫QsG(ξ,τ)bα1(ξ)dξ-

(τ∈?Qs)

(17)

其中，

由于式(17)仍含有區(qū)域積分，需要做進一步處理。對于式(17)等號右端第三項，采用徑向積分方法[14]可得：

(18)

式中：rα(p,ξ)為p點和ξ之間的距離，對于二維問題，α=2，對于三維問題，α=3。

再處理式(17)等號右端第四項，記：

(19)

(20)

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

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

(21)

綜上所述，便可得到單胞函數(shù)所滿足的邊界積分方程。通過對單胞Qs的邊界進行網(wǎng)格剖分，并進行邊界單元插值，計算各項積分，得到線性代數(shù)方程組，求解可得單胞問題的數(shù)值解。

3.2 算法流程

基于統(tǒng)計多尺度邊界元算法的隨機復(fù)合材料熱傳導(dǎo)性能的算法流程為：

(1) 根據(jù)給定的概率分布模型P，生成一個隨機樣本單胞Q(ωs)(ωs∈P)[15]，并對其進行網(wǎng)格剖分[16]；

(2) 利用邊界元算法求解單胞問題(式(8)和式(12))；

(4) 根據(jù)求解出的期望均勻化系數(shù)，求解均勻化問題(式(11))，得到宏觀溫度解θ0(x)；

(5) 對應(yīng)樣本單胞Q(ωs)，基于已求出的單胞函數(shù)解和宏觀均勻化溫度解，通過多尺度計算公式(式(7)和式(13))，可以分別確定單胞εQs內(nèi)的溫度場、熱流密度場，還可以進一步計算出熱流密度極值。

4 算例分析

為了研究隨機復(fù)合材料的微觀構(gòu)造對材料宏觀導(dǎo)熱性能的影響，考慮三種微觀孔洞分布狀態(tài)：圓球孔洞中心均勻分布，圓球孔洞中心以單胞的形心點為中心正態(tài)分布，長橢球孔洞(長軸為中軸和短軸的2倍)的傾角沿x3方向正態(tài)分布，且中心在單胞內(nèi)均勻分布，如圖2所示。

由于孔洞分布的隨機性，即使孔洞的分布模型完全相同，其數(shù)值計算結(jié)果也會因樣本的不同而有所差別。為了得到較為準(zhǔn)確的預(yù)測結(jié)果，必須進行大量樣本計算，再計算出統(tǒng)計意義下的結(jié)果。本文的計算結(jié)果均為隨機抽樣50次統(tǒng)計所得。

當(dāng)基體材料為陶瓷，導(dǎo)熱系數(shù)為4.41 W/mK時，將采用統(tǒng)計多尺度邊界元算法預(yù)測的期望均勻化導(dǎo)熱系數(shù)與采用Hashin-Shtrikman(HS)上下界法的結(jié)果進行比較，結(jié)果如表1～表3所示。

表1 圓球孔洞均勻分布不同體積分?jǐn)?shù)時均勻化導(dǎo)熱系數(shù)的多尺度解與HS界的比較

表2 圓球孔洞正態(tài)分布不同體積分?jǐn)?shù)時均勻化導(dǎo)熱系數(shù)的多尺度解與HS界的比較

表3 橢球孔洞傾角正態(tài)分布不同體積分?jǐn)?shù)時均勻化導(dǎo)熱系數(shù)的多尺度解與HS界的比較

從表1～表3可以看出：對于不同體積分?jǐn)?shù)和分布狀態(tài)的多孔隨機分布復(fù)合材料，其多尺度數(shù)值解均在HS界內(nèi)，且對于相同分布狀態(tài)的復(fù)合材料結(jié)構(gòu)，其宏觀均勻化系數(shù)隨孔洞體積分?jǐn)?shù)的增大而發(fā)生較大變化，表明孔洞的體積分?jǐn)?shù)、結(jié)構(gòu)形態(tài)以及分布狀態(tài)對復(fù)合材料的宏觀等效導(dǎo)熱性能具有重要影響。

多尺度有限元網(wǎng)格數(shù)目與多尺度邊界元網(wǎng)格數(shù)目的比較如表4所示，可以看出：邊界元網(wǎng)格的數(shù)目遠(yuǎn)少于有限元網(wǎng)格的數(shù)目，即采用多尺度邊界元算法相比有限元方法可減少相當(dāng)大的計算量，尤其是在大量取樣計算過程中，采用邊界元算法進行計算所節(jié)約的計算量是相當(dāng)可觀的。

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

當(dāng)選取陶瓷作為基體材料時，其導(dǎo)熱系數(shù)為21.16W/mK時，取樣一次，分別采用有限元方法和邊界元算法計算出的均勻化系數(shù)如表5所示(兩種算法的計算精度基本一致)。

表5 圓球孔洞均勻分布時均勻化導(dǎo)熱系數(shù)計算結(jié)果

從表4～表5可以看出：多尺度邊界元算法僅用較少的時間就能得到滿意的結(jié)果，是一種高效率、高精度的數(shù)值算法。

選取不同樣本數(shù)量時的期望均勻化系數(shù)計算結(jié)果如圖3所示，可以看出：計算結(jié)果的分散性隨著樣本數(shù)量的增加而減小。

5 結(jié) 論

(1) 本文采用多尺度邊界元算法預(yù)測了多孔隨機復(fù)合材料的等效導(dǎo)熱系數(shù)，并將數(shù)值計算結(jié)果與理論值進行對比，表明多尺度邊界元算法能夠有效預(yù)測隨機復(fù)合材料的宏觀導(dǎo)熱性能。

(2) 多孔隨機復(fù)合材料的熱傳導(dǎo)性能不僅依賴于宏觀條件，還與復(fù)合材料的細(xì)觀結(jié)構(gòu)有關(guān)，包括孔洞的體積分?jǐn)?shù)、結(jié)構(gòu)形態(tài)以及分布狀態(tài)等，統(tǒng)計多尺度邊界元算法能夠較為準(zhǔn)確地捕捉多孔隨機復(fù)合材料的細(xì)觀結(jié)構(gòu)信息。

(3) 統(tǒng)計多尺度邊界元算法在求解具有大量孔洞隨機分布復(fù)合材料的熱傳導(dǎo)問題時，僅以較少的計算時間就能得到滿意的結(jié)果，是一種高效率、高精度的數(shù)值算法。

[1]EshelbyJD.Thedeterminationoftheelasticfieldofanellipsoidalinclusion,andrelatedproblems[J].ProceedingsoftheRoyalSocietyofLondon, 1957, 241(2): 376-396.

[2]HashinZ,ShtrikmanS.Avariationalapproachtothetheoryoftheelasticbehaviourofmultiphasematerials[J].JournaloftheMechanicsandPhysicsofSolids, 1963, 11(2): 127-140.

[3]HillRW.Theelasticbehaviorofacrystallineaggregate[J].ProceedingsofthePhysicalSociety, 1952, 65(5): 349-354.

[4]MoriT,TanakaK.Averagestressinmatrixandaverageelasticenergyofmaterialswithmisfittinginclusions[J].ActaMetallurgica, 1973, 21(5): 571-574.

[5]KangGuozheng,ShaoXuejiao,GuoSujuan.EffectofinterfacialbondingonuniaxialratchettingofSiCP/6061Alcomposites:finiteelementanalysiswith2-Dand3-Dunitcells[J].MaterialsScienceandEngineering:A, 2008, 487(1/2): 431-444.

[6] 崔俊芝, 曹禮群. 基于雙尺度漸進分析的有限元算法[J]. 計算數(shù)學(xué), 1998, 20(1): 89-103.CuiJunzhi,CaoLiqun.Finiteelementmethodbasedontwo-scaleasymptoticanalysis[J].MathematicaNumericaSinica, 1998, 20(1): 89-103.(inChinese)

[7]YangZihao,CuiJunzhi,WuYatao,etal.Second-ordertwo-scaleanalysismethodfordynamicthermo-mechanicalproblemsinperiodicstructure[J].InternationalJournalofNumericalAnalysis&Modeling, 2015, 12(1): 144-161.

[8]CuiJunzhi,YuXingang.Atwo-scalemethodforidentifyingmechanicalparametersofcompositematerialswithperiodicconfiguration[J].ActaMechanicaSinica, 2006, 22(6): 581-594.

[9]BensoussanA,LionsJL,PapanicolaouG.Asymptoticanalysisforperiodicstructures[M].Amsterdam:North-HollandPublishingCompany, 1978.

[10]GuanXiaofei,LiuXian,JiaXin,etal.Astochasticmultiscalemodelforpredictingmechanicalpropertiesoffiberreinforcedconcrete[J].InternationalJournalofSolidsandStructures, 2015, 56/57: 280-289.

[11]YangZihao,CuiJunzhi.Thestatisticalsecond-ordertwo-scaleanalysisfordynamicthermo-mechanicalperformancesofthecompositestructurewithconsistentrandomdistributionofparticles[J].ComputationalMaterialsScience, 2013, 69: 359-373.

[12]LiYouyun,CuiJunzhi.Themulti-scalecomputationalme-thodforthemechanicsparametersofthematerialswithrandomdistributionofmulti-scalegrains[J].CompositesScienceandTechnology, 2005, 65(9): 1447-1458.

[13] 姚振漢, 王海濤. 邊界元法[M]. 北京: 高等教育出版社, 2010.YaoZhenhan,WangHaitao.Boundaryelementmethod[M].Beijing:HigherEducationPress, 2010.(inChinese)

[14] 高效偉, 王靜, 彭海峰, 等. 高等邊界單元法——理論與程序[M]. 北京: 科學(xué)出版社, 2015.GaoXiaowei,WangJing,PengHaifeng,etal.Advancedboundaryelementmethod:theoryandprogramm[M].Beijing:SciencePress, 2015.(inChinese)

[15]YuYan,CuiJunzhi,HanFei.Aneffectivecomputergenerationmethodforthecompositeswithrandomdistributionoflargenumbersofheterogeneousgrains[J].CompositesScienceandTechnology, 2008, 68(12): 2543-2550.

[16] 韓非. 隨機復(fù)合材料力學(xué)性能預(yù)測的二階雙尺度方法[D]. 西安: 西北工業(yè)大學(xué), 2010.HanFei.Thesecond-ordertwo-scalemethodforpredictingmechanicalperformanceofrandomcompositematerials[D].Xi’an:NorthwesternPolytechnicalUniversity, 2010.(inChinese)

(編輯：馬文靜)

Research on Statistical Multiscale Boundary Element Algorithm for Random Composite Materials

Li Yiqiang1, Yang Zihao1, Yang Zhiqiang2

(1.School of Natural and Applied Sciences, Northwestern Polytechnical University, Xi’an 710072, China)(2.School of Astronautics, Harbin Institute of Technology, Harbin 150001, China)

Random uncertainty is the inherent attribute of natural and engineering materials; the heat conduction behavior is one of the most important problems in servicing process of thermal protection system, and the research on the heat conduction problem of random composite materials can provide theoretical basis for integrative design of material and structure in thermal protection system. A novel statistical multiscale analysis method based on two-scale asymptotic expansions is proposed to predict heat conduction performances of random porous materials. The statistical multiscale formulations and statistical multiscale boundary element algorithm are brought forward. Besides, the validity of the proposed method by comparison with theoretical results is verified. Finally, the effect of microstructures on the macroscopic thermal properties for the porous ceramic materials is investigated. Numerical results prove the accuracy and efficiency of our method for multiscale simulation of heat conduction problem in random porous materials.

random composite materials; heat conduction performance; statistical multiscale boundary element algorithm; effective conduction parameter

2017-02-24；

2017-03-29

國家自然科學(xué)基金(11471262,11501449)

李義強，liyiqiang@mail.nwpu.edu.cn

1674-8190(2017)02-199-07

TB33

A

10.16615/j.cnki.1674-8190.2017.02.012

李義強(1984-)，男，博士研究生。主要研究方向：邊界元方法、并行計算。

楊自豪(1987-)，男，博士，講師。主要研究方向：復(fù)合材料多尺度分析方法。

楊志強(1984-)，男，博士，副教授。主要研究方向：復(fù)合材料多尺度分析方法。