李慶，楊曉翔

(福州大學機械工程及自動化學院，福建福州 350116)

0 引言

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

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

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

1 復合材料細觀力學研究基礎

1.1 代表性體積單元

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

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

1.2 周期性邊界條件

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

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

圖2所示為二維代表性體積單元的幾何表述，其外框為一邊長為L的正方形，頂點A0為坐標原點，頂點A1和A2分別位于坐標軸x和y上，則對邊邊界上點的自由度應滿足以下約束關系:

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

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

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

由邊界虛位移引起的虛應變能變化為:

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

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

因為RVE對邊邊界形狀相同，l1=l2，l3=l4，則式(7)可寫為:

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

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

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

對邊邊界上點的應力應滿足以下約束關系:

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

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

2 有限元建模

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

圖3 RVE模型Fig.3 RVE model

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

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

橡膠基體本構關系采用近乎不可壓縮的Mooney－Rivlin模型［11－12］

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

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

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

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

3 結果與討論

3.1 變形場和應力場分析

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

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

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

圖5 RVE的von－Mises應力等值線圖Fig.5 The distribution of von － Mises stress in RVE

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

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

從圖6中可以看出應力與應變呈非線性關系，符合橡膠材料的非線性彈性的特性.所以，本文研究所建立的模型可以用來模擬炭黑填充橡膠復合材料的力學行為.

3.2 炭黑顆粒含量和形狀對復合材料有效彈性模量的影響

圖6 RVE的平均應力－應變曲線Fig.6 The stress－ strain behavior of RVE

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

圖7 不同顆粒體積分數(shù)下的復合材料應力－應變關系曲線Fig.7 Influence of volume fraction of particles on the stress－ strain behavior

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

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

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

4 結語

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

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

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

3)在相同的炭黑填料體積分數(shù)下，方形粒子模型對炭黑填充橡膠材料有效彈性模量的預測結果明顯高于圓形粒子模型;

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

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

［2］傅政.橡膠材料性能與設計應用［M］.北京:化學工業(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］夏勇.輪胎膠料在較大變形范圍內準靜態(tài)力學性能的研究——測試、表征以及細觀數(shù)值本構模型［D］.合肥:中國科學技術大學，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］張少實，莊茁.復合材料與粘彈性力學［M］.北京:機械工業(yè)出版社，2005.

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