瀝青混合料三維離散元虛擬單軸蠕變?cè)囼?yàn)*

張德育 黃曉明 高英

(東南大學(xué)交通學(xué)院，江蘇南京210096)

隨著公路交通量的增加和重載車輛的增多，車轍已成為瀝青路面的主要病害之一，嚴(yán)重影響路面行車安全和舒適性，公路瀝青路面的車轍已經(jīng)成為人們普遍關(guān)注的問(wèn)題.對(duì)于瀝青混合料的永久變形性能，各國(guó)道路工作者提出了多種研究方法，其中以現(xiàn)象試驗(yàn)法最為普遍，但這種方法不僅需要消耗大量時(shí)間和成本，重復(fù)性和再現(xiàn)性也較差［1-2］.單軸蠕變?cè)囼?yàn)作為常用的瀝青混合料高溫永久變形性能試驗(yàn)亦存在這種缺點(diǎn).

瀝青混合料由集料、瀝青膠結(jié)料及空隙構(gòu)成，集料在瀝青混合料中呈現(xiàn)不連續(xù)的顆粒特性，空隙結(jié)構(gòu)亦導(dǎo)致瀝青混合料呈現(xiàn)不連續(xù)特性，因此，瀝青混合料是一種多相復(fù)合材料，其材料性質(zhì)較為復(fù)雜，實(shí)際的應(yīng)力、應(yīng)變場(chǎng)并不連續(xù)，集料、瀝青膠結(jié)料及空隙在混合料中扮演著不同的力學(xué)角色，呈現(xiàn)出不均勻和離散的力學(xué)特性.長(zhǎng)期以來(lái)，道路工程領(lǐng)域一直沿用連續(xù)均質(zhì)的力學(xué)方法分析瀝青混合料的永久變形性能并指導(dǎo)材料設(shè)計(jì)，這與瀝青混合料的不連續(xù)特征不符，難以真實(shí)地反映瀝青混合料的材料力學(xué)特性.

近年來(lái)，離散元方法逐漸被引入到瀝青混合料性能分析中.離散元方法可以反映瀝青混合料的不連續(xù)特性，從細(xì)觀角度研究瀝青混合料的永久變形性能，建立的虛擬試驗(yàn)的重復(fù)性和再現(xiàn)性好，經(jīng)濟(jì)方便，能夠克服傳統(tǒng)試驗(yàn)的不足.但在以往的研究中，幾乎都是采用二維(2D)離散元模型對(duì)瀝青混合料進(jìn)行分析［3-8］，而二維結(jié)構(gòu)與三維(3D)結(jié)構(gòu)相差較大，二維離散元模型具有一定的局限性.因此，文中采用離散元軟件PFC3D生成了包含集料、瀝青砂漿及空隙的瀝青混合料三維離散元虛擬試件，建立了瀝青混合料三維離散元虛擬單軸蠕變?cè)囼?yàn)，并與二維虛擬試驗(yàn)及室內(nèi)試驗(yàn)結(jié)果進(jìn)行了比較，驗(yàn)證了三維離散元虛擬試驗(yàn)方法的正確性.

1 三維離散元虛擬試件的建立

為了兼顧離散元模型的計(jì)算效率和瀝青砂漿粘彈性參數(shù)的獲取，將粒徑大于1.18 mm的集料視為粗集料，粒徑小于1.18 mm的細(xì)集料和瀝青膠結(jié)料視為均質(zhì)的瀝青砂漿，因此，文中建立的瀝青混合料離散元虛擬試件由集料、瀝青砂漿和空隙構(gòu)成.首先根據(jù)瀝青混合料的油石比、空隙率及集料的級(jí)配(見表1)等計(jì)算出各檔集料的體積，然后將各檔集料逐次投放到由墻構(gòu)成的直徑為100 mm，高為150 mm的瀝青混合料虛擬試件尺寸空間中，并通過(guò)循環(huán)保證生成的集料球單元不重疊.生成的具有級(jí)配特征的集料如圖1所示.

表1 AC-16瀝青混合料級(jí)配Table 1 Gradation of AC-16 asphalt mixture

圖1 具有級(jí)配特征的集料Fig.1 Graded aggregates

待具有級(jí)配特征的集料生成后，在瀝青混合料試件尺寸空間內(nèi)生成148200個(gè)規(guī)則排列(每個(gè)單元與四周的6個(gè)單元相接觸)且半徑為1mm的離散單元(如圖2所示).遍歷新生成的離散單元，逐一判斷新單元與原集料球單元的位置是否重疊，如重疊則將新單元視為集料單元，否則視為瀝青砂漿單元，并將每個(gè)集料范圍內(nèi)的離散單元定義為剛性“聚粒”，以節(jié)省計(jì)算時(shí)間.判斷結(jié)束后刪除原有集料單元，生成如圖3所示的三維瀝青混合料虛擬試件.

圖2 規(guī)則排列的離散單元Fig.2 Rugularly-packed discrete elements

圖3 三維瀝青混合料虛擬試件Fig.3 3D virtual specimen of asphalt mixture

由于空隙對(duì)瀝青混合料的永久變形性能影響較大，為了在瀝青混合料離散元虛擬試件中生成空隙，在瀝青砂漿單元中隨機(jī)刪除5928個(gè)單元作為4%的空隙.三維離散元虛擬試件中的空隙分布如圖4所示.由于室內(nèi)試件及現(xiàn)場(chǎng)取芯試件內(nèi)部空隙分布的復(fù)雜性，文中將虛擬試件內(nèi)部的空隙分布視為隨機(jī)分布.

圖4 瀝青混合料空隙分布Fig.4 Distribution of air voids in asphalt mixture

2 微觀接觸模型及參數(shù)

2.1 微觀接觸模型

圖5 瀝青混合料微觀組成成分的接觸行為Fig.5 Interactions among microscale components in asphalt mixture

在瀝青混合料離散元模型中，瀝青混合料內(nèi)部微觀組成成分之間有4種接觸，分別為集料內(nèi)部單元、相鄰集料之間、瀝青砂漿內(nèi)部單元、瀝青砂漿和集料之間的接觸，如圖5所示.根據(jù)PFC3D中不同接觸模型的特點(diǎn)，采用接觸剛度模型描述集料內(nèi)部單元之間的接觸行為，采用接觸剛度模型和滑移模型描述相鄰集料之間的接觸行為，采用伯格斯接觸模型和接觸粘結(jié)模型描述瀝青砂漿內(nèi)部單元之間以及瀝青砂漿單元和集料之間的接觸行為.

2.2 模型參數(shù)

PFC3D中微觀接觸模型的參數(shù)很難直接從室內(nèi)試驗(yàn)得到，需要建立微觀模型參數(shù)與宏觀試驗(yàn)參數(shù)之間的關(guān)系，通過(guò)室內(nèi)試驗(yàn)得到宏觀試驗(yàn)參數(shù)，進(jìn)而根據(jù)已建立的轉(zhuǎn)化關(guān)系得到微觀模型參數(shù).

考慮到瀝青砂漿的粘彈性，采用接觸粘結(jié)模型和微觀伯格斯接觸模型(如圖6所示)來(lái)描述瀝青砂漿內(nèi)部單元之間及瀝青砂漿與集料之間的接觸行為.圖6 中 E1、η1、E2、η2為宏觀伯格斯模型參數(shù)，Kmn、Cmn、Kkn、Ckn為微觀伯格斯接觸模型法向參數(shù)，Kms、Cms、Kks、Cks為微觀伯格斯接觸模型切向參數(shù)，m1、m2表示顆粒單元，fs為摩擦系數(shù).

圖6 伯格斯模型Fig.6 Burger’s model

由于文中主要研究瀝青混合料的蠕變特性，不考慮材料的破壞，接觸粘結(jié)強(qiáng)度取較大的值以避免單元粘結(jié)的破壞.瀝青砂漿的粘彈性參數(shù)可以由DSR動(dòng)態(tài)剪切試驗(yàn)獲取，文中采用的瀝青砂漿粘彈性參數(shù)分別為 E1=4.64 MPa，η1=83.34 MPa·s，E2=0.82MPa，η2=10.00 MPa·s.根據(jù)已有的研究成果［9］，微觀伯格斯接觸模型參數(shù)(見圖6(b)和圖6(c))與宏觀伯格斯模型參數(shù)(見圖6(a))之間的轉(zhuǎn)化關(guān)系為

式中，L為相鄰單元的球心距，υ為瀝青砂漿泊松比.

將集料視為彈性體，采用接觸剛度模型來(lái)描述集料內(nèi)部單元的接觸行為，采用接觸剛度模型和滑動(dòng)模型來(lái)描述集料之間的接觸行為.由于離散單元之間的排列方式為規(guī)則的矩形排列，根據(jù)已有文獻(xiàn)的研究成果［10］，集料的接觸剛度與宏觀彈性模量的轉(zhuǎn)化關(guān)系為

式中:kn、ks分別為接觸剛度模型的法向剛度和切向剛度;E 為集料的宏觀彈性模量，文中取 55.5GPa［4，6，11］;R為集料單元的半徑;υ'為集料泊松比.

3 時(shí)溫等效原理

粘彈性離散元模型能夠更好地模擬瀝青混合料的粘彈性行為，而粘彈性模型的計(jì)算相當(dāng)耗費(fèi)時(shí)間［9］，以 CPU 2.8 GHz、RAM 2.0 GB 計(jì)算機(jī)進(jìn)行真實(shí)時(shí)間下的虛擬試驗(yàn)將會(huì)花費(fèi)數(shù)年的時(shí)間.因此，如果計(jì)算中采用真實(shí)加載時(shí)間，將會(huì)對(duì)虛擬試驗(yàn)的進(jìn)行造成很大的困難，甚至無(wú)法進(jìn)行.瀝青混合料作為粘彈性材料，其變形依賴于時(shí)間和溫度，為了減少計(jì)算時(shí)間，文中采用基于時(shí)溫等效原理的方法對(duì)虛擬試驗(yàn)的計(jì)算效率進(jìn)行優(yōu)化.

時(shí)溫等效原理是將材料對(duì)于外界荷載的響應(yīng)進(jìn)行時(shí)間和溫度的等效轉(zhuǎn)化，因此可以提高粘彈性材料的加載溫度進(jìn)而減少加載時(shí)間，再將高溫下的結(jié)果通過(guò)時(shí)溫等效原理向低溫條件下進(jìn)行轉(zhuǎn)化，轉(zhuǎn)化關(guān)系為

式中:ε(T，t)為真實(shí)溫度T和加載時(shí)間t下的蠕變應(yīng)變;ε(Tr，tr)為參考溫度Tr和縮減時(shí)間tr下的蠕變應(yīng)變;αT為移位因子.

在瀝青混合料的研究中，時(shí)溫等效原理常用來(lái)建立主曲線，進(jìn)而可以將一定時(shí)間、溫度范圍內(nèi)的試驗(yàn)測(cè)定結(jié)果拓延到更廣泛的時(shí)間溫度空間中去，大大減少材料研究中的試驗(yàn)工作，并在相當(dāng)程度上降低對(duì)試驗(yàn)裝置的技術(shù)要求.同樣的，可以利用時(shí)溫等效原理減少粘彈性離散元模型的計(jì)算時(shí)間.瀝青混合料粘彈性伯格斯模型中真實(shí)溫度和時(shí)間下的蠕變應(yīng)變可以表達(dá)為

式中，σ0為軸向應(yīng)力.

當(dāng)溫度升高到參考溫度時(shí)，粘彈性參數(shù)也會(huì)隨之改變.假定 E1、E2、η1r、η2r為參考溫度下的伯格斯模型參數(shù)，則在縮減時(shí)間下的蠕變應(yīng)變?yōu)?/p>

瀝青混合料在真實(shí)溫度和參考溫度下的伯格斯模型參數(shù)是不同的，但不同溫度下的蠕變應(yīng)變?cè)诓煌募虞d時(shí)間下可以進(jìn)行等效.通過(guò)求解式(11)－(14)，可以得到參考溫度下的伯格斯模型參數(shù):

在以往的文獻(xiàn)中，離散元粘彈性模擬分析都是采用真實(shí)時(shí)間下的伯格斯模型參數(shù)，因此計(jì)算時(shí)間會(huì)很長(zhǎng).如果采用參考溫度和縮減時(shí)間下的伯格斯模型參數(shù)，計(jì)算時(shí)間就取決于縮減時(shí)間.當(dāng)移位因子取較大值時(shí)，縮減時(shí)間較真實(shí)時(shí)間大幅度縮短，因此模型的計(jì)算時(shí)間將會(huì)大大減少.經(jīng)過(guò)試算，當(dāng)移位因子不大于104時(shí)，計(jì)算結(jié)果較為穩(wěn)定(如圖7所示)，這與You等［12］的研究結(jié)論較為一致.因此，文中取104作為移位因子的取值，那么虛擬試驗(yàn)的加載時(shí)間僅為真實(shí)時(shí)間的1/104，計(jì)算時(shí)間從數(shù)年減少到了數(shù)天.

圖7 不同移位因子下的瀝青混合料軸向應(yīng)變Fig.7 Axial strain of asphalt mixture with different shifting factors

4 虛擬單軸蠕變?cè)囼?yàn)

4.1 三維虛擬試驗(yàn)結(jié)果

在上述建立的虛擬試件基礎(chǔ)上，對(duì)虛擬試件施加0.7 MPa的靜載.在PFC3D中，通常是通過(guò)控制墻的移動(dòng)速度來(lái)對(duì)試件進(jìn)行加載，因此本虛擬試驗(yàn)通過(guò)在試件的上下表面各設(shè)定一面墻，固定試件下表面的墻，通過(guò)PFC3D中的fish語(yǔ)言編寫伺服控制程序，不斷調(diào)整試件上表面墻的移動(dòng)速度使軸向應(yīng)力達(dá)到恒定的值.由于文中采用了基于時(shí)溫等效的計(jì)算優(yōu)化方法且移位因子取為104，當(dāng)室內(nèi)試驗(yàn)加載時(shí)間為100 s時(shí)，虛擬試驗(yàn)僅需加載0.01 s，計(jì)算時(shí)間大大減少.

由于瀝青路面車轍通常發(fā)生在高溫條件下，尤其是接近瀝青軟化點(diǎn)的時(shí)候車轍量會(huì)迅速增加［13］.因此，文中的試驗(yàn)溫度選擇具有代表性的60℃.瀝青混合料試件的直徑為100 mm，高為150 mm，最佳油石比為4.8%，空隙率為4%，混合料級(jí)配如表1所示.分別進(jìn)行相同條件下的三維離散元虛擬試驗(yàn)和室內(nèi)試驗(yàn)，并進(jìn)行比較分析，結(jié)果如圖8所示.

圖8 三維虛擬試驗(yàn)與室內(nèi)試驗(yàn)結(jié)果比較Fig.8 Comparing of results of 3D virtual test and laboratory measurement

從圖8可以看出，三維離散元虛擬試驗(yàn)結(jié)果與真實(shí)試驗(yàn)的結(jié)果較為吻合，但不完全一致，其原因可能是虛擬試件中集料的形狀與真實(shí)試件中集料的形狀有所差別，虛擬試件內(nèi)部空隙率的分布與真實(shí)試件也有所不同，而且室內(nèi)試驗(yàn)的加載條件也沒(méi)有虛擬試驗(yàn)理想.從總體上看，三維離散元虛擬試驗(yàn)?zāi)軌蜉^為準(zhǔn)確地預(yù)估瀝青混合料的永久變形性能，可以作為瀝青混合料永久變形性能分析的輔助手段.

4.2 三維和二維虛擬試驗(yàn)結(jié)果比較

為了比較三維和二維離散元模型預(yù)估瀝青混合料永久變形性能的準(zhǔn)確性，文中還進(jìn)行了相同試驗(yàn)條件下的二維虛擬單軸蠕變?cè)囼?yàn).二維離散元模型直接由三維離散元模型相互垂直的兩個(gè)切片獲得(如圖9所示).三維和二維離散元模型對(duì)比結(jié)果如圖10(a)所示，不同荷載下的三維和二維離散元模型對(duì)比結(jié)果如圖10(b)所示.

圖9 瀝青混合料二維離散元模型的獲取Fig.9 Illustration of obtaining 2D DE model of asphalt mixture

圖10 三維和二維離散元模型計(jì)算結(jié)果比較Fig.10 Comparing of results of 3D and 2D DE models

從圖10(a)可以看出:加載初期軸向應(yīng)變的三維預(yù)測(cè)值大于二維預(yù)測(cè)值;隨著加載的進(jìn)行，三維預(yù)測(cè)值則小于二維預(yù)測(cè)值，且差異隨之增大.從圖10(b)可以看出，軸向應(yīng)變的二維預(yù)測(cè)值較三維預(yù)測(cè)值及室內(nèi)試驗(yàn)值要大，且隨著荷載的增大差異也隨之增大.分析其原因可能是二維模型低估了集料的嵌鎖作用.三維離散元模型在預(yù)測(cè)瀝青混合料高溫永久變形性能方面較二維離散元模型更為理想，可以準(zhǔn)確地預(yù)估瀝青混合料的高溫永久變形性能.

5 結(jié)語(yǔ)

文中建立了包括瀝青混合料三維離散元虛擬試件的生成、瀝青混合料微觀組成成分微觀接觸模型的選取、微觀接觸模型參數(shù)的獲取以及加載方式的實(shí)現(xiàn)等在內(nèi)的完整的瀝青混合料三維虛擬單軸蠕變?cè)囼?yàn)方法，提出了基于時(shí)溫等效原理的離散元粘彈性模型計(jì)算優(yōu)化方法，將瀝青混合料粘彈性離散元模型的計(jì)算時(shí)間從數(shù)年減少到了數(shù)天，從而解決了三維離散元粘彈性模型的計(jì)算時(shí)間問(wèn)題.

二維虛擬試驗(yàn)結(jié)果大于三維虛擬試驗(yàn)及室內(nèi)試驗(yàn)結(jié)果，且隨著荷載的增大差異也隨之增大，說(shuō)明三維離散元模型在預(yù)測(cè)瀝青混合料高溫永久變形性能方面較二維離散元模型更為理想，所建立的瀝青混合料三維離散元虛擬單軸蠕變?cè)囼?yàn)?zāi)軌驕?zhǔn)確地預(yù)估瀝青混合料高溫永久變形性能，可以作為瀝青混合料永久變形性能分析的輔助手段.

［1］譚積青.粗集料形態(tài)與瀝青混合料級(jí)配組成虛擬力學(xué)試驗(yàn)的基礎(chǔ)方法研究［D］.廣州:華南理工大學(xué)土木與交通學(xué)院，2006.

［2］胡霞光.瀝青混合料微觀力學(xué)分析綜述［J］.長(zhǎng)安大學(xué)學(xué)報(bào):自然科學(xué)版，2005，25(2):6-9.Hu Xia-guang.Review on asphalt mixture micromechanics analysis［J］.Journal of Chang＇an University:Natural Science Edition，2005，25(2):6-9.

［3］Chang G K，Meegoda J N.Micromechanical simulation of hot mix asphalt［J］.Journal of Engineering Mechanics，1997，123(5):495-503.

［4］Buttlar W G，You Z P.Discrete element modeling of asphalt concrete:microfabric approach ［J］.Transportation Research Record，2001，1757:111-118.

［5］王端宜，張肖寧，王紹懷.用虛擬試驗(yàn)方法評(píng)價(jià)瀝青混合料的級(jí)配類型［J］.華南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2003，31(2):48-51.Wang Duan-yi，Zhang Xiao-ning，Wang Shao-huai.Evaluation on grading type of asphalt mixture with virtual test method［J］.Journal of South China University of Technology:Natural Science Edition，2003，31(2):48-51.

［6］You Z P，Buttlar W G.Discrete element modeling to predict the modulus of asphalt concrete mixtures［J］.Journal of Materials in Civil Engineering，2004，16(2):140-146.

［7］Abbas A，Masad E，Paapagiannakis T，et al.Micromechanical modeling of the viscoelastic behavior of asphalt mixtures using the discrete element method［J］.International Journal of Geomechanics，2007，7(2):131-139.

［8］Kim H，Wagoner W P，Buttlar W G.Simulation of fracture behavior in asphalt concrete using a heterogeneous cohesive zone discrete element model［J］.Journal of Materials in Civil Engineering，2008，20(8):552-563.

［9］Liu Y，Dai Q L，You Z P.Viscoelastic model for discrete element simulation of asphalt mixtures［J］.Journal of Engineering Mechanics，2009，135(4):324-333.

［10］Thornton C.The conditions for failure of a face-centered cubic array of uniform rigid spheres［J］.Geotechnique，1979，29(4):441-459.

［11］You Z P.Development of a micromechanical modeling approach to predict asphalt mixture stiffness using the discrete element method［D］.Urbana-Champaign:Department of Civil and Environmental Engineering，University of Illinois at Urbana-Champaign，2003.

［12］You Z P，Liu Y，Dai Q L.Three-dimensional microstructural-based discrete element viscoelastic modeling of creep compliance tests for asphalt mixtures［J］.Journal of Materials in Civil Engineering，2011，23(1):79-87.

［13］黃曉明，范躍武，趙永利，等.高速公路瀝青路面高溫車轍的調(diào)查與分析［J］.公路交通科技，2007，24(5):16-20.Huang Xiao-ming，F(xiàn)an Yue-wu，Zhao Yong-li，et al.Investigation and test of expressway asphalt pavement hightemperature performance［J］.Journal of Highway and Transportation Research and Development，2007，24(5):16-20.