李慶，楊曉翔

(福州大學(xué)機(jī)械工程及自動(dòng)化學(xué)院，福建福州 350116)

0 引言

橡膠是一種無(wú)定形的高聚物，又可稱(chēng)為彈性體.橡膠具有優(yōu)異的彈性、良好的柔軟性，較高的耐磨性、電絕緣性、致密性以及耐腐蝕、耐溶劑、耐高低溫等特殊性能，廣泛應(yīng)用于能源建設(shè)、交通運(yùn)輸、土木建筑、醫(yī)療衛(wèi)生、軍事國(guó)防以及航天、航空、宇宙開(kāi)發(fā)等多個(gè)領(lǐng)域.雖然橡膠具有優(yōu)良的特性，但是單用生膠(包括天然橡膠及合成橡膠)并不能制得符合各種使用要求的橡膠制品.工業(yè)上使用的橡膠材料都是由本體材料和各種配合劑構(gòu)成的細(xì)觀復(fù)合材料，其中炭黑填料粒子是橡膠工業(yè)中最重要的補(bǔ)強(qiáng)劑，它能提高橡膠的模量、拉伸強(qiáng)度、撕裂強(qiáng)度、耐磨性、扯斷伸長(zhǎng)率及黏度等，會(huì)降低壓縮永久變形和口型膨脹等性能，所以，炭黑的填充明顯改善了橡膠的靜態(tài)與動(dòng)態(tài)性能［1］.

炭黑填料對(duì)改善橡膠材料的力學(xué)性能起著至關(guān)重要的作用，它對(duì)橡膠性能的影響主要取決于填料的粒徑、含量、形態(tài)結(jié)構(gòu)和表面性質(zhì)等結(jié)構(gòu)因素［2］，所以，為指導(dǎo)相關(guān)材料設(shè)計(jì)，縮短研發(fā)周期，必須建立炭黑顆粒填充橡膠復(fù)合材料的細(xì)觀結(jié)構(gòu)與宏觀性能的內(nèi)在聯(lián)系和規(guī)律的定量關(guān)系.對(duì)此，Govindjee［3］，Bergstrom＆Boyce［4］以及夏勇［5］等國(guó)內(nèi)外很多學(xué)者進(jìn)行了探索研究，并建立多種細(xì)觀數(shù)值模型對(duì)炭黑填充橡膠進(jìn)行力學(xué)行為分析，主要有單顆粒夾雜二維平面應(yīng)力模型和三維軸對(duì)稱(chēng)模型，多顆粒夾雜二維平面應(yīng)力模型以及采用若干個(gè)相同的菱形十二面體建立的三維數(shù)值模型.就炭黑填充橡膠復(fù)合材料細(xì)觀數(shù)值模型的發(fā)展而言，目前僅處在初步階段，實(shí)際考察的因素僅有體積分?jǐn)?shù)一種，仍有許多不完善之處.

本研究基于細(xì)觀結(jié)構(gòu)的周期性假設(shè)，建立聯(lián)系復(fù)合材料宏觀和細(xì)觀特征的細(xì)觀數(shù)值模型，對(duì)炭黑填充橡膠復(fù)合材料進(jìn)行單向拉伸的數(shù)值模擬，比較分析圓形和方形炭黑填料粒子模型的變形場(chǎng)和應(yīng)力場(chǎng)，以及炭黑填料的含量對(duì)于復(fù)合材料有效彈性模量的影響規(guī)律.

1 復(fù)合材料細(xì)觀力學(xué)研究基礎(chǔ)

1.1 代表性體積單元

如圖1所示，對(duì)于顆粒增強(qiáng)型復(fù)合材料通常假定宏觀材料具有周期性的細(xì)觀結(jié)構(gòu)，即采用代表性體積單元(RVE)進(jìn)行細(xì)觀力學(xué)研究，圖1(b)為單顆粒夾雜RVE示意圖，圖1(c)為多顆粒夾雜RVE示意圖.對(duì)于RVE的選取，它應(yīng)該選取足夠大以包含充分的細(xì)觀結(jié)構(gòu)組成信息，而同時(shí)RVE的尺寸又要足夠小以反映其細(xì)觀非均質(zhì)性［6－7］.

圖1 RVE選取示意圖Fig.1 Extracting diagram of RVE

1.2 周期性邊界條件

代表性體積單元是宏觀復(fù)合材料連續(xù)統(tǒng)一體的代表，需要在鄰近代表性體積單元的邊界施加周期性邊界條件［5，8－9］以確保變形場(chǎng)的協(xié)調(diào)性，主要包括位移連續(xù)條件和應(yīng)力連續(xù)條件.

1.2.1 位移連續(xù)條件:RVE對(duì)邊邊界形狀保持相同

圖2所示為二維代表性體積單元的幾何表述，其外框?yàn)橐贿呴L(zhǎng)為L(zhǎng)的正方形，頂點(diǎn)A0為坐標(biāo)原點(diǎn)，頂點(diǎn)A1和A2分別位于坐標(biāo)軸x和y上，則對(duì)邊邊界上點(diǎn)的自由度應(yīng)滿足以下約束關(guān)系:

圖2 二維RVE的幾何表述Fig.2 Geometrical frame of two －dimensional RVE

式中:uAi代表頂點(diǎn)Ai的位移矢量，i=0，1，2.若分別用l1、l2、l3、l4代表RVE中x=0、x=L、y=0、y=L時(shí)的邊界，用uli代表RVE邊界li上點(diǎn)的位移矢量，則約束方程(1)和(2)可寫(xiě)為:

則RVE邊界l2和l4上的節(jié)點(diǎn)(不包含頂點(diǎn))稱(chēng)為被約束節(jié)點(diǎn)，而邊界l1和l3上的節(jié)點(diǎn)稱(chēng)為保留節(jié)點(diǎn)，被約束節(jié)點(diǎn)的位置由保留節(jié)點(diǎn)的位置來(lái)確定.

由邊界虛位移引起的虛應(yīng)變能變化為:

式中:σl代表RVE邊界的柯西應(yīng)力張量，nl代表其外法線單位向量.由lRVE=l1+l2+l3+l4，則式(5)可展開(kāi)為:

把位移約束方程式(3)和式(4)代入式(6)，可得:

因?yàn)镽VE對(duì)邊邊界形狀相同，l1=l2，l3=l4，則式(7)可寫(xiě)為:

根據(jù)變形體的虛功原理，外力在虛位移上所做的外力虛功dW恒等于內(nèi)應(yīng)力在虛應(yīng)變上所做的虛變形功(即虛應(yīng)變能變化dV)，則有:

式中:Pli代表外部作用于RVE邊界li上的分布載荷，fAi代表外部作用于RVE頂點(diǎn)Ai上的集中載荷.根據(jù)式(8)和式(9)，可得出外部載荷表達(dá)式為:

1.2.2 應(yīng)力連續(xù)條件:RVE對(duì)邊邊界上應(yīng)力矢量大小相等，方向相反

對(duì)邊邊界上點(diǎn)的應(yīng)力應(yīng)滿足以下約束關(guān)系:

將式(17)和式(18)分別代入到式(15)和式(16)，可得:

由以上推導(dǎo)可得，周期性邊界條件約束的結(jié)果不僅使被約束節(jié)點(diǎn)的位置可以由保留節(jié)點(diǎn)的位置來(lái)確定，還使RVE所受的外部載荷可以由分配到頂點(diǎn)上的集中載荷來(lái)等效，因此，具有周期性的二維RVE模型的外加載荷僅僅需要作用于頂點(diǎn)A0，A1和A2上.

2 有限元建模

形如圖1(b)所示的單顆粒夾雜代表性體積單元只適用于研究顆粒分布均勻的填充復(fù)合材料，對(duì)于顆粒隨機(jī)分布的復(fù)合材料，單顆粒夾雜RVE已經(jīng)不再適用，必須建立符合顆粒隨機(jī)分布真實(shí)情況的RVE.所以，本研究采用多顆粒夾雜的代表性體積單元，如圖1(c)所示，考慮復(fù)合材料中顆粒隨機(jī)分布特征，通過(guò)建立具有隨機(jī)分布形態(tài)的代表性體積單元的數(shù)值模型，來(lái)模擬分析炭黑顆粒填充橡膠復(fù)合材料的力學(xué)行為.

圖3 RVE模型Fig.3 RVE model

對(duì)于多顆粒夾雜隨機(jī)分布的二維代表性體積單元建模時(shí)，將橡膠基體和炭黑粒子之間的界面假設(shè)為理想連接狀態(tài)，不發(fā)生脫落現(xiàn)象［4］.如圖3所示，采用RSA算法［10］通過(guò)MATLAB隨機(jī)函數(shù)程序在給定區(qū)間［0，L］產(chǎn)生兩個(gè)均勻分布的偽隨機(jī)數(shù)，用來(lái)確定顆粒形心的位置坐標(biāo)，分別建立圓形粒子和方形粒子的RVE模型，RVE邊長(zhǎng)均為L(zhǎng)=24 μm.其中，圓形顆粒半徑為r=1 μm，顆粒個(gè)數(shù)m=10，顆粒所占體積分?jǐn)?shù)φr=5.45%;方形顆粒邊長(zhǎng)為b=2 μm，顆粒個(gè)數(shù)n=8，顆粒所占體積分?jǐn)?shù)φb=5.56%.

模型中炭黑顆?？醋鳛榫€性彈性材料，彈性模量Ef=200 MPa，泊松比設(shè)為 υf=0.3［4］.

橡膠基體本構(gòu)關(guān)系采用近乎不可壓縮的Mooney－Rivlin模型［11－12］

其中:C10和C01是同溫度有關(guān)的材料常數(shù)，與初始剪切模量μ0的關(guān)系為［12－13］:

設(shè)初始彈性模量Em=2 MPa，泊松比υm=0.5，則有初始彈性模量與剪切模量關(guān)系式［12］

根據(jù)式(21)和式(22)，并取材料常數(shù)C10與C01的比值為4［13］，可計(jì)算得到Mooney－Rivlin模型中的材料常數(shù)C10=0.266 7 MPa，C01=0.066 7 MPa.

在Abaqus/Standard中進(jìn)行網(wǎng)格劃分時(shí)，橡膠基體采用平面應(yīng)變8節(jié)點(diǎn)四邊形二次完全積分的雜交單元CPE8H①ABAQUS Inc.ABAQUS analysis user’s manual.The State of Rhode Island and Providence Plantations:ABAQUS Inc，1978.，炭黑顆粒采用CPE8.特別提出的是，要保證模型對(duì)邊網(wǎng)格相同，節(jié)點(diǎn)一一對(duì)應(yīng)，以便于施加周期性邊界條件的約束方程.在ABAQUS軟件中，通過(guò)添加約束方程*Equation來(lái)實(shí)現(xiàn)方程式(1)和方程式(2)所描述的約束.為給定RVE有限元模型充分的約束，還需要約束頂點(diǎn)A0在x、y方向的位移，約束頂點(diǎn)A1在y方向的位移，約束頂點(diǎn)A2在x方向的位移.然后通過(guò)在頂點(diǎn)A1施加x方向位移載荷進(jìn)行單軸拉伸的模擬.

3 結(jié)果與討論

3.1 變形場(chǎng)和應(yīng)力場(chǎng)分析

基于本研究提出的方法和參數(shù)對(duì)圖3所示的幾何模型沿x方向施加拉伸位移載荷，對(duì)RVE施加周期性邊界條件約束后進(jìn)行有限元計(jì)算.圖4給出了拉伸位移為15 μm時(shí)的變形圖，從圖4中可以看出橡膠基體產(chǎn)生的變形較大，而炭黑顆粒基本上沒(méi)有變形.通過(guò)分析發(fā)現(xiàn)RVE對(duì)邊邊界的形狀在變形過(guò)程中始終保持完全相同，符合周期性結(jié)構(gòu)位移連續(xù)的真實(shí)情況.

圖5給出了拉伸位移為15 μm時(shí)的von－Mises應(yīng)力等值線圖，并且沿x和y方向分別周期排列2個(gè).從兩幅圖中可以看出模型中應(yīng)力分布并不均勻，距離填充顆粒越近的橡膠基體處等值線越密集，應(yīng)力梯度越大;距離填充顆粒越遠(yuǎn)的橡膠基體處等值線越稀疏，應(yīng)力梯度越小.將方形粒子與圓形粒子模型的應(yīng)力等值線圖相比較可以看出，在相同的載荷條件下，方形粒子模型中的應(yīng)力集中現(xiàn)象更加明顯，應(yīng)力值相對(duì)更高，如圖5所示，方形粒子模型中最大應(yīng)力值為29.996 MPa，而圓形粒子模型中的最大應(yīng)力值只有5.066 MPa.分別觀察兩幅圖中相鄰RVE共同邊界處的應(yīng)力等值線分布情況，發(fā)現(xiàn)它們?cè)诩虞d過(guò)程中始終構(gòu)成連續(xù)閉合曲線，RVE中對(duì)邊邊界的應(yīng)力等值線分布相同，滿足RVE對(duì)應(yīng)邊界的對(duì)應(yīng)節(jié)點(diǎn)應(yīng)力一致條件，符合周期性結(jié)構(gòu)應(yīng)力連續(xù)的真實(shí)情況.

圖4 RVE的變形圖Fig.4 Deformation of RVE

圖5 RVE的von－Mises應(yīng)力等值線圖Fig.5 The distribution of von － Mises stress in RVE

圖6給出了圓形和方形粒子模型的平均應(yīng)力－應(yīng)變曲線，其中，復(fù)合材料的真實(shí)應(yīng)力(MPa)由加載方向的約束反力除以變形后的表面積得到;真實(shí)應(yīng)變由如下公式得到:

式中:u為加載的位移(μm)，L為RVE邊長(zhǎng)(μm).

從圖6中可以看出應(yīng)力與應(yīng)變呈非線性關(guān)系，符合橡膠材料的非線性彈性的特性.所以，本文研究所建立的模型可以用來(lái)模擬炭黑填充橡膠復(fù)合材料的力學(xué)行為.

3.2 炭黑顆粒含量和形狀對(duì)復(fù)合材料有效彈性模量的影響

圖6 RVE的平均應(yīng)力－應(yīng)變曲線Fig.6 The stress－ strain behavior of RVE

基于本研究提出的方法和參數(shù)，保持炭黑顆粒大小不變，增加顆粒數(shù)目，分別建立不同炭黑體積含量的代表性體積單元，通過(guò)平面應(yīng)變模型進(jìn)行單軸拉伸的有限元模擬，計(jì)算橡膠復(fù)合材料在不同炭黑含量下的平均應(yīng)力－應(yīng)變曲線，如圖7所示.

圖7 不同顆粒體積分?jǐn)?shù)下的復(fù)合材料應(yīng)力－應(yīng)變關(guān)系曲線Fig.7 Influence of volume fraction of particles on the stress－ strain behavior

采用應(yīng)力－應(yīng)變曲線在真應(yīng)變?yōu)?.002處的正割模量來(lái)近似炭黑填充橡膠材料的有效彈性模量＜E＞［4］，即

其中:ε=0.002，σ為應(yīng)力－應(yīng)變曲線上ε=0.002時(shí)的應(yīng)力值(MPa)，所以，有效彈性模量＜E＞的單位也為MPa.由此方法，可以得到不同顆粒含量下的炭黑填充橡膠材料的有效彈性模量＜E＞，如圖8所示.同時(shí)，圖8還給出了Bergstrom＆Boyce通過(guò)實(shí)驗(yàn)得到的7%、15%和25%炭黑N600填充氯丁橡膠的有效彈性模量的實(shí)驗(yàn)結(jié)果.將圓形粒子和方形粒子的有限元模型對(duì)不同含量炭黑填充的橡膠復(fù)合材料有效彈性模量的預(yù)測(cè)結(jié)果與實(shí)驗(yàn)結(jié)果進(jìn)行比較分析，發(fā)現(xiàn)它們所呈現(xiàn)的趨勢(shì)一致.即，炭黑顆粒的填充使得橡膠材料的彈性模量明顯增大，炭黑填充橡膠材料的有效彈性模量隨著炭黑顆粒所占體積分?jǐn)?shù)的增加而增大;表明本研究基于周期性邊界條件所建立的二維多顆粒夾雜RVE模型是合理的，能夠用于炭黑顆粒增強(qiáng)橡膠基復(fù)合材料有效性能的模擬分析.另外，從圖8中還可以看出，在相同的炭黑填料體積分?jǐn)?shù)下，方形粒子模型對(duì)炭黑填充橡膠材料的有效彈性模量的預(yù)測(cè)結(jié)果明顯高于圓形粒子模型.

圖8 不同炭黑顆粒含量下橡膠復(fù)合材料的有效彈性模量Fig.8 Influence of volume fraction of particles on the normalized Young’s modulus

4 結(jié)語(yǔ)

通過(guò)建立具有隨機(jī)分布形態(tài)的代表性體積單元的二維數(shù)值模型，對(duì)炭黑顆粒填充橡膠復(fù)合材料的力學(xué)行為進(jìn)行非線性有限元分析.研究表明，該方法通過(guò)周期性邊界條件的約束使得模型在整個(gè)加載過(guò)程符合周期性結(jié)構(gòu)的位移連續(xù)和應(yīng)力連續(xù)的真實(shí)情況，保證了變形場(chǎng)的協(xié)調(diào)性，可較好地模擬顆粒增強(qiáng)橡膠基復(fù)合材料的等效力學(xué)性能.結(jié)果顯示:

1)炭黑的填充明顯增大了橡膠材料的彈性模量;

2)炭黑填充橡膠材料的有效彈性模量隨著炭黑顆粒所占體積分?jǐn)?shù)的增加而增大;

3)在相同的炭黑填料體積分?jǐn)?shù)下，方形粒子模型對(duì)炭黑填充橡膠材料有效彈性模量的預(yù)測(cè)結(jié)果明顯高于圓形粒子模型;

4)本研究中圓形粒子和方形粒子的有限元模型對(duì)不同含量的炭黑填充橡膠材料有效彈性模量的預(yù)測(cè)結(jié)果與實(shí)驗(yàn)結(jié)果吻合較好，證實(shí)了本研究基于周期性邊界條件所建立的二維多顆粒夾雜RVE模型是合理的，為炭黑顆粒增強(qiáng)橡膠基復(fù)合材料其他性能的分析提供了依據(jù).

［1］王艷秋.橡膠材料基礎(chǔ)［M］.北京:化學(xué)工業(yè)出版社，2006.

［2］傅政.橡膠材料性能與設(shè)計(jì)應(yīng)用［M］.北京:化學(xué)工業(yè)出版社，2003.

［3］Govindjee S.An evaluation of strain amplification concepts via Monte Carlo simulations of an ideal composite［J］.Rubber Chemistry and Technology，1997，70:25－37.

［4］Bergstrom J S，Boyce M C.Mechanical behavior of particle filled elastomers［J］.Rubber Chemistry and Technology，1999，72:633－656.

［5］夏勇.輪胎膠料在較大變形范圍內(nèi)準(zhǔn)靜態(tài)力學(xué)性能的研究——測(cè)試、表征以及細(xì)觀數(shù)值本構(gòu)模型［D］.合肥:中國(guó)科學(xué)技術(shù)大學(xué)，2004.

［6］Gitman I M，Askes H，Sluys L J.Representative volume:existence and size determination［J］.Engineering Fracture Mechanics，2007，74:2 518 －2 534.

［7］Hassani B，Hinton E.A review of homogenization and topology optimization II－analytical and numerical solution of homogenization equations［J］.Computers and Structures，1998，69:719 － 738.

［8］Smit R J M，Brekelmans W A M，Meijer H E H.Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi－ level finite element modeling［J］.Computer Methods in Applied Mechanics and Engineering，1998，155:181 －192.

［9］Smit R J M，Brekelmans W A M，Meijer H E H.Prediction of the large－strain mechanical response of heterogeneous polymer systems:local and global deformation behavior of a representative volume element of voided polycarbonate［J］.Journal of the Mechanics and Physics of Solids，1999，47:201 －221.

［10］Segurado J，Llorca J.A numerical approximation to the elastic properties of sphere － reinforced composites［J］.Journal of the Mechanics and Physics of Solids，2002，50:2 107 －2 121.

［11］Boyce M C，Arruda E M.Constitutive models of rubber elasticity:a review［J］.Rubber Chemistry and Technology，2000，73(3):504－552.

［12］張少實(shí)，莊茁.復(fù)合材料與粘彈性力學(xué)［M］.北京:機(jī)械工業(yè)出版社，2005.

［13］Gent A N.Engineering with rubber:how to design rubber components［M］.Lübeck:Hanser Gardner Publications，2001.