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

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

0 引言

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

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

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

文中以高填充分數(shù)的固體推進劑為背景，研究了球形顆粒填充細觀結(jié)構(gòu)的約束不完整性，提出了顆?？臻g分布差異的子域形心表征方法。

1 體積填充/質(zhì)量分數(shù)及粒徑分布的約束不完整性

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

隨機樣本的二重性表明Vf是隨機變量。計算了不同粒徑規(guī)格下顆?？傮w積vp與顆?？倲?shù)目np。vp可擬合為np的線性函數(shù)，但均存在局部的“漲落”。

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

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

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

圖2 整體線性關(guān)系與局部“漲落”

2 顆?？臻g分布不確定性及其表征

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

圖3 細觀構(gòu)型的局部差異

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

填充顆粒數(shù)目及粒徑均一致、填充分數(shù)較高時，構(gòu)型I、II結(jié)構(gòu)差異的直接描述是非常困難的?？紤]到質(zhì)量的分布是構(gòu)型的直觀反映，文中發(fā)展了子域形心表示方法。

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

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

3 Matlab/Abaqus統(tǒng)一的分析平臺及算例

3.1 Matlab/Abaqus統(tǒng)一的分析平臺

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

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

3.2 算例——顆?？臻g分布不確定性對初始模量的影響分析

圖4 分析平臺的結(jié)構(gòu)

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

圖5 拉伸方向的應(yīng)變云圖

圖6 差異對比

4 結(jié)論

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

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

2)在給定顆粒數(shù)目的強化約束條件下，子域形心法結(jié)合圖像表示可以有效表征顆?？臻g分布的不確定性，隨機生成的模型會導(dǎo)致初始等效模量的正態(tài)分布。

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

[2]宋玉普.多種混凝土材料的本構(gòu)關(guān)系和破壞準(zhǔn)則[M].北京:中國水利水電出版社，2002.

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

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

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

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

[7]姜芳，寧建國.有界面脫粘時顆粒增強金屬基復(fù)合材料的彈塑性性能分析[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].應(yīng)用力學(xué)學(xué)報，2008，25(2):342－346.