袁 斌，霍宇翔，巨能攀，宋浩燃

(地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護國家重點實驗室(成都理工大學)，四川 成都 610059)

顆粒材料作為建筑材料廣泛應(yīng)用于交通土建工程中，而顆粒形狀作為顆粒材料一個及其重要的屬性參數(shù)，其對顆粒宏細觀力學性質(zhì)的影響一直是顆粒材料研究的一個熱點問題。為了更好地了解顆粒形狀對顆粒材料力學性質(zhì)的影響，學者們進行了大量的室內(nèi)試驗和數(shù)值模擬試驗研究。常曉林等[1]利用三維變形體離散元單元法研究了顆粒的長短徑比在宏觀和細觀層面對堆石體力學性能的影響；王蘊嘉等[2]利用PFC3d研究了堆石料顆粒形狀對堆積密度及強度的影響；張成功等[3]利用顆粒料軟件研究了顆粒形狀對顆粒材料圓柱塌落的影響；王鵬程等[4-5]利用PFC2d對碎石集料雙軸壓縮試驗進行了數(shù)值模擬，研究了顆粒形狀對不良級配碎石集料剪切特性的影響；鄒德高等[6]利用三維激光掃描技術(shù)量化了堆石料三維形狀特征，并研究了其對顆粒破碎的影響。已有研究成果表明，不同顆粒形狀有著不同的歷史成因，顆粒之間的嵌擠和摩擦受顆粒形狀的影響，使整體強度發(fā)生變化[7-10]。對于顆粒形態(tài)的量化，國內(nèi)外學者很早就開展了一系列研究。Barrett[11]提出顆粒形態(tài)特征包括顆粒形狀、棱角性以及顆粒表面紋理等三個獨立的特性，這三個特性分別代表顆粒在大、中、小三個尺寸上的變化規(guī)律。棱角度作為一種中尺度形態(tài)表征參數(shù)，最早由Wadell[12]于1932年提出的圓度(Roundness)進行評價并被廣泛應(yīng)用于工程實踐中。Masad[13]等基于等效橢球提出了棱角度AI(Angular index)以評價顆粒表面棱角起伏程度，Zhou[14]等基于球諧函數(shù)對該等效橢球體的計算方法進行了改進。

但已有研究對于顆粒形狀的描述大多還處于二維平面層面，且形狀較為規(guī)則，缺乏精細模型的三維模擬。本文通過PFC3d軟件建立了直接剪切試驗條件，利用軟件內(nèi)置的clump算法生成了不規(guī)則碎石顆粒精細三維模型，同時基于文獻[15]對于顆粒形狀的研究，改進了棱度指數(shù)的概念，并將其作為形狀評價指標，定量研究了顆粒形狀對顆粒材料剪切特性的影響。

1 單顆粒三維模型及三維形狀指標

1.1 顆粒三維模型的建立

為了獲取原始顆粒的三維形狀信息，利用CT掃描獲得顆粒逐層二維CT影像[16-18]，再將CT影像導(dǎo)入MATLAB中，利用對影像進行二值化、邊界檢測處理后生成的像素信息矩陣提取顆粒輪廓點的相對坐標，從而得到顆粒表面的點云集。將點云導(dǎo)入三維建模軟件建立顆粒三維模型，將其輪廓進行Delaunay 三角化后，利用PFC3d中的clump 算法對其進行填充得到單顆顆粒(圖1)。

圖1 試驗顆粒模型Fig.1 Particle model of the experiment

1.2 顆粒三維形狀指標

棱角度AI由Masad[13]等基于顆粒等效橢球提出，其計算公式如下：

(1)

式中：t——檢測步；

rp——顆粒在球坐標上的極半徑；

rEE——等效橢球體在極坐標的極半徑。

該等效橢球體為體積等效橢球體，即橢球體體積等于顆粒體積，并且該橢球體長短徑比與原顆粒一致。

本研究考慮到在假定顆粒不會破損的情況下，顆粒體的運動主要為平移和旋轉(zhuǎn)，其旋轉(zhuǎn)受外力以及其本身的轉(zhuǎn)動慣量張量控制，因此顆粒等效橢球體的重心體積主慣性力矩應(yīng)與顆粒一致，即在其長軸、中軸、短軸方向上的轉(zhuǎn)動慣量應(yīng)與顆粒體一致。對此改進了顆粒等效橢球體的計算公式。

以模型形心為原點建立空間直角坐標系，假設(shè)模型輪廓由n個三角形面組成，連接原點和每個輪廓點，則可將模型分為n個四面體。

(2)

圖2 第i個四面體Fig.2 The tetrahedron i

由于計算轉(zhuǎn)動張量需要用三重積分：

(3)

但四面體的邊界在三維直角坐標系中表達太復(fù)雜，所以導(dǎo)入?yún)?shù)方程對其進行簡化：

(4)

其邊界條件Ω為：

(5)

將式(4)代入式(3)進行三重積分換元：

(6)

(7)

式中：detAi——矩陣Ai的特征值。

將式(7)代入式(6)得到換元后的三重積分公式：

(8)

假設(shè)顆粒為均質(zhì)材料且密度為ρ，則第i個四面體對原點的慣性張量Bi計算式如下：

(9)

式(9)中：

顆粒慣性張量B：

(10)

假設(shè)過原點任意軸的方向向量(單位向量)為(α，β，γ)，則該軸的轉(zhuǎn)動慣量I表達式為：

(11)

利用式(11)可求出過原點任意旋轉(zhuǎn)軸的轉(zhuǎn)動慣量，最大轉(zhuǎn)動慣量旋轉(zhuǎn)軸方向?qū)?yīng)等效橢球體短軸方向，最小轉(zhuǎn)動慣量旋轉(zhuǎn)軸方向?qū)?yīng)等效橢球體長軸方向，中軸方向同時垂直于長軸和短軸，可用長軸方向向量與短軸方向向量叉乘得出，并將中軸方向向量代入式(11)求出相應(yīng)轉(zhuǎn)動慣量。

等效橢球體長半軸長a、中半軸長b、短半軸長c可由以下方程組求解：

(12)

式中：Imin，Imax——顆粒繞過形心任意轉(zhuǎn)軸的最小、最大轉(zhuǎn)動慣量；

Imid——顆粒繞等效橢球體中軸的轉(zhuǎn)動慣量；

M——橢球體的質(zhì)量。

計算等效橢球體各軸長及方向后，將直角坐標系轉(zhuǎn)化為球坐標系，采用式(1)計算顆粒棱度AI。

為了方便表述，將利用轉(zhuǎn)動張量矩陣得到的等效橢球體簡稱為張量橢球體，將利用體積等效得到的等效橢球體簡稱為體積橢球體。

用2種方法分別求出了顆粒b的等效橢球體，并作出三維圖形進行對比(圖3)。分別計算了2個等效橢球體與實體顆粒的體積之比(V/VT)、表面積之比(S/ST)、長短徑比之比(EI/EIT)、重心體積主慣性力矩之比(I1/I1T、I2/I2T、I3/I3T)，結(jié)果如表1所示。

圖3 顆粒b等效橢球體對比圖Fig.3 Comparison diagram of the equivalent ellipsoid of the particle b(注：自左向右依次為三維圖形的俯視圖、左視圖、正視圖以及三維視圖，其中黑色線框為顆粒b輪廓面，藍色線框為張量橢球體，紅色線框為球諧函數(shù)橢球體。)

表1 兩種等效橢球體參數(shù)比較

張量橢球體因與顆粒重心體積主慣性力矩等效，所以其I1/I1T、I2/I2T、I3/I3T值都為1，顆粒與其張量橢球體體積的差距在4%以內(nèi)，表面積和長短徑比的差距在20%以內(nèi)。體積等效橢球體因與顆粒體積等效且長短徑比也相同，所以其V/VT、EI/EIT值都為1，顆粒與其體積等效橢球體表面積的差距也在20%以內(nèi)，略小于張量橢球體，3個方向上慣性力矩的差距都超過了10%，最大差距超過30%。

綜上所述，2種等效橢球體與實體顆粒體積的差距、表面積的差距相差不大，雖然體積等效橢球體的長短徑比與實體顆粒一致，但張量橢球體更加符合顆粒的力學行為。

圖 4為5種不同形狀的顆粒及其棱度指數(shù)，可以看出顆粒AI值越小，顆粒就越接近于一個橢球體，其輪廓也越規(guī)則。表2為5種顆粒各形狀指標，其中SI為球度即顆粒等體積球體的表面積與顆粒實際表面積的比值，EI為伸長率即顆粒短軸長與長軸長的比值，e為扁平率即顆粒短軸長與中軸長的比值。

圖4 不同形狀顆粒及其棱度指數(shù)Fig.4 Different particle shapes and its AI

表2 5種顆粒各形狀指標

2 PFC3d直剪實驗數(shù)值模擬

2.1 參數(shù)標定

為了標定數(shù)值模擬的參數(shù)，進行了粒徑為 1～2 cm的砂巖碎石在圍壓為0.1，0.2，0.3 MPa 條件下的室內(nèi)中型直接剪切試驗，砂巖的抗壓強度為 30 MPa。試樣的高度為16 cm、寬度為14 cm、長度為16 cm。試樣密度為2 650 kg/m3，孔隙比約為0.65。

同時，利用PFC3d[19]建立了碎石集料的中型直接剪切試驗條件，用10個墻體模擬直剪盒。為了盡可能模擬顆粒間的鑲嵌咬合作用，在室內(nèi)試驗試樣中選取10個形狀較為典型的顆粒進行掃描建模，并導(dǎo)入PFC 中利用clump算法對其進行填充后作為生成顆粒的剛性簇模板，再由這10個模板等量地隨機生成粒徑范圍為1～2 cm形狀不規(guī)則的碎石，最終的數(shù)值模擬模型如圖5所示。

圖5 混合顆粒離散元模型Fig.5 Discrete element model of the mixed particle

模擬過程中顆粒之間的接觸模型采用線性接觸剛度模型，需要定義接觸的切法向剛度及摩擦系數(shù)。表3為本文所選取碎石顆粒的主要計算參數(shù)。為了與室內(nèi)試驗進行對比，所有數(shù)值模型的終止剪切應(yīng)變均為12%，剪切速率也與室內(nèi)試驗保持一致。

圖6為室內(nèi)試驗以及數(shù)值模擬計算得到的應(yīng)力-應(yīng)變對比曲線，與室內(nèi)試驗結(jié)果對比，可以看出數(shù)值模型的曲線要比室內(nèi)試驗曲線略高。在室內(nèi)試驗過程

表3 顆粒計算參數(shù)

圖6 室內(nèi)試驗和數(shù)值計算得到的應(yīng)力-應(yīng)變對比曲線Fig.6 Comparison of stress-strain of the laboratory tests and numerical calculations

中，粒間摩擦會產(chǎn)生細小粉塵降低粒間摩擦系數(shù)以及發(fā)生顆粒破碎的現(xiàn)象，而這一過程在數(shù)值模型中并沒有得到體現(xiàn)。因此，數(shù)值模型得到的碎石剪切強度要比室內(nèi)試驗得到的結(jié)果略高?？傮w上來說，數(shù)值模型試驗計算結(jié)果基本可以描述曲線的主要形態(tài)，驗證了模型及計算參數(shù)選取的合理性。

2.2 單一顆粒形狀直剪試驗數(shù)值模擬

在標定完參數(shù)后，對圖 4所示的 5種顆粒逐一進行單一顆粒形狀直剪試驗數(shù)值模擬，采用的接觸模型及相關(guān)參數(shù)、直剪盒尺寸、顆粒粒徑、剪切速率與參數(shù)標定時一致。

圖7為數(shù)值模擬試驗得到的在0.3 MPa圍壓下不同顆粒碎石集料的應(yīng)力-應(yīng)變曲線，可以看出，不同形狀顆粒粒組的應(yīng)力-應(yīng)變曲線在形態(tài)上有很大的不同。從整體上看，曲線為非線性，初始剪切階段線性較為明顯，而隨著剪切應(yīng)變的增加，顆粒集料剪應(yīng)力在應(yīng)變達到10% 后均沒有太大的變化。

圖7 不同顆粒形狀下應(yīng)力-應(yīng)變數(shù)值模擬試驗曲線Fig.7 Stress-strain of the numerical simulation test under different particle shapes

根據(jù)摩爾-庫倫準則整理得到不同形狀顆粒的抗剪強度線(圖 8)，抗剪強度隨圍壓的增加顯著增加。顆粒內(nèi)摩擦角即抗剪強度線的斜率隨顆粒形狀的變化也明顯不同。

圖8 不同顆粒形狀下圍壓與抗剪強度關(guān)系Fig.8 Relationship between the confining pressure and shear strength under different particle shapes

圖9為 0.3 MPa圍壓下5個顆粒各個形狀參數(shù)與抗剪強度的關(guān)系，從圖9中可看出，顆粒棱度AI較其他3個形狀指標與抗剪強度有更好的對應(yīng)關(guān)系。

圖9 不同形狀指標與抗剪強度關(guān)系Fig.9 Relationship between the different shape indexes and shear strength

在驗證了模型的正確性后，對另外11種塊石顆粒進行建模，分別計算了顆粒的棱度指標。再分別導(dǎo)入PFC3d生成clump模板進行直剪試驗后，用試驗結(jié)果對比顆粒棱度進行分析。

在圍壓相同的情況下，抗剪強度隨棱度的增大而增大，在棱度大于0.3后，抗剪強度增加的幅度逐漸減小(圖 10)，可用對數(shù)函數(shù)進行擬合：

(13)

圖10 棱度指數(shù)對抗剪強度的影響Fig.10 Effect of AI on the shear strength

計算出每個塊體的內(nèi)摩擦角后，與顆粒棱度指標進行對比，內(nèi)摩擦角隨棱度的增加而增加，增加幅度隨棱度的增加而減小(圖 11)，其關(guān)系也可用對數(shù)函數(shù)進行擬合，擬合關(guān)系式為：

φ=7.621lnAI+63.653

(14)

3 結(jié)論

(1)提出了用塊體轉(zhuǎn)動張量確定塊體的等效橢球體的方法，將其與球諧函數(shù)確定方法進行比較，論證了其合理性，并在此基礎(chǔ)上改進了顆粒三維棱度指標的計算公式。

(2)由碎石顆粒直剪試驗應(yīng)力-應(yīng)變曲線可得出，隨著應(yīng)變的不斷增加，剪應(yīng)力不斷增加，但增加的速率逐漸趨近于零，有明顯的極值。曲線的形態(tài)和極值跟顆粒形狀有很大的關(guān)系。

(3)級配不良碎石集料的抗剪強度主要受其顆粒棱度SI的影響，隨著顆粒棱度的增加，其抗剪強度和內(nèi)摩擦角都明顯增加，其關(guān)系可以用對數(shù)函數(shù)進行擬合，內(nèi)摩擦角與顆粒棱度自然對數(shù)擬合直線的斜率為7.621、截距為63.653。在棱度較低的時候，抗剪強度和內(nèi)摩擦角增加較快，當棱度較大時，抗剪強度指標增加幅度較小。

(4)本研究中的室內(nèi)試驗和數(shù)值模擬試驗樣品顆粒級配不良，對于級配良好顆粒形狀的評價還應(yīng)結(jié)合其他指標深入研究。