顆粒摩擦對散粒堆積體拱效應的影響*

2024-01-10 10:50:48戴北冰鄧林杰陳智剛

中山大學學報(自然科學版)(中英文) 2023年6期

戴北冰，鄧林杰，陳智剛

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

2.南方海洋科學與工程廣東省實驗室（珠海），廣東珠海 519082

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

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

散粒堆積體內應力的分布一直是學者們感興趣的研究對象（戴北冰等，2022；Dai et al.，2022）。學者們在沙堆中觀察到了一種有趣且違反直覺的現象，即沙堆錐頂正下方的底部垂直應力會出現局部最小值，即堆積體底部形成“M”形應力分布（Wittmer et al.，1996；Watson，1996；Cates et al.，1999），這種現象稱作應力凹陷現象。應力凹陷與拱效應緊密相關，堆積體內部的應力拱使應力傳遞產生了屏蔽效應，導致應力分布局部化和不均勻化，從而引起沙堆頂部正下方底部應力的局部下降（Fang et al.，2022）。因此，研究應力凹陷現象對理解沙堆中的拱效應至關重要。眾多學者已圍繞顆粒形狀（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）等影響因素，對應力凹陷現象開展了較為深入的研究，發現：隨著顆粒摩擦的增大，沙堆底部壓力局部下降越明顯，應力凹陷現象增強。Fang et al.（2022）發現，不同顆粒粒徑形成的沙堆內摩擦角不一致，內摩擦角大的沙堆底部應力凹陷越明顯，這也間接佐證了顆粒粒徑會影響沙堆的應力凹陷。Goldenberg et al.（2005）和Zhou et al.（2009）通過數值模擬發現：顆粒圓度越低，沙堆底部壓力局部下降越明顯。因此，顆粒基本性質對散粒堆積體底部的應力凹陷現象起著舉足輕重的作用。

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

目前，學者們雖然在關于顆粒摩擦對土拱效應影響方面做了許多工作，但尚未深刻闡明顆粒摩擦對散粒堆積體內部細觀力學響應的影響規律以及不同顆粒摩擦影響下拱效應產生的細觀力學機制。本研究旨在利用DEM 模擬，研究顆粒摩擦對散粒堆積體宏細觀力學行為的影響，通過將堆積體內部大于平均接觸力的接觸力沿錐角方向進行投影、統計顆粒接觸方位的概率密度分布來分析散粒堆積體內部拱效應的優勢發揮方位，基于統計分析研究堆積體內部強弱力鏈占比，對不同顆粒摩擦下散粒堆積體內力鏈網絡進行可視化分析，建立宏觀力學特征與細觀力學指標的聯系，以揭示散粒堆積體內拱效應形成的細觀力學機制。

1 顆粒設計以及堆積體建模

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

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

圖1 顆粒形狀Fig.1 Particle shape

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

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

2 分析與討論

2.1 堆積體休止角計算及分析

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 顆粒間接觸力投影

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

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

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

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

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

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

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

圖14 給出了沿m→向（即錐面方向）投影值峰值所在的方位角度θ與摩擦系數的關系。可以看出，無論是法向接觸力、切向接觸力，還是總接觸力，投影峰值方位角度θ隨著μ的增大而增大，直至達到一個穩定值。如果不考慮摩擦系數較低的情況（μ≤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分布。由于散粒堆積體基本呈錐形，可以近似認為顆粒接觸方位在環向是對稱分布，因此這里只考慮顆粒接觸以豎直軸作為參考基準的方位分布規律。γ為顆粒接觸方位相對于豎直軸的方位角，取值范圍為0°～90°。從圖15可以看出，隨著γ的增大，顆粒接觸方位的概率密度f基本呈先增加后降低的趨勢，但概率密度峰值fmax的方位γ因摩擦系數不同而不同。μ= 0.01 時，fmax出現在γ= 35～40°；μ= 0.05時，fmax出現在γ= 30～35°；μ= 0.10 時，fmax出現在γ= 20～30°；μ= 0.30～2.00 時，峰值出現在γ=15～25°的時候。可以看出，峰值方位角隨著摩擦系數增大而降低，但當μ> 0.30時，峰值出現的方位（γ= 15～25°）不隨顆粒摩擦系數的增加而顯著改變，且該方位范圍與法向接觸力投影峰值的方位（θ≈15°）以及切向接觸力投影峰值的方位（θ≈25°）非常接近（見圖14）。需指出的是，這里的θ方位亦以豎直軸為參考。這些觀察充分證明了在顆粒摩擦系數較高（μ ≥0.3）的情況下，散粒堆積體內部拱效應的優勢發揮方位出現在偏離豎直軸15°～25°的方位。

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

2.5 堆積體中力鏈統計

圖16 統計了堆積體中弱力鏈的出現概率。弱力鏈是指接觸力小于平均接觸力的力鏈。Pn是指法向接觸力弱力鏈的概率，Pt是指切向接觸力弱力鏈的概率。由圖16 可知，Pn隨著顆粒摩擦系數的增加而微弱增加，并在μ大于0.10后，維持在0.68 左右。Pt在摩擦系數較低時接近1.0，而在μ大于0.30 后，隨著μ的增加而降低。可以發現，不論是法向接觸力還是切向接觸力，堆積體內弱力鏈的數量總是多于強力鏈的數量；切向接觸力中弱力鏈的占比要高于法向接觸力中弱力鏈的占比；隨著摩擦系數的增大，部分切向接觸力的弱力鏈會逐漸轉化為強力鏈。

圖16 弱力鏈統計Fig.16 The statistics of weak force chains

圖17 給出了μ= 0.05，0.10，0.30，0.50 時堆積體的力鏈網絡圖。圖中力鏈粗細（或顏色）表征接觸力的大小。可以發現，平均接觸力隨著μ增加而增大。當μ= 0.05，0.10 時，散粒堆積體中的強力鏈數量偏少，且呈水平攤開的分布，難以觀察到成“拱”現象。從圖17 中可見，當μ ≥0.3 時，隨顆粒摩擦系數的增大，強力鏈在水平和高度方向上逐漸擴展，強力鏈也逐漸與堆積體中心軸呈一定傾斜角度，成“拱”效應也越來顯著。