高填充細觀模型不確定性及其表征*

姚 東，高 波，楊月誠，史宏斌

(1第二炮兵工程大學，西安 710025;2中國航天科技集團公司第四研究院第41所，國體火箭發動機燃燒、熱結構與內流場國防科技重點實驗室，西安 710025)

0 引言

針對顆粒填充/增強(韌)復合材料的細觀模型，高希光[1]、宋玉普[2]、邵軍超[3]等分別在陶瓷基復合材料、混凝土材料、顆粒增強金屬基復合材料上開展了大量工作。宋衛東[4]研究了圓形顆粒填充的分子動力學方法，趙玖玲[5]基于最小碰撞時間探討了橢球體形狀參數的控制。宋衛東[4]、趙玖玲[5]等采用SEM結合數字圖像處理技術分別對鎢合金、混凝土及固體推進劑開展了細觀模型的研究。侯善芹[6]等在雙尺度分析的框架下對一定條件下顆粒數目、胞元尺寸等的影響進行了研究。姜芳[7]、Brassart[8]等基于Mori-Tanaka理論和Eshelby等效夾雜理論考慮了顆粒與基體之間粘接界面的脫粘，魏高峰[9]等則借助有限覆蓋技術和改進的Voronoi單元對同類問題進行了研究。

從現有文獻報道來看，建立細觀模型一般有近似法(順序法、并序法等)以及直接法兩大類。對于近似方法，如果不考慮顆粒填充過程中重疊檢測、取向等細節，填充、構造的最終效果均可以等效為對既有顆粒樣本的放置。很明顯，細觀結構存在以下兩個層次的不確定性:滿足指定規則時，顆粒總數不確定;滿足指定規則且顆粒總數一定時，顆粒分布、取向等不確定。對于具體的計算過程而言，還存在計算胞元的尺寸問題。

侯善芹[6]針對30%含量顆粒材料開展的均勻化分析和數值研究表明，顆粒隨機分布時宏觀等效模量是正態分布的，且所取體胞尺寸與顆粒尺寸之比越大，模量越集中在樣本均值的小領域內;并指出胞元尺寸L存在與體積分數有關的臨界值，大于該值后初始模量趨于穩定。Drugan針對大小為R的顆粒的研究表明，該臨界值為 L/R=5。趙玖玲[5]采用雙尺度方法針對復合固體推進劑的研究則表明，對于給定的體胞，填充顆粒越多，正態分布的初始模量概率密度的帶寬越窄。總體來看，相比于均勻化等宏細觀統一分析方法的全面工作，細觀結構不確定性的系統研究尚未見報道。

文中以高填充分數的固體推進劑為背景，研究了球形顆粒填充細觀結構的約束不完整性，提出了顆粒空間分布差異的子域形心表征方法。

1 體積填充/質量分數及粒徑分布的約束不完整性

設單相增強顆粒半徑服從正態分布，即r～N(μ，σ2)，則體積填充或質量分數的要求為:

隨機樣本的二重性表明Vf是隨機變量。計算了不同粒徑規格下顆粒總體積vp與顆粒總數目np。vp可擬合為np的線性函數，但均存在局部的“漲落”。

1)相同均值μ下，滿足相同體積分數所需的顆粒數(擬合值)隨σ的增大而減小。設滿足Vf=vp/v的np線性擬合結果為，對于 v=1、μ 為 0.1、σ 分別為0.01/0.05的情況，滿足體積分數0.6所需分別約140/75。

圖1 線性擬合與局部“漲落”

2)局部“漲落”導致滿足一定體積分數的顆粒數目不唯一。以Vf=vp/v=0.60為例，考慮有效位數取舍，r ～ N(0.01，0.0012)時顆粒數目從 138000 ～140000均可滿足，r～ N(0.1，0.052)時顆粒數目則在60～90之間。

圖2 整體線性關系與局部“漲落”

2 顆粒空間分布不確定性及其表征

設細觀構型I、II存在圖2所示的局部差異:I、II中填充顆粒數目、尺寸均一致，除了微元dΩ、dΩ'在I、II中分布位置有差異外，其余顆粒的分布、微元內部的顆粒分布、兩微元的形狀均一致。設 I、II參考坐標重合后dΩ、dΩ'形心之間的矢量為D。

圖3 細觀構型的局部差異

如不考慮熱應力等其它因素，上述差異可以通過dΩ與dΩ'的位置互換來實現。在dΩ或dΩ'足夠小時，構型I、II的應力/應變僅在其相應位置及附近領域內存在差異，從而對整體的性能產生一定的、有限的影響。

填充顆粒數目及粒徑均一致、填充分數較高時，構型I、II結構差異的直接描述是非常困難的。考慮到質量的分布是構型的直觀反映，文中發展了子域形心表示方法。

設存在獨立于胞元的網格將胞元劃分為n個互不重疊的、單聯通的子域，m個顆粒的球心落入子域Φi內，設其形心設為r(i)，則有:

式中:ljk為Φi內顆粒k的形心在整體坐標軸j方向的分量，vt為顆粒的體積(t=1，2，…，m)。當各子域均只包含1個顆粒時，r(i)即為描述顆粒位置的坐標。

3 Matlab/Abaqus統一的分析平臺及算例

3.1 Matlab/Abaqus統一的分析平臺

由于增強顆粒的形貌、大小、取向等均存在隨機性，相比于三維四向、針刺4D編織等結構幾乎固定不變的三維 /二維編織復合材料，運用漸進展開與攝動求解顆粒增強復合材料的數學要求較高且效率較低。能量法在特定邊界條件下探討胞元變形能與均質等效體變性能的關系，結合有限元技術建立了簡潔、高效的等效性能預測方法。

在Matlab平臺下開發了順序法填充生成細觀模型的算法，集成了約束強化、顆粒生成、三維填充、子域劃分 /斷層掃描等功能;Python腳本文件導入顆粒尺寸及填充位置數據文件，生成Abaqus的.cae文件并提取預制集合的輸出信息;返回Matlab平臺下按能量等效對等效初始模量進行求解。

3.2 算例——顆粒空間分布不確定性對初始模量的影響分析

圖4 分析平臺的結構

初始模量 E 呈正態分布，μE為 19.73EB，σE為0.20EB，EB為基體模量。取背景網格為2 ×2 ×2(分胞元為8份)，[rij]規模為8 × 3。對 E 分別為19.96EB、19.62EB的模型2、模型19，均按式(3)處理，圖像對比見圖6。

圖5 拉伸方向的應變云圖

圖6 差異對比

4 結論

文中以固體推進劑為背景，研究了顆粒填充細觀結構的不確定性，提出了兩個層次不確定性的表征方法。算例表明:

1)對于粒徑服從正態分布的球形顆粒，填充分數是不完整的約束，隨機樣本的二重性決定Vf是顆粒數目的隨機變量，顆粒總體積vp可擬合為顆粒總數目np的線性函數，但存在局部的“漲落”。對于固體推進劑，由于其顆粒組分體積分數高、粒徑小，np的不確定性更為明顯。

2)在給定顆粒數目的強化約束條件下，子域形心法結合圖像表示可以有效表征顆粒空間分布的不確定性，隨機生成的模型會導致初始等效模量的正態分布。

[1]高希光，宋迎東，孫志剛.陶瓷基復合材料高精度宏細觀統一本構模型研究[J].航空動力學報，2008，23(9):1617－1622.

[2]宋玉普.多種混凝土材料的本構關系和破壞準則[M].北京:中國水利水電出版社，2002.

[3]邵軍超，劉越.顆粒增強金屬基復合材料力學行為有限元模擬研究現狀[J].材料導報，2007，21(9):111－115.

[4]宋衛東，王靜，劉海燕.顆粒增強復合材料真實結構有限元建模[J].北京理工大學學報，2009，29(6):501－505.

[5]趙玖玲.復合固體推進劑力學損傷與細觀損傷演化過程數值模擬[D].西安:第二炮兵工程學院，2010.

[6]侯善芹，劉書田.顆粒增強復合材料彈性性能的統計特征分析[J].科學技術與工程，2008，18(15):4078－4082.

[7]姜芳，寧建國.有界面脫粘時顆粒增強金屬基復合材料的彈塑性性能分析[J].材料工程，2006(z1):366－373.

[8]Brassart L，Inglis H M，Delannay L，et al. An extended Mori-Tanaka method homogenization scheme for finite strain modeling of debonding in particle-reinforced elastomers[J]. Computational Materials Science，2009，45(3):611－616.

[9]魏高峰，馮偉，高洪芬.基于位移插值的Voronoi單元有限元方法[J].應用力學學報，2008，25(2):342－346.