基于三維離散元法的瀝青混合料斷裂過程模擬*

陳俊 黃曉明

(1.河海大學(xué) 土木與交通學(xué)院，江蘇南京210098;2.東南大學(xué) 交通學(xué)院，江蘇南京210096)

瀝青混凝土的斷裂性能一直以來都是道路工作者關(guān)注的熱點.在過去的幾十年間，人們采用了多種研究手段對瀝青混凝土的斷裂性能進行測試和分析，其中間接拉伸試驗和小梁彎曲斷裂試驗等室內(nèi)宏觀測試是最為普遍的研究方法.此類方法能容易獲得瀝青混凝土的斷裂強度和極限應(yīng)變，并能以此為基礎(chǔ)，進行基于斷裂性能的瀝青混凝土組成材料的優(yōu)選.但是，由于瀝青混凝土的力學(xué)性能在很大程度上取決于其內(nèi)部集料和空隙的體積含量、排列分布特征，而宏觀尺度的室內(nèi)試驗卻不能充分考慮混凝土的細(xì)觀結(jié)構(gòu)，因此此種方法對于深入理解瀝青混合料材料斷裂本質(zhì)特征、準(zhǔn)確把握材料的斷裂機理不能起到積極作用.

近年來，伴隨著圖像處理技術(shù)和計算機科學(xué)的快速發(fā)展，瀝青混凝土內(nèi)部空間細(xì)觀結(jié)構(gòu)的獲取引起人們越來越多的關(guān)注.在獲得瀝青混凝土結(jié)構(gòu)的基礎(chǔ)上，人們采用諸如有限元、離散元和邊界元之類的計算機分析軟件，試圖建立細(xì)觀結(jié)構(gòu)與宏觀性能的關(guān)系.如Sadd等［1-2］曾采用有限元方法從細(xì)部結(jié)構(gòu)入手預(yù)測瀝青混合料的宏觀模量.在此基礎(chǔ)上，Dai等［3］還運用有限元研究了瀝青混合料模量的各向異性問題.然而，隨著有限元方法在材料領(lǐng)域的廣泛應(yīng)用，其難以模擬粗集料和瀝青砂漿的界面形狀、不適合處理剛?cè)峤佑|問題和模擬材料大變形的缺陷逐漸被人們所認(rèn)識.

作為一種新發(fā)展起來的數(shù)值模擬方法，離散元法憑借其能夠處理應(yīng)力不連續(xù)、大變形問題等方面的優(yōu)勢，在道路材料領(lǐng)域的應(yīng)用越來越廣泛［4-7］.總體而言，離散元法在瀝青混合料性能模擬方面的研究大體可分為兩個階段.第一階段為2008年前，這期間Rothenburg等［8-12］對二維離散元方法應(yīng)用于瀝青混合料性能分析的可行性進行了研究，其主要研究思路是:采用對瀝青混合料試件的截面進行數(shù)碼照相;運用圖像處理技術(shù)，把數(shù)碼圖像轉(zhuǎn)化為數(shù)字試件，并導(dǎo)入到離散元軟件內(nèi)，生成瀝青混合料的二維數(shù)字試件;最后，進行瀝青混合料的力學(xué)模擬試驗.這期間研究工作的特點是:①以二維離散元分析為主，沒有開展三維分析;②瀝青材料假定為彈性，忽略其粘彈性特征.第二階段是2008年后，You等［13］開展了瀝青混合料三維結(jié)構(gòu)重構(gòu)的研究，研究方法是對瀝青混合料試件每隔一定高度h(≥8mm)進行切割，并對切割后的截面數(shù)碼照相，獲得n個二維數(shù)字試件，然后在離散元軟件內(nèi)把n個二維試件重構(gòu)為三維試件.重構(gòu)過程中，在高度為h的空間內(nèi)采用的是同一個二維試件.這樣做雖然能建立三維試件，但是試件內(nèi)部集料的形狀、空間布局等等都與實際試件差距較大，導(dǎo)致了集料與集料的嵌鎖作用不能得以體現(xiàn).另外，關(guān)于瀝青材料的粘彈性特征，Liu等［14］建立了Burgers模型中各元件微觀參數(shù)的確定方法，但并沒有應(yīng)用于三維結(jié)構(gòu)的分析.由此可見，離散元方法應(yīng)用于瀝青混合料性能模擬分析仍有諸多問題尚待解決，表現(xiàn)為:①如何在瀝青混合料三維結(jié)構(gòu)性能模擬中考慮瀝青材料的粘彈性特征;②二維模型的建立方法受照相技術(shù)、圖像處理過程的影響較大，而且試件的成型只是為了獲得其截面的數(shù)碼相片，造成了較多人力和物力的浪費;且混合料二維模型并不能考慮集料的嵌鎖作用和空隙的分布特征;③瀝青混合料真實三維空間結(jié)構(gòu)的獲取需要采用具有較大功率的CT掃描設(shè)備，而目前這些設(shè)備在我國道路材料研究領(lǐng)域還相當(dāng)缺乏，如何運用離散元技術(shù)直接建立能夠考慮集料顆粒不規(guī)則形狀、集料級配特征和空隙特征的混合料試件仍需研究.

為此，文中采用三維離散元顆粒流軟件(PFC3D)，通過自行編寫的顆粒生成程序，結(jié)合顆粒隨機投放算法，建立了可以考慮集料顆?？臻g不規(guī)則形狀、集料級配特征和空隙大小的混合料小梁試件的三維模型;在－10和15℃溫度條件下進行了考慮粘彈性特征的小梁試件中點彎曲試驗?zāi)M，并與實際試驗結(jié)果比較，驗證了離散元模型的正確性.

1 瀝青混合料的三維離散元模型

1.1 集料顆粒的空間不規(guī)則形狀

大量研究表明，集料尤其是粗集料顆粒的不規(guī)則形狀對瀝青混合料宏觀力學(xué)行為具有顯著的影響［15］，在建立瀝青混合料三維離散元模型時，應(yīng)當(dāng)充分考慮集料的不規(guī)則形狀.為此，參照Lu等［16］關(guān)于不規(guī)則顆粒生成的研究成果，直接在PFC3D內(nèi)以多個球形單元按照一定的算法相互重疊形成“clump”模擬顆粒.圖1為不同粒徑粗集料的典型顆粒形狀.

圖1 典型顆粒的不規(guī)則形狀Fig.1 Representative irregular shape of aggregate particles

對于上述生成的任一顆粒，其體積按照下式計算

式中:n、m分別為一個集料顆粒內(nèi)所包含的球形單元數(shù)和重疊個數(shù);Vi為集料顆粒內(nèi)第i個球形單元的體積;Voverlapj為第j個重疊部分的體積，

式中:R1，j、R2，j分別為形成第 j個重疊的兩個球形單元的半徑;h1，j、h2，j分別為兩個球形單元重疊部分球冠的高度.將式(1)和式(2)編制子程序嵌入到PFC3D內(nèi)，可以方便地計算出各個集料顆粒的體積，在集料密度設(shè)定后，即可獲得每個顆粒的質(zhì)量，這就為按照集料質(zhì)量級配建立混合料離散元模型奠定了基礎(chǔ).需要說明的是，在集料顆粒生成的過程中，顆粒的棱角豐富程度、集料的形狀(扁平率等)都可以通過顆粒生成算法內(nèi)的形狀參數(shù)加以改變.

1.2 集料混合料的三維離散元模型

考慮到瀝青混合料中2.36mm以下的集料與瀝青膠結(jié)料組成的瀝青砂漿材料的力學(xué)特征比較穩(wěn)定，兼顧到離散元計算效率的要求，文中的瀝青混合料離散元模型認(rèn)為是由2.36mm以上的粗集料顆粒和瀝青砂漿、空隙所組成.把上述集料顆粒按照一定的級配，投放于長、寬和高分別為300 mm、50 mm和50mm的瀝青混合料小梁試件的空間內(nèi).投放算法步驟如下:①按照集料粒徑由大至小逐級投放;②每投放完畢一個集料顆粒，判斷該顆粒與已投放顆粒是否重疊、該顆粒是否在設(shè)定的投放空間區(qū)域內(nèi);若滿足要求，則可以進行下一集料顆粒的投放，否則重復(fù)本步驟;③判斷該顆粒投放后，該檔集料質(zhì)量是否達到設(shè)定的質(zhì)量，若沒達到，則開始該檔下一集料的投放，若達到，則開始投放粒徑較小的下一檔集料.

采用上述投放算法，按表1中的級配進行粗集料顆粒的投放，形成如圖2所示的具有級配特征的粗集料混合物，其中13.2～16mm為120g(8個集料顆粒)、9.5 ～13.2mm 為360g(145 個顆粒)、4.75 ～9.5mm 為 480 g(810 個顆粒)、2.36 ～4.75 mm 為240g(4505個顆粒).值得注意的是，粗集料投放過程中，各檔集料的質(zhì)量和粗集料總體的質(zhì)量可以按照研究任務(wù)的不同任意設(shè)定，也就是說上述算法充分考慮了集料的級配特征和粗集料的體積含量.在集料混合物所在的長方體空間內(nèi)，按照單元矩形排列的方式(每個球形單元與周圍的6個相鄰單元相接觸)布滿半徑為1mm的球形單元，共93750個.判斷新布滿的單元與原有集料顆粒的位置關(guān)系，若重疊，則將新單元視為集料單元，若不重疊，則視為瀝青砂漿單元.刪除原有集料顆粒后，就建立了包含粗集料和瀝青砂漿的瀝青混合料梁式試件，如圖3所示.

表1 集料級配Table 1 Gradation of the aggregate blend

圖2 具有級配特征的集料顆?；旌衔颋ig.2 Aggregate blend with the given gradati

圖3 瀝青混合料梁式試件的離散元模型Fig.3 Discrete element model of beam-type specimem of asphalt mixture

為了考慮瀝青混合料內(nèi)部的空隙，在圖3所示三維模型的瀝青砂漿內(nèi)部，按照空隙的大小隨機刪除一定數(shù)量的球形單元，以模擬空隙在混合料內(nèi)部的存在.這樣處理的理由是，根據(jù)Masad等［17］的研究，瀝青混合料截面上的單個空隙面積在0.56～3mm2之間，與文中球形單元的面積相當(dāng).圖4為在如圖3所示的混合料三維結(jié)構(gòu)內(nèi)隨機刪除了2812個砂漿單元以模擬3%的空隙率的模型.

圖4 瀝青混合料內(nèi)部的空隙分布Fig.4 Illustration of air void in asphalt concrete

2 離散元微觀力學(xué)參數(shù)

2.1 宏觀與微觀力學(xué)參數(shù)

在上述離散元模型內(nèi)，涉及到瀝青砂漿內(nèi)部單元、集料內(nèi)部單元以及集料單元與砂漿單元的接觸與粘結(jié)問題，文中以PFC3D內(nèi)的線彈性剛度模型模擬集料內(nèi)部單元的接觸狀態(tài)，以Burgers模型表征瀝青砂漿間單元的粘彈性接觸.由于PFC3D內(nèi)材料性能的輸入是微觀參數(shù)，因此對于規(guī)則排列的單元可以由宏觀力學(xué)參數(shù)推導(dǎo)得到微觀參數(shù).

對于線彈性模型，采用下式確定微觀結(jié)構(gòu)內(nèi)單元線彈性接觸剛度來作為PFC3D內(nèi)參數(shù)的輸入:

式中:kn為單元的法向接觸剛度;R為單元的半徑;E為材料宏觀模量.

關(guān)于Burgers模型內(nèi)參數(shù)的確定，Liu等［14］建立了模型內(nèi)各元件微觀參數(shù)與宏觀參數(shù)的關(guān)系:

式中:L為瀝青混合料離散元模型內(nèi)相鄰單元的半徑之和;E1、η1分別為Maxwell內(nèi)彈簧的勁度和粘壺的黏度;E2、η2分別為Kelvin內(nèi)彈簧勁度和粘壺的黏度;Kmn、Cmn分別為Maxwell內(nèi)彈簧微觀勁度和粘壺的微觀黏度;Kkn、Ckn分別為Kelvin內(nèi)彈簧微觀勁度和粘壺的微觀黏度.E1、η1、E2和 η2可以根據(jù)下式計算得到:

式中:E*為復(fù)數(shù)彈性模量;ω為試驗頻率;θ為滯后角.

2.2 宏觀力學(xué)參數(shù)的測試

上述微觀參數(shù)的計算過程中，需要確定材料(尤其是瀝青砂漿)的Burgers模型的各元件參數(shù)，以及瀝青砂漿材料的抗拉強度.為此，在室內(nèi)首先成型瀝青砂漿的圓柱體試件，瀝青砂漿的級配按照表1所示的2.36mm粒徑以下的細(xì)集料級配確定，假定瀝青膜厚度為8 μm，計算得到的瀝青用量為14%.需要說明的是，瀝青砂漿試件成型時，由于其具有較大的瀝青用量，因而往往不需要外力擊實就能自密成型.在SPT材料試驗機上，進行－10和15℃下頻率為10Hz的動態(tài)壓縮試驗，并按照式(8)－(11)計算得到4個元件的參數(shù)，結(jié)果見表2.

表2 －10和15℃時的Burgers模型參數(shù)Table 2 Parameters in Burger’s model at －10 and 15℃

此外，按照上述級配和瀝青用量成型瀝青砂漿的馬歇爾試件，在UTM試驗機上進行材料的劈裂試驗，測試得到的－10和15℃抗拉強度，分別為5.0和 2.1 MPa.對于集料的模量按照 You［10］的研究成果，取55.5GPa.考慮到瀝青混合料彎曲斷裂時集料顆粒斷裂的情形很少，因此集料內(nèi)部單元之間的抗拉強度可以取大值，確保其不破壞即可.

3 瀝青混合料彎曲斷裂模擬與驗證

對圖3所示的瀝青混合料離散元模型，約束兩端底部的豎向自由度，對小梁的中點位置施加速率為50 mm/min的豎向荷載，進行三維結(jié)構(gòu)的彎曲斷裂模擬.圖5給出了三維結(jié)構(gòu)和其中某個二維截面的加載示意圖，以及二維截面斷裂后的形狀.

圖5 小梁彎曲斷裂的離散元模擬結(jié)果Fig.5 DEM simulation results of beams bending fracture

圖6給出了梁跨中底部不同拉應(yīng)變時，小梁中部60mm區(qū)域內(nèi)單元之間粘結(jié)失效情況，以清晰反映試件的斷裂過程.從圖中可以看出，隨著梁底拉應(yīng)變的增大，梁中部局部位置單元之間出現(xiàn)了粘結(jié)失效，由此可以清晰地看出梁內(nèi)裂紋的擴展情況.不僅如此，對上述的粘結(jié)失效，按照失效位置的不同，分為瀝青砂漿內(nèi)單元之間的粘結(jié)失效和砂漿與集料界面位置的粘結(jié)失效兩類，采用PFC分別采集了圖5(b)小梁內(nèi)兩類粘結(jié)失效的數(shù)量，如圖7所示.

圖6 不同拉應(yīng)變時小梁內(nèi)部的粘結(jié)失效情況Fig.6 Failed bonds in the beam at different levels of tensile strain

圖7 不同部位的粘結(jié)失效數(shù)量Fig.7 Number of failed bonds at different locations

從圖7中清晰地看出，粘結(jié)失效(裂紋)在瀝青砂漿和砂漿與集料界面上的數(shù)量隨著梁底拉應(yīng)變的變化情況.從上述分析和采集的數(shù)據(jù)可以看出，三維離散元模擬可以很好地模擬出梁式試件的斷裂過程，尤其是裂紋在試件內(nèi)的出現(xiàn)和擴展的數(shù)量、發(fā)生部位情形，這為從細(xì)觀結(jié)構(gòu)角度分析混合料斷裂特征提供了良好的手段.

在上述小梁彎曲斷裂的模擬過程中，采用PFC每隔10個計算時步采集梁底60mm范圍內(nèi)的拉應(yīng)變和豎向壓力，并按式(12)計算梁底拉應(yīng)力σt:

式中:L'為小梁試件的跨徑，mm;P為對小梁施加的最大壓力，N;w為小梁的平均寬度，mm;h為小梁平均高度，mm.同時，為了比較分析模擬結(jié)果與室內(nèi)實際試驗結(jié)果的差距，采用的室內(nèi)成型瀝青混合料小梁試件的級配和瀝青用量與斷裂模擬相同，試件的空隙率也按照三維離散元模型的3%加以設(shè)計.采用與小梁斷裂彎曲模擬相同的加載條件分別在－10和15℃溫度條件下進行室內(nèi)斷裂試驗，并采集了小梁彎曲斷裂過程中梁底拉應(yīng)力和拉應(yīng)變的關(guān)系曲線.圖8給出了不同溫度下小梁彎曲過程中跨中梁底拉應(yīng)力與拉應(yīng)變的實測結(jié)果與離散元模擬結(jié)果.

圖8 離散元模擬結(jié)果與實測結(jié)果的比較Fig.8 Comparisons between DEM simulation results and experimental ones

由圖8可以看出，在相同混合料級配、瀝青用量和空隙率時，三維離散元的模擬結(jié)果曲線與實測結(jié)果整體相差無幾，只是在曲線的局部位置處應(yīng)力和應(yīng)變值略有差異.其原因可能是:盡管室內(nèi)成型的試件與文中建立的三維離散元模型具有相同的瀝青用量、級配和空隙大小，但是它們內(nèi)部的細(xì)觀結(jié)構(gòu)尤其是粗集料和空隙的分布卻不盡相同，由此引起了斷裂曲線存在一定的差異.通過上述分析可見，文中所建立的三維離散元模型可以較好地模擬瀝青混合料彎曲斷裂尤其是裂紋出現(xiàn)和擴展的情形，其模擬結(jié)果基本上與實測結(jié)果保持一致，說明文中所建立的三維離散元模型能夠反映混合料的斷裂特征，三維離散元模擬可以作為研究瀝青混合料彎曲斷裂的輔助手段.

4 結(jié)語

文中運用離散元顆粒流程序PFC3D內(nèi)的“Fish”語言，編寫了集料顆?？臻g不規(guī)則形狀生成的程序，并結(jié)合顆粒隨機投放算法和計算機編程技術(shù)，建立了可以考慮粗集料形狀、級配特征和空隙大小的瀝青混合料小梁試件的三維離散元模型.三維離散元模型的建立方法為從微細(xì)觀結(jié)構(gòu)角度研究瀝青混合料的宏觀性能奠定了基礎(chǔ).

文中還根據(jù)材料宏觀力學(xué)性能與微觀力學(xué)參數(shù)的關(guān)系，結(jié)合瀝青砂漿的室內(nèi)性能試驗，獲得了瀝青砂漿粘彈性參數(shù)和各類微觀參數(shù)作為離散元力學(xué)試驗?zāi)M的輸入;并運用三維離散元方法，實現(xiàn)了小梁試件中點彎曲斷裂試驗的模擬，通過觀察模擬過程的現(xiàn)象，以及與室內(nèi)實際試驗結(jié)果的比較分析，驗證了離散元彎曲斷裂模擬的正確性.小梁彎曲斷裂模擬為研究瀝青混合料斷裂性能提供了輔助手段.

［1］Sadd M H，Dai Q L.A comparison of micro-mechanical modeling of asphalt materials using finite elements and doublet mechanics［J］.Mechanics of Materials，2005，37(6):641-662.

［2］Dai Q L，Sadd M H，You Z P.A micromechanical finite element model for linear and damage-coupled viscoelastic behaviour of asphalt mixture［J］.International Journal for Numerical and Analytical Methods in Geomechanics，2006，30(11):1135-1158.

［3］Dai Q L，You Z P.Prediction of creep stiffness of asphalt mixture with micromechanical finite-element and discreteelement models［J］.Journal of Engineering Mechanics，2007，133(2):163-169.

［4］Cundall P A，Strack O D L.A discrete numerical model for granular assemblies［J］.Géotechnique，1979，29(1):47-65.

［5］陳俊，黃曉明.基于離散元法的瀝青混合料虛擬疲勞試驗方法［J］.吉林大學(xué)學(xué)報:工學(xué)版，2010，40(2):435-440.Chen Jun，Huang Xiao-ming.Virtual fatigue test of asphalt mixture based on discrete element method［J］.Journal of Jilin University:Engineering and Technology Edition，2010，40(2):435-440.

［6］Chen Jun，Huang Xiao-ming.Virtual fracture test of asphalt mixture based on discrete element method［J］.Journal of Southeast University:English Edition，2009，25(4):518-522.

［7］陳俊，黃曉明.集料分布特征對混合料疲勞性能的影響分析［J］.建筑材料學(xué)報，2009，12(4):442-447.Chen Jun，Huang Xiao-ming.Research on influence of distribution characteristics of aggregate on fatigue performance of asphalt mixture ［J］.Journal of Building Materials，2009，12(4):442-447.

［8］Rothenburg L，Bathurst R J.Numerical simulation of idealized granular assemblies with plane elliptical particles［J］.Computers and Geotechnics，1991，11:315-329.

［9］Buttlar W G，You Z P.Discrete element modeling of asphalt concrete:a micro-fabric approach［J］.Transportation Research Record:Journal of the Transportation Research Board，2001，1757:111-118.

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

［11］Kim H，Buttlar W G.Discrete fracture modeling of asphalt concrete［J］.International Journal of Solids and Structures，2009，46(13):2593-2604.

［12］Kim H，Wagoner M P.Micromechanical fracture modeling of asphalt concrete using a single-edge notched beam test［J］.Materials and Structures，2009，42(5):677-689.

［13］You Z P，Adhikari S，Dai Q L.Three-dimensional discrete element models for asphalt mixtures［J］.Journal of Engineering Mechanics，2008，134(12):1053-1062.

［14］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.

［15］Pan T，Tutumluer E，Carpenter S H.Effect of coarse aggregate morphology on permanent deformation behavior of hot mix asphalt［J］.ASCE Journal of Transportation Engineering，2006，132(7):580-589.

［16］Lu M，Mcdowell G R.The importance of modelling ballast particle shape in the discrete element method ［J］.Granular Matter，2007，9:69-80.

［17］Masad E，Jandhyala V K，Dasgupta N，et al.Characterization of air void distribution in asphalt mixes using X-ray computed tomography［J］.Journal of Materials in Civil Engineering，2002，14(2):122-129.