鐘汝能，鄭勤紅，，向 泰，姚 斌

(1.云南師范大學能源與環(huán)境科學學院，云南昆明，650500；2.云南省光電信息技術重點實驗室，云南昆明，650500)

1 引言

在材料科學、電子電氣和電磁工程領域，常常遇到大量基質顆粒隨機分布于基體中的混合物質，如金屬基復合材料、聚合物基復合材料、含有雜質的礦物質或自然堆放的粒狀農產(chǎn)品等等，研究這類混合物的介電性能具有較好的價值[1]，然而，技術成本高一直是新材料研發(fā)或物料微波熱加工器件設計中的瓶頸問題，建立有效的數(shù)值仿真模型先行進行性能預測已被廣大研究者采用[2，3]。

眾所周知，混合物的等效介電性能依賴于各組成成分的介電常數(shù)、比例、形狀和位置分布等，且基質顆粒在基體中的分布都是非均勻體系，其形狀、大小、位置的存在狀態(tài)并非完全規(guī)則[4]。對此，學者們采用不同的技術方法開展了大量的模擬研究[5，6]，其中，Cheng采用ANSYS軟件對顆粒無序填充兩相復合材料的有效介電常數(shù)進行了數(shù)值模擬[7]；Rakesh等采用COMSOL多物理場仿真軟件(COMSOL Multiphysics)模擬分析了金屬核殼納米/聚合物介質復合電容器的有效介電性能[8]；Zhao等采用蒙特卡羅法和有限元法模擬計算了立方晶格三維結構中雙組分復合材料的復介電常數(shù)[4]；Jebbor使用有限元法和MATLAB軟件仿真了二維情形下球形顆粒隨機分布的復合材料介電性能[9]；浦毅杰基于蒙特卡羅方法，利用ANSYS的APDL參數(shù)化設計語言，構建了隨機分布納米顆粒增強陶瓷基復合材料性能數(shù)值分析模型[10]。前人的研究成果為工程的實際應用提供了借鑒，然而，在早期的研究中還面臨著一些可以改進的問題：①大部分研究主要針對雙組分混合物而開展，且在建立模型時或多或少假設了基質顆粒在基體中的分布是有序的重復單元，即顆粒的分布狀態(tài)主要基于結構抽象化和顆粒規(guī)則化兩種類型，顆粒的形狀多以球體、立方體、橢球或層狀等單一結構形狀為主[11]，這些簡化后的計算模型與混合體中顆粒的真實存在狀態(tài)并非完全吻合；部分學者開展了顆粒無序分布對混合物性能的影響研究[12]，但研究內容側重于顆粒的結構形狀或者分布狀態(tài)中的某一方面。②當前，針對多組分混合物和變溫變頻條件下混合物的介電性能模擬分析研究成果鮮見報道，在常用有限元商用軟件中，COMSOL的多物理場耦合能力使變溫變頻條件下混合物的性能分析成為可能，然而，該軟件中針對三維情形下顆粒隨機生成及投放的功能缺失，需用戶自行編制腳本程序來實現(xiàn)。

本文采用COMSOL的App開發(fā)器，編制了基于蒙托卡羅方法的顆粒隨機生成和投放程序，構建了基質顆粒的空間位置、形狀和體積在基體中隨機分布的立方體多胞有限元模型，通過對數(shù)值結果的比較分析驗證了方法的可行性和有效性，以期為混合物性能的深入模擬分析提供一種借鑒方法。

2 仿真原理和模型

2.1 顆粒填充混合物等效介電特性的仿真原理

當混合物中顆粒的尺寸遠遠小于電磁波波長時，可用準靜態(tài)條件下的平行板電容器模型分析混合物的等效介電特性[6]。如圖1所示，平行板電容器上、下極板的電勢分別為φ1=1V、φ2=0V，側面邊界條件為?φ/?n=0，極板間為混合物。

圖1 顆粒填充雙組分混合物模擬模型

在有限元方法中，每個網(wǎng)格單元所擁有的靜電能為

(1)

式中，εk和Vk分別表示第k個單元的介電特性和體積?；旌象w所諸存的靜電能為

(2)

由電磁理論可知，圖1所示電容器中所存儲的靜電能可等效地表示為[13]

(3)

式中，φ1-φ2為兩個極板之間的電勢差，d為兩個極板之間的距離，S為極板的面積，ε0=8.85×10-12F/m，εeff為混合物等效介電特性。

令Wc=We，由式(2)、(3)即可獲得圖1所示電容器中混合物的等效介電特性(εeff)。

在分析微波頻段下導電型混合物的介電特性時，物質的介電特性εk=1-j(σ/ωε0)，其中σ為電導率，ω為角頻率(ω=2πF)，F(xiàn)為微波頻率。

為準確模擬計算各向同性混合物的介電特性，可生成同一條件下的多個模擬模型，取其所有模型計算結果的平均值作為數(shù)值結果，即:

(4)

2.2 隨機顆粒填充混合物的模型構建

2.2.1 模型構建思路

設定一個大的立方體(或長方體)代表基體，假設立方基體中可以劃分出很多相似的小立方單元，每個單元由小立方基體和中心在立方對稱中心附近、邊界不超出單元邊界的介質顆粒構成。與傳統(tǒng)模型不同的是該單元不再是重復單元，即每個單元中顆粒的位置、大小和形狀是隨機分布的；對于同一個顆粒，它可能隨機出現(xiàn)在所劃分出的任意一個小立方單元中，顆粒的中心也沒有固定在小立方對稱中心，而是結合顆粒的結構參數(shù)隨機分布在小立方對稱中心附近。

2.2.2 顆粒生成及投放步驟

1)在COMSOL軟件的“App開發(fā)器”中創(chuàng)建“Model method”，在“聲明”中定義不同的數(shù)組、變量分別儲存各顆粒的位置坐標值、體積值、體積比和結構參數(shù)。

2)將基體劃分為M×N×P個相似的小單元，M、N、P根據(jù)混合物的基體結構進行設置，文中選取M=N=P=10。

3)根據(jù)實際需求，使用“錄制”方法生成不同結構形狀顆粒的建模代碼備用，根據(jù)需要將代碼嵌入到主程序中，設置結構參數(shù)變化范圍，使用隨機函數(shù)(Math.random)分步生成顆粒的隨機位置坐標參數(shù)和隨機結構參數(shù)。

4)將每組隨機數(shù)與小立方單元中的填充顆粒建立對應關系，將隨機值賦予擬生成顆粒，并通過程序判斷語句，確保顆粒不溢出小單元、顆粒間可以接觸但不重疊。

5)設定顆粒數(shù)量，編制程序將介電性能相同的顆粒設定為同一個“累積”，累積的數(shù)量根據(jù)混合物中顆粒的種類確定。

6)運行method方法，實現(xiàn)顆粒的隨機生成及投放。

2.2.3 有限元模型構建和賦值

執(zhí)行method程序，在指定區(qū)域生成立方體模型和顆粒模型，使用布爾運算將立方體區(qū)域減去顆粒區(qū)域即得到基體區(qū)域，形成“聯(lián)合體”，完成模型構建。作為示例，圖2顯示了不同形狀顆粒填充混合物的三維模型圖，其中圖2(A)為球狀顆粒模型圖，圖2(B)為球狀和立方狀顆?；旌夏Ｐ蛨D，圖2(C)為球狀、立方狀和橢球狀顆?；旌夏Ｐ湍Ｐ蛨D。

圖2 不同形狀顆粒填充混合物的三維模型圖、顆粒網(wǎng)格剖分示意圖

使用COMSOL網(wǎng)格剖分功能中的“物理場控制網(wǎng)格”對模型進行自由四面體網(wǎng)格剖分，并采用“中等”優(yōu)化級別進行單元質量優(yōu)化，圖2(D)顯示了圖2(B)所示模型中顆粒網(wǎng)格的剖分情況。在材料屬性中選取基體區(qū)域進行屬性賦值，選擇“累積”對不同的基質顆粒進行物理屬性賦值。

可以看出，通過本程序可以實現(xiàn)任意不同結構形狀的顆粒在不同體積分數(shù)下的隨機化混合生成和投放，該模型更符合實際工程中基質顆粒物質在基體中的真實狀態(tài)，且能實現(xiàn)多種組成成分混合物的模擬模型構建。

3 結果和討論

3.1 隨機模型的有效性

采用球狀顆粒數(shù)值計算得到ε1/ε2為1/10時混合物等效介電常數(shù)εeff與體積分數(shù)fv的關系(見圖3)，每個體積分數(shù)計算50個拓撲結構后取其平均值作為最終結果。將計算結果與Maxwell Garnet (MG)、Bruggeman、JayaSundere Smith (JS)等經(jīng)典模型[14]進行比較，結果表明，體積分數(shù)計算范圍fv∈(0，50%)的數(shù)值結果介于JS與Bruggeman之間，基質顆粒的體積分數(shù)較小時(fv<20%)四種模型結果比較接近，基質顆粒的體積分數(shù)較大時(fv>20%)數(shù)值計算結果介于JS模型和Bruggeman模型之間，但與MG模型相差較大。對于MG、Bruggeman和JS模型來說，MG主要適用于基質的體積分數(shù)相對較小時的球狀顆粒情形，體積分數(shù)增大后會逐漸形成拓撲微結構，此時計算值會偏?。籅ruggeman模型在計算兩相無明顯包圍的對稱有效介質時能獲得較好的結果[15]；而JS模型考慮了混合物中各組成成分之間相互作用。所以數(shù)值計算結果在基質的體積分數(shù)較小時與3種經(jīng)典模型吻合較好，而在體積分數(shù)較高時偏離MG模型，與Bruggeman模型和JS模型接近，這種變化趨勢是正確的，說明所建立的隨機模型是有效的。

圖3 隨機模型數(shù)值結果與經(jīng)典理論對比情

圖4 隨機模型數(shù)值結果的分布情況

3.2 數(shù)值方法的有效性

3.3 顆粒結構形狀對混合物介電性能的影響

以聚偏氟乙烯雙組分復合材料為對象(ε1=1、ε2=2.61+0.12j)，采用隨機模型計算了顆粒為不同結構形狀時混合體的等效介電常數(shù)，其中，A、B、C、D、E分別表示基質顆粒為球狀、立方狀、長方狀、圓柱狀和橢球狀。圖5(a)結果表明，當顆粒的體積分數(shù)較小時(fv<20%)，顆粒結構形狀對等效介電特性的影響不明顯；當顆粒的體積分數(shù)相對較大時(fv>20%)，等效介電常數(shù)大小與顆粒結構形狀的對應關系為：εD>εA+B>εB>εA+C>εC>εA+D>εA>εD+E>εA+E>εE。分析可得，顆粒結構形狀為D(圓柱狀)和A+B(球狀+立方狀)時計算結果相對較大，此時數(shù)值結果與實測值[17]更加接近；顆粒結構為E(橢球狀)時計算結果相對較小，這是由于橢球的特殊規(guī)則結構導致了粒子內偶極矩方向稍稍偏離[18]；其余結構形狀對應的等效介電常數(shù)介于εD和εE之間。為便于比較，圖5(b)顯示了ε1/ε2=1/10時4種不同結構形狀下εeff與fi對應關系，圖示表明，當基質顆粒和基體的介電常數(shù)比相對較大時，上述趨勢更加明顯。

圖5 顆粒結構形狀對介電常數(shù)的影響

3.4 雙組分混合物的等效介電性能分析

以聚偏氟乙烯雙組分復合材料為對象(ε1=1、ε2=2.61+0.12j)，采用球狀顆粒計算得到了混合體等效介電常數(shù)εeff與填充比fv的關系如圖6所示，將數(shù)值結果與顆粒規(guī)則分布、位置隨機分布和實測值[17]進行了對比。結果表明，由于隨機模型考慮了顆粒的空間位置和體積的隨機性，與顆粒為簡立方分布、面心立方分布、位置隨機分布時相比，數(shù)值結果更接近實測值，說明該隨機模型具有更好的準確性。但精確地講，無論是有序分布還是隨機分布，計算結果比實測值都小一些，這可能與模擬顆粒的結構形狀單一、顆粒的高階多極矩貢獻及多重散射效應引起的附加損耗有關[19]。

圖6 采用隨機模型計算的介電常數(shù)數(shù)值結果與文獻值的比較(聚偏氟乙烯復合材料)

為評價隨機模型在高介電常數(shù)比情況時的準確性，以BaTiO3/PVDF復合材料為對象(ε′BaTiO3=2000、ε′PVDF=9.6)，采用球狀顆粒和立方顆粒計算得到體積分數(shù)為5%、10%、15%、20%時，得到復合材料等效介電常數(shù)分別為11.92、14.44、16.57、19.78，計算結果與實驗值相吻合[20]。

3.5 多組分混合物的等效介電性能分析

以硅基氮化物陶瓷復合材料(SiO2-AIN-BN)為對象，使用隨機模型計算三組分復合材料的等效介電常數(shù)，基體為SiO2(εSiO2=4.2-0.002j)，設定球狀顆粒為第一基質相AIN(εAIN=8.8-0.076j)，立方狀顆粒為第二基質相BN(εBN=5-0.040j)。選取各組分的體積分數(shù)比fSiO2:fAIN:fBN=7:1:2，得到SiO2-AIN-BN復合材料的等效介電常數(shù)為εeff=4.7286-0.0267j，計算結果與實測值[21]相吻合。

4 結論

1)基于蒙托卡羅法的顆粒隨機生成、投放程序，可以實現(xiàn)基質顆粒的空間位置和體積在基體中的隨機分布，隨機模型更符合現(xiàn)實工程材料中顆粒的真實存在狀態(tài)，模型適用于二維、三維情形下混合體的介電性能分析。

2)基于能量法原理和有限元理論的數(shù)值計算結果對比分析表明，基質顆粒的空間位置和結構形狀對數(shù)值結果有影響，選擇長徑比相對較大的圓柱狀顆粒(D)或球狀、立方狀組合顆粒群(A+B)進行數(shù)值仿真計算，能獲得與實測值更加接近的數(shù)值結果。

3)在隨機模型中，將不同的基質顆粒設定為相應的結構形狀和“累積”，可自由實現(xiàn)多種組分的顆粒生成和介電參數(shù)賦值并進行仿真計算，這為多組分(各組分的形狀各異)混合體的介電性能研究提供一種可借鑒方法。

4)在隨機模型的基礎上，充分利用COMSOL的多物理場耦合功能，可對變溫變頻條件下顆粒填充混合物的傳熱傳質過程進行模擬仿真分析，相關研究將在下一步工作中開展。