顆粒棱度指標(biāo)的改進(jìn)及其對(duì)剪切特性的影響

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

(地質(zhì)災(zāi)害防治與地質(zhì)環(huán)境保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室(成都理工大學(xué))，四川 成都 610059)

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

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

1 單顆粒三維模型及三維形狀指標(biāo)

1.1 顆粒三維模型的建立

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

圖1 試驗(yàn)顆粒模型Fig.1 Particle model of the experiment

1.2 顆粒三維形狀指標(biāo)

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

(1)

式中：t——檢測(cè)步；

rp——顆粒在球坐標(biāo)上的極半徑；

rEE——等效橢球體在極坐標(biāo)的極半徑。

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

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

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

(2)

圖2 第i個(gè)四面體Fig.2 The tetrahedron i

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

(3)

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

(4)

其邊界條件Ω為：

(5)

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

(6)

(7)

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

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

(8)

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

(9)

式(9)中：

顆粒慣性張量B：

(10)

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

(11)

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

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

(12)

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

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

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

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

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

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

圖3 顆粒b等效橢球體對(duì)比圖Fig.3 Comparison diagram of the equivalent ellipsoid of the particle b(注：自左向右依次為三維圖形的俯視圖、左視圖、正視圖以及三維視圖，其中黑色線框?yàn)轭w粒b輪廓面，藍(lán)色線框?yàn)閺埩繖E球體，紅色線框?yàn)榍蛑C函數(shù)橢球體。)

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

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

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

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

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

表2 5種顆粒各形狀指標(biāo)

2 PFC3d直剪實(shí)驗(yàn)數(shù)值模擬

2.1 參數(shù)標(biāo)定

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(13)

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

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

φ=7.621lnAI+63.653

(14)

3 結(jié)論

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

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

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

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