顆粒摩擦對(duì)散粒堆積體拱效應(yīng)的影響*

戴北冰， 鄧林杰， 陳智剛

1.中山大學(xué)土木工程學(xué)院，廣東 珠海 519082

2.南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室（珠海），廣東 珠海 519082

3.重慶建工第一市政工程有限責(zé)任公司，重慶 400020

散粒材料在自然界和人類生產(chǎn)生活中普遍存在（Terzaghi，1936；Karl，1943）。散體顆粒在重力或其他外力因素作用下形成的堆狀顆粒集合體稱為散粒堆積體（Nedderman，1992；賈海莉等，2003）。散粒堆積體的自然休止角以及內(nèi)部的拱效應(yīng)是其重要的物理力學(xué)特征，在農(nóng)業(yè)、工業(yè)、制藥等領(lǐng)域有著十分重要的應(yīng)用背景和研究價(jià)值（Dai et al.，2017； Dai，2018）。

散粒堆積體內(nèi)應(yīng)力的分布一直是學(xué)者們感興趣的研究對(duì)象（戴北冰等，2022；Dai et al.，2022）。學(xué)者們?cè)谏扯阎杏^察到了一種有趣且違反直覺的現(xiàn)象，即沙堆錐頂正下方的底部垂直應(yīng)力會(huì)出現(xiàn)局部最小值，即堆積體底部形成“M”形應(yīng)力分布（Wittmer et al.，1996；Watson，1996；Cates et al.，1999），這種現(xiàn)象稱作應(yīng)力凹陷現(xiàn)象。應(yīng)力凹陷與拱效應(yīng)緊密相關(guān)，堆積體內(nèi)部的應(yīng)力拱使應(yīng)力傳遞產(chǎn)生了屏蔽效應(yīng)，導(dǎo)致應(yīng)力分布局部化和不均勻化，從而引起沙堆頂部正下方底部應(yīng)力的局部下降（Fang et al.，2022）。因此，研究應(yīng)力凹陷現(xiàn)象對(duì)理解沙堆中的拱效應(yīng)至關(guān)重要。眾多學(xué)者已圍繞顆粒形狀（Zuriguel et al.，2007；Zuriguel et al.，2008；Zhao et al.，2019）、顆粒摩擦（Geng et al.，2001；Silbert et al.，2002；Goldenberg et al.，2005；Horabik et al.，2017）等影響因素，對(duì)應(yīng)力凹陷現(xiàn)象開展了較為深入的研究，發(fā)現(xiàn)：隨著顆粒摩擦的增大，沙堆底部壓力局部下降越明顯，應(yīng)力凹陷現(xiàn)象增強(qiáng)。Fang et al.（2022）發(fā)現(xiàn)，不同顆粒粒徑形成的沙堆內(nèi)摩擦角不一致，內(nèi)摩擦角大的沙堆底部應(yīng)力凹陷越明顯，這也間接佐證了顆粒粒徑會(huì)影響沙堆的應(yīng)力凹陷。Goldenberg et al.（2005）和Zhou et al.（2009）通過數(shù)值模擬發(fā)現(xiàn)：顆粒圓度越低，沙堆底部壓力局部下降越明顯。因此，顆?；拘再|(zhì)對(duì)散粒堆積體底部的應(yīng)力凹陷現(xiàn)象起著舉足輕重的作用。

底部應(yīng)力凹陷和拱效應(yīng)的產(chǎn)生本質(zhì)上是由力鏈的發(fā)展和傳遞實(shí)現(xiàn)的（Meng et al.，2018；Zaidi，2020），而顆粒摩擦是其重要的影響因素之一。許多學(xué)者研究了顆粒摩擦對(duì)力鏈傳遞的影響規(guī)律，發(fā)現(xiàn)散粒堆積體內(nèi)強(qiáng)弱力鏈的相對(duì)位置和數(shù)量會(huì)隨著顆粒摩擦的改變而改變（孫其誠等，2008；張煒等，2022），堆積體內(nèi)強(qiáng)力鏈傳遞了大部分荷載（Ma et al.，2016；戴 北 冰 等，2019）。張 煒 等（2022）在研究鐵粉末壓制過程中力鏈演化規(guī)律時(shí)，發(fā)現(xiàn)隨著顆粒間摩擦系數(shù)增大，整體力鏈數(shù)目變少，力鏈方向系數(shù)變大。Zhu et al.（2013）發(fā)現(xiàn)，對(duì)于較大的顆粒摩擦，強(qiáng)力鏈傾向于在更陡的方位上分布，形成更穩(wěn)定的拱形以承擔(dān)更多的沙堆重量。Bouchaud et al.（2002）通過光彈實(shí)驗(yàn)分析了顆粒集合體在外力作用下的排列和相互作用，以及力鏈的傳遞路徑。Majmudar et al.（2005）測量了光彈顆粒集合在純剪切和各向同性壓縮作用下法向力和切向力的傳遞路徑。由此可見，顆粒摩擦對(duì)散粒體堆積體中的力鏈傳遞以及拱效應(yīng)的發(fā)揮機(jī)制有著顯著影響。

目前，學(xué)者們雖然在關(guān)于顆粒摩擦對(duì)土拱效應(yīng)影響方面做了許多工作，但尚未深刻闡明顆粒摩擦對(duì)散粒堆積體內(nèi)部細(xì)觀力學(xué)響應(yīng)的影響規(guī)律以及不同顆粒摩擦影響下拱效應(yīng)產(chǎn)生的細(xì)觀力學(xué)機(jī)制。本研究旨在利用DEM 模擬，研究顆粒摩擦對(duì)散粒堆積體宏細(xì)觀力學(xué)行為的影響，通過將堆積體內(nèi)部大于平均接觸力的接觸力沿錐角方向進(jìn)行投影、統(tǒng)計(jì)顆粒接觸方位的概率密度分布來分析散粒堆積體內(nèi)部拱效應(yīng)的優(yōu)勢(shì)發(fā)揮方位，基于統(tǒng)計(jì)分析研究堆積體內(nèi)部強(qiáng)弱力鏈占比，對(duì)不同顆粒摩擦下散粒堆積體內(nèi)力鏈網(wǎng)絡(luò)進(jìn)行可視化分析，建立宏觀力學(xué)特征與細(xì)觀力學(xué)指標(biāo)的聯(lián)系，以揭示散粒堆積體內(nèi)拱效應(yīng)形成的細(xì)觀力學(xué)機(jī)制。

1 顆粒設(shè)計(jì)以及堆積體建模

本文使用DEM 程序PFC3D 來模擬沙堆的形成，數(shù)值模擬所用的材料是一種理想的彈性顆粒材料，如圖1 所示。顆粒形狀系數(shù)Ra（即子顆粒的直徑與顆粒的長度之比）為0.6，顆粒粒徑（與該顆粒體積相等的球體直徑）在0.01～0.015 m 之間分布，相關(guān)模型參數(shù)如表1所示（戴北冰等，2019）。數(shù)值模擬采用了8 種不同的顆粒摩擦系數(shù)(μ= 0.01，0.05，0.10，0.30，0.50，0.70，1.00、2.00)，以研究顆粒摩擦對(duì)散粒堆積體拱效應(yīng)的影響。所有模擬中，重力加速度均采用10 m/s2，底部墻體的摩擦系數(shù)為0.50。數(shù)值模擬中顆粒粒徑保持一致，不會(huì)影響所發(fā)現(xiàn)的顆粒摩擦影響規(guī)律。

表1 離散元數(shù)值模擬參數(shù)Table 1 Discrete element numerical simulation parameters

圖1 顆粒形狀Fig.1 Particle shape

本數(shù)值模擬通過點(diǎn)源釋放法來形成散粒堆積體，顆粒釋放工具如圖2所示。釋放器上方豁口是邊長為0.25 m的正方形，釋放器下方豁口是邊長為0.075 m 的正方形，釋放器高度為2.5 m，組成釋放器的所有墻體的摩擦系數(shù)均為0。在本數(shù)值模擬中，散粒堆積體內(nèi)一共有約2×104個(gè)顆粒，由釋放器分10次下落形成，每次下落都保證釋放器底部離堆積體最高處（第一次為水平面）為一固定高度0.5 m。

圖2 漏斗型釋放器Fig.2 Funnel type releaser

2 分析與討論

2.1 堆積體休止角計(jì)算及分析

圖3給出了不同顆粒摩擦條件下所形成的堆積體數(shù)值模型?？梢钥闯?，μ為0.01 和0.05 時(shí)，難以形成具有明顯坡角的堆積體；μ為0.10～0.50時(shí)，可以觀察到形成的堆積體有明顯的坡角，且坡角隨著μ的增大而增大；μ超過0.50后，坡角基本不變。圖4給出了自然休止角α與顆粒摩擦系數(shù)μ的關(guān)系。從圖3 和圖4 可以發(fā)現(xiàn)，在顆粒間摩擦系數(shù)μ小于0.50 時(shí)，自然休止角α隨著顆粒摩擦系數(shù)的增加而增加，在顆粒間摩擦系數(shù)μ> 0.50后，自然休止角α趨于一個(gè)穩(wěn)定值。

圖3 不同摩擦系數(shù)形成的散粒堆積體Fig.3 Granular heaps on different inter-particle friction

圖4 不同顆粒摩擦系數(shù)下堆積體的自然休止角Fig.4 α of granular heaps with different inter-particle friction

2.2 顆粒堆積體底部接觸力分析

圖5給出了不同摩擦條件下堆積體底部接觸力的云圖，X、Y坐標(biāo)的單位為m，色卡的數(shù)值單位為N?？梢园l(fā)現(xiàn)，μ為0.01 時(shí)，堆積體底部接觸力的峰值出現(xiàn)在中心位置；μ為0.05～2.00時(shí)，堆積體底部接觸力的峰值不在中心位置，而是出現(xiàn)在偏離中心的位置。進(jìn)一步觀察可以發(fā)現(xiàn)，隨著顆粒摩擦系數(shù)增大，堆積體底部中心區(qū)域的顏色由紅色（μ= 0.01）變 為 黃 色（μ= 0.05），然 后 是 綠 色（μ≥0.10），最后是藍(lán)色（μ≥0.70）。紅色代表應(yīng)力水平相對(duì)較高，而藍(lán)色代表應(yīng)力水平相對(duì)較低。因此，從中心區(qū)域顏色的變化趨勢(shì)可以判斷中心區(qū)域應(yīng)力水平相比較于整個(gè)底部區(qū)域應(yīng)力峰值水平的差距。隨著摩擦系數(shù)的增大差距增大，并在μ> 0.5 后，這種差距幾乎趨于穩(wěn)定，意即底部應(yīng)力凹陷現(xiàn)象隨著摩擦系數(shù)增加而變得顯著，并逐步趨于穩(wěn)定。

圖5 堆積體底部接觸力云圖Fig.5 Cloud chart of contact force at the bottom of granular heaps

圖6為堆積體底部法向接觸力和切向接觸力沿半徑方向的分布，x/R是參考位置到圓心距離與堆積體底部圓半徑的比值。可以發(fā)現(xiàn)，顆粒摩擦系數(shù)μ= 0.01時(shí)，底部法向和切向接觸力分布均未出現(xiàn)應(yīng)力凹陷現(xiàn)象。當(dāng)μ= 0.05時(shí)，底部法向和切向接觸力分布則出現(xiàn)了微弱的應(yīng)力凹陷現(xiàn)象。隨著μ繼續(xù)增大，堆積體底部法向和切向接觸力的峰值明顯地增大，并在μ> 0.50后變化不大，而堆積體底部中心處法向和切向接觸力的峰值則表現(xiàn)出先增大（μ< 0.5）后減小（μ> 0.5）的趨勢(shì)。

圖6 底部法向和切向接觸力沿半徑方向的分布Fig.6 Distribution of normal contact force and tangential contact force at the bottom along the radius

圖7-8給出了堆積體底部峰值接觸力、中心接觸力以及中心接觸力相對(duì)于峰值的減少程度隨顆粒摩擦系數(shù)的變化。從圖7可以發(fā)現(xiàn)，隨著顆粒摩擦系數(shù)的增大，堆積體底部接觸力峰值先增大，后逐漸趨于一個(gè)穩(wěn)定值；堆積體底部中心處的接觸力則先增大后減小。圖8顯示中心接觸力的減小程度隨著摩擦系數(shù)的增大而增大，亦有逐步趨于穩(wěn)定的趨勢(shì)，這進(jìn)一步佐證了底部接觸力云圖的觀察發(fā)現(xiàn)，即底部應(yīng)力凹陷現(xiàn)象隨著顆粒摩擦系數(shù)的增大而增強(qiáng)，并在顆粒摩擦較大的情況下趨于穩(wěn)定。

圖7 峰值接觸力與中心接觸力隨顆粒摩擦系數(shù)的變化Fig.7 The change of the peak contact force Fp and the center contact force Fo with inter-particle friction

圖8 中心接觸力相對(duì)于峰值的減少程度Fig.8 The decreasing degree of the center contact force relative to the peak value

圖9 給出了底部接觸力峰值位置xp/R與μ的關(guān)系。可以發(fā)現(xiàn)，隨著μ增大，峰值位置從底部中心逐漸往外遷移，并在μ> 0.50后，峰值位置維持在xp/R= 0.21左右。

圖9 底部接觸力峰值位置與顆粒摩擦系數(shù)的關(guān)系Fig.9 The relationship between the position of the peak contact force at the bottom and the inter-particle friction coefficient

2.3 顆粒間接觸力投影

散粒堆積體的形成與穩(wěn)定取決于其內(nèi)部顆粒間的相互作用，力鏈的傳遞對(duì)拱效應(yīng)的形成起著關(guān)鍵作用。為進(jìn)一步分析顆粒間接觸力的傳遞對(duì)拱效應(yīng)形成的影響以及拱效應(yīng)的優(yōu)勢(shì)發(fā)揮方位，將堆積體內(nèi)部大于平均接觸力的接觸力在如圖10所示的a?、i?、m?（其中m→為沿錐角方向，i→為與m→垂直的環(huán)向方向，a→為與i→、m→垂直的方向）三個(gè)方向進(jìn)行投影。

圖10 投影方式示意圖Fig.10 Schematic diagram of projection mode

圖11-13 為法向接觸力、切向接觸力、總接觸力沿a?、i?、m→三個(gè)方向的投影值，縱坐標(biāo)F為接觸力的投影值大小，橫坐標(biāo)為投影角度θ。

圖11 法向接觸力沿a→、i→、m→三個(gè)方向的投影值Fig.11 Projection value of normal contact force along a→、i→、m→

圖12 切向接觸力沿a→、i→、m→三個(gè)方向的投影值Fig.12 Projection value of tangential contact force along a→、i→、m→

從圖11-13 可以發(fā)現(xiàn)，沙堆內(nèi)部顆粒法向接觸力和切向接觸力在各個(gè)方向上的投影值大小隨著顆粒間摩擦系數(shù)的增大而增大，且在i→方向（環(huán)向）的投影幾乎為零，這也間接證明了應(yīng)力在環(huán)向是對(duì)稱的。法向接觸力在a→方向的投影值在θ= 15°左右時(shí)為0，而在θ= 105°左右達(dá)到最大。顯然，投影最大值和最小值的方位剛好差90°。由于m→與a→垂直，法向接觸力在m→方向的投影規(guī)律與a→方向的相反，投影值在θ= 105°左右為零點(diǎn)，而在θ= 15°左右時(shí)達(dá)到最大，即沿著m→與a→投影峰值（或者零點(diǎn)）的方位剛好相差90°。切向接觸力在a→方向的投影值在θ= 25°左右時(shí)為0，在θ=115°左右時(shí)達(dá)到最大，而在m→方向的投影值則在θ= 115°左右時(shí)為0，在θ= 25°左右時(shí)達(dá)到最大。如圖13 所示，由于切向接觸力較法向接觸力小很多，因此總接觸力的投影規(guī)律主要受法向接觸力控制，與圖11中法向接觸力投影規(guī)律幾乎相同。這里需注意的是，對(duì)于μ= 0.01 和0.05 的情況，法向接觸力和切向接觸力投影值均很小，這主要是因?yàn)轭w粒在這兩種摩擦狀態(tài)下很難形成有效的堆積。

圖13 總接觸力沿a→、i→、m→三個(gè)方向的投影值Fig.13 Projection value of total contact force along a→、i→、m→

圖14 給出了沿m→向（即錐面方向）投影值峰值所在的方位角度θ與摩擦系數(shù)的關(guān)系?？梢钥闯?，無論是法向接觸力、切向接觸力，還是總接觸力，投影峰值方位角度θ隨著μ的增大而增大，直至達(dá)到一個(gè)穩(wěn)定值。如果不考慮摩擦系數(shù)較低的情況（μ≤0.1），法向接觸力和總接觸力的投影峰值方位接近θ= 15°，而切向接觸力的投影峰值方位接近θ= 25°。

圖14 法向接觸力、切向接觸力和總接觸力沿m→向投影峰值所在的方位角Fig.14 The orientation angle when conical surface force,tangential contact force and total contact force reaches the peak value along m→

2.4 顆粒接觸方位向量

圖15為顆粒接觸方位的概率密度f分布。由于散粒堆積體基本呈錐形，可以近似認(rèn)為顆粒接觸方位在環(huán)向是對(duì)稱分布，因此這里只考慮顆粒接觸以豎直軸作為參考基準(zhǔn)的方位分布規(guī)律。γ為顆粒接觸方位相對(duì)于豎直軸的方位角，取值范圍為0°～90°。從圖15可以看出，隨著γ的增大，顆粒接觸方位的概率密度f基本呈先增加后降低的趨勢(shì)，但概率密度峰值fmax的方位γ因摩擦系數(shù)不同而不同。μ= 0.01 時(shí)，fmax出現(xiàn)在γ= 35～40°；μ= 0.05時(shí)，fmax出現(xiàn)在γ= 30～35°；μ= 0.10 時(shí)，fmax出現(xiàn)在γ= 20～30°；μ= 0.30～2.00 時(shí)，峰值出現(xiàn)在γ=15～25°的時(shí)候?？梢钥闯?，峰值方位角隨著摩擦系數(shù)增大而降低，但當(dāng)μ> 0.30時(shí)，峰值出現(xiàn)的方位（γ= 15～25°）不隨顆粒摩擦系數(shù)的增加而顯著改變，且該方位范圍與法向接觸力投影峰值的方位（θ≈15°）以 及 切 向 接 觸 力 投 影 峰 值 的 方 位（θ≈25°）非常接近（見圖14）。需指出的是，這里的θ方位亦以豎直軸為參考。這些觀察充分證明了在顆粒摩擦系數(shù)較高（μ ≥0.3）的情況下，散粒堆積體內(nèi)部拱效應(yīng)的優(yōu)勢(shì)發(fā)揮方位出現(xiàn)在偏離豎直軸15°～25°的方位。

圖15 顆粒接觸方位概率密度分布Fig.15 Probability density distribution of particle contact orientation

2.5 堆積體中力鏈統(tǒng)計(jì)

圖16 統(tǒng)計(jì)了堆積體中弱力鏈的出現(xiàn)概率。弱力鏈?zhǔn)侵附佑|力小于平均接觸力的力鏈。Pn是指法向接觸力弱力鏈的概率，Pt是指切向接觸力弱力鏈的概率。由圖16 可知，Pn隨著顆粒摩擦系數(shù)的增加而微弱增加，并在μ大于0.10后，維持在0.68 左右。Pt在摩擦系數(shù)較低時(shí)接近1.0，而在μ大于0.30 后，隨著μ的增加而降低?？梢园l(fā)現(xiàn)，不論是法向接觸力還是切向接觸力，堆積體內(nèi)弱力鏈的數(shù)量總是多于強(qiáng)力鏈的數(shù)量；切向接觸力中弱力鏈的占比要高于法向接觸力中弱力鏈的占比；隨著摩擦系數(shù)的增大，部分切向接觸力的弱力鏈會(huì)逐漸轉(zhuǎn)化為強(qiáng)力鏈。

圖16 弱力鏈統(tǒng)計(jì)Fig.16 The statistics of weak force chains

圖17 給 出 了μ= 0.05，0.10，0.30，0.50 時(shí) 堆 積體的力鏈網(wǎng)絡(luò)圖。圖中力鏈粗細(xì)（或顏色）表征接觸力的大小?？梢园l(fā)現(xiàn)，平均接觸力隨著μ增加而增大。當(dāng)μ= 0.05，0.10 時(shí)，散粒堆積體中的強(qiáng)力鏈數(shù)量偏少，且呈水平攤開的分布，難以觀察到成“拱”現(xiàn)象。從圖17 中可見，當(dāng)μ ≥0.3 時(shí)，隨顆粒摩擦系數(shù)的增大，強(qiáng)力鏈在水平和高度方向上逐漸擴(kuò)展，強(qiáng)力鏈也逐漸與堆積體中心軸呈一定傾斜角度，成“拱”效應(yīng)也越來顯著。

圖17 散粒堆積體力鏈網(wǎng)絡(luò)圖Fig.17 Force chains network diagram of granular heaps

3 結(jié) 論

1）顆粒摩擦?xí)绊懚逊e體休止角；摩擦系數(shù)越大，休止角越大；當(dāng)顆粒摩擦系數(shù)大于0.50 后，堆積體的休止角趨于一個(gè)穩(wěn)定值。

2）隨著顆粒摩擦系數(shù)增大，堆積體底部中心接觸力呈先增大后減小的趨勢(shì)；底部接觸力峰值先增大后趨于一個(gè)穩(wěn)定值，其位置從底部中心xp/R= 0 逐漸往外遷移，最終維持在xp/R= 0.21 左右；底部中心處接觸力值相對(duì)于峰值的減小程度ρ從0%增加到72.7%，堆積體底部應(yīng)力凹陷現(xiàn)象增強(qiáng)。

3）在μ≥0.30 時(shí)，法向接觸力和總接觸力在m→向的投影峰值方位接近θ= 15°，切向接觸力的投影峰值方位接近θ= 25°，且與接觸方位向量的概率密度峰值方位接近，證明散粒堆積體內(nèi)部拱效應(yīng)的優(yōu)勢(shì)發(fā)揮方位出現(xiàn)在偏離豎直軸15°～25°的方位上。

4）在摩擦系數(shù)較低時(shí)（μ≤0.10），強(qiáng)力鏈的數(shù)量偏少且呈水平攤開分布，難以形成明顯的拱效應(yīng)；當(dāng)摩擦系數(shù)較高時(shí)（μ≥0.30），強(qiáng)力鏈在豎直和水平方向擴(kuò)展，并形成傾斜傳遞，促進(jìn)拱效應(yīng)的產(chǎn)生。