隨機復合材料統計多尺度邊界元算法研究

李義強，楊自豪，楊志強

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

隨機復合材料統計多尺度邊界元算法研究

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

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

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

隨機復合材料；熱傳導性能；統計多尺度邊界元算法；等效導熱系數

0 引 言

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

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

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

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

1 隨機復合材料微結構表征

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

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

(1)

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

(2)

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

2 統計多尺度分析模型

對于隨機分布復合材料，考慮如下具有混合邊界條件的穩態熱傳導問題：

(3)

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

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

(4)

由于ξ=x/ε，存在如下鏈式法則：

(5)

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

O(ε)=f

(6)

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

(7)

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

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

(8)

通過式(9)可得與樣本ωs有關的均勻化熱傳導系數：

(9)

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

(10)

(11)

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

(12)

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

(13)

3 統計多尺度邊界元算法

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

3.1 單胞問題邊界元算法

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

(14)

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

(15)

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

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

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

(τ∈?Qs)

(16)

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

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

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

(τ∈?Qs)

(17)

其中，

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

(18)

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

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

(19)

(20)

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

再采用徑向積分方法將上述區域積分轉換為邊界積分，即

(21)

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

3.2 算法流程

基于統計多尺度邊界元算法的隨機復合材料熱傳導性能的算法流程為：

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

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

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

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

4 算例分析

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

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

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

表1 圓球孔洞均勻分布不同體積分數時均勻化導熱系數的多尺度解與HS界的比較

表2 圓球孔洞正態分布不同體積分數時均勻化導熱系數的多尺度解與HS界的比較

表3 橢球孔洞傾角正態分布不同體積分數時均勻化導熱系數的多尺度解與HS界的比較

從表1～表3可以看出：對于不同體積分數和分布狀態的多孔隨機分布復合材料，其多尺度數值解均在HS界內，且對于相同分布狀態的復合材料結構，其宏觀均勻化系數隨孔洞體積分數的增大而發生較大變化，表明孔洞的體積分數、結構形態以及分布狀態對復合材料的宏觀等效導熱性能具有重要影響。

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

表4 有限元網格與邊界元網格數比較

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

表5 圓球孔洞均勻分布時均勻化導熱系數計算結果

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

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

5 結 論

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

(2) 多孔隨機復合材料的熱傳導性能不僅依賴于宏觀條件，還與復合材料的細觀結構有關，包括孔洞的體積分數、結構形態以及分布狀態等，統計多尺度邊界元算法能夠較為準確地捕捉多孔隨機復合材料的細觀結構信息。

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

[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]. 計算數學, 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]. 北京: 科學出版社, 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] 韓非. 隨機復合材料力學性能預測的二階雙尺度方法[D]. 西安: 西北工業大學, 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

國家自然科學基金(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-)，男，博士，講師。主要研究方向：復合材料多尺度分析方法。

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