2維任意形狀顆粒間的滾動阻力模型

摘 要：本文將2維圓形顆粒間的滾動阻力模型推廣到非圓形顆粒體系，用于模擬準靜態(tài)剪切實驗，得到了滿足率無關(guān)和客觀性的剪切、剪脹曲線。給出的滾動阻力模型對Hertz接觸模型和線性接觸模型是普適的，改變顆粒形狀也是適用的，所得剪切、剪脹曲線和已有的實驗符合；作為對比，將文獻中幾種常用的滾動阻力模型用于相同的體系，單一速度下的剪切、剪脹曲線和實驗符合，但不同速度下的剪切、剪脹曲線卻不重合即不滿足準靜態(tài)顆粒體系下的率無關(guān)，且所得有效摩擦系數(shù)偏小。滾動阻力對剪切強度、剪脹-剪縮影響顯著，是必須考慮的相互作用，所給的任意形狀顆粒間的滾動速度滿足客觀性且其特殊情形和圓形顆粒間滾動速度一致。

關(guān)鍵詞：準靜態(tài)顆粒體系；非圓形顆粒；滾動阻力；率無關(guān)；客觀性

中圖分類號：O469

文獻標識碼： A

顆粒物質(zhì)在日常生活中無處不在，如農(nóng)業(yè)生產(chǎn)中的糧食、工業(yè)生產(chǎn)中的藥品、大氣污染中的粉塵、巖土工程中的土體、斷層間的斷層泥等等。顆粒物質(zhì)，即便是最簡單的無粘性的干顆粒體系，表現(xiàn)出來的力學(xué)性質(zhì)也十分復(fù)雜[1]。單個顆粒的微觀結(jié)構(gòu)如何影響體系的宏觀性質(zhì)也知之甚少[2]。緊致、準靜態(tài)顆粒流性質(zhì)的研究有著十分重要的意義，如體系從堵塞態(tài)到流動的轉(zhuǎn)變以及剪切過程中的滯滑和擾動與斷層失穩(wěn)、地震觸發(fā)有著密切聯(lián)系[3]。堵塞態(tài)到剪切帶的形成過程中力學(xué)性質(zhì)和結(jié)構(gòu)的演化一直是研究的重點，如用二維光彈材料對剪切體系中力鏈和接觸網(wǎng)絡(luò)的各向異性及力鏈結(jié)構(gòu)的長程相關(guān)性的研究[4-5]，及磁共振和X射線成像得到的三維環(huán)剪體系內(nèi)顆粒微觀結(jié)構(gòu)和速度場之間的關(guān)系[6]，還有根據(jù)直剪過程中聲速的變化來探測剪切帶的形成[7]等等。

用實驗的手段探測剪切流內(nèi)部結(jié)構(gòu)需用光彈材料[4-5]或磁共振和X射線成像[6]，費用較大，特別是對于三維顆粒體系，

獲取顆粒間接觸力和配位數(shù)等微觀量也是比較困難的，而離散元方法可以很方便地獲取實驗上難以取得的一些微觀量，成為研究顆粒體系從堵塞態(tài)到流化相變過程的有效手段。1979年Cundall和Strack關(guān)于離散元的工作[8]發(fā)表之后，離散元法就成為模擬顆粒體系的有效工具，得到了迅速的發(fā)展。顆粒間的接觸模型也由最直觀和簡單的線性模型[9-10]發(fā)展為復(fù)雜些和更接近理論的Hertz-Mindlin-Deresiewicz模型[11]。由于模型的復(fù)雜性，離散元中多用基于它的簡化模型，如半自鎖的法向力模型等等[12]；還有更為復(fù)雜的黏彈塑性滯回模型來模擬真實的加卸載過程[13]。同時，發(fā)現(xiàn)與球形或圓形顆粒相比橢球形或橢圓形顆粒更能真實反映顆粒物質(zhì)特性[14]，3維的非球形或2維的非圓形顆粒的離散元模擬通常用“multi-sphere（MS）models”法，即用相互重疊的不同半徑的球形或圓形組合出任意形狀的顆粒，來模擬真實的體系[12，15-16]。

非球形或非圓形顆粒體系的離散元模擬直接用球形或圓形顆粒間的法向和切向接觸模型來表示，只需將相應(yīng)的半徑換成接觸點的曲率即可[11，17]。然而，球形或圓形顆粒間滾動阻力模型卻一直有著不同的甚至相互矛盾的模型。用來計算滾動阻力的模型有

這些不同的滾動阻力模型的提出都有一定的理論依據(jù)和清晰的數(shù)學(xué)定義，用來研究不同的體系都有一定的合理性，但不同文獻對滾動速度的定義卻是不同的。如對于球形或圓形顆粒間的滾動速度，文獻[13]vr=-aijn×（ωi-ωj），而文獻[24]vr=-aijn×ωi-ωj-0.5（ai-aj）/ai+aj）vt，后者的定義比前者多了和相對平動速度有關(guān)的項。對于顆粒氣體或者沙漏、堆積等問題，顆粒間的滾動阻力對體系的影響要比準靜態(tài)的緊致顆粒體系要弱，用不同的滾動阻力模型也可以得到和實驗吻合較好的模擬結(jié)果[11-12，17-20]。但對于直剪體系，

顆粒間的滾動在剪切帶的形成中起重要作用[21-23]，不同的滾動阻力模型對模擬結(jié)果影響是顯著的，對于球形顆粒的剪切系統(tǒng)，不同的滾動阻力都能得到和實驗吻合的剪切曲線，但有些滾動阻力模型卻不滿足準靜態(tài)體系的率無關(guān)，即不同速度的剪切曲線不能重合[25]。正確滾動阻力模型引入前提是要有正確的滾動速度，文獻[13]中所用的滾動速度是符合客觀性，用于模擬滾動阻力是滿足率無關(guān)的[26]，但文獻[13]給出的只是球形或圓形顆粒間的滾動速度，非球形或非圓形顆粒間的滾動速度雖然文獻[24]中給出了但卻是不恰當?shù)模愿鶕?jù)非球形或非圓形顆粒間滾動速度來確定滾動力矩的離散元模型很少有人涉及?，F(xiàn)有的非球形或非圓形顆粒直剪模擬或用商業(yè)軟件如PFC3D[27]，或直接將顆粒間滾動阻力忽略[28]，非球形或非圓形顆粒直剪體系滾動阻力模型討論的也較少。其可能原因是球形顆粒的滾動定義尚未統(tǒng)一，非球形顆粒間滾動速度則更為復(fù)雜。

本文將文獻[13]中圓形顆粒間滾動速度向2維任意形狀推廣得到非圓形顆粒間的滾動速度，用以計算非圓形顆粒間的滾動阻力，對非圓形顆粒剪切體系進行了離散元模擬，得到了符合實驗結(jié)果[29-32]的滿足率無關(guān)的剪切、剪脹曲線。所得出的任意形狀間的滾動速度滿足客觀性，滾動速度表達式對3維任意形狀顆粒也適用。由于本文涉及到的模型較多，含大量的公式和變量，限于篇幅不能每個公式及變量作詳細的解釋，僅給出簡要的說明，列出了公式的出處供進一步了解。

1 非圓形顆粒的離散元模型

1.1 非圓形顆粒的生成及剪切體系

離散元模擬中任意形狀顆粒的描述是一個難點[33]，大多數(shù)離散元模擬采用的是圓形或球形顆粒，非圓形或非球形顆粒也多采用能用解析方程表示的形狀如二維的橢圓三維的橢球[11-12，14]。然而真實顆粒的形狀很少能用簡單的解析方程來描述，用解析的方法判斷接觸也是極為困難的，本文用五個相互重疊的半徑不同的圓來描述橢圓顆粒[33]，所逼近的橢圓方程為（x/a）2+（y/b）2=1，a=1.5，b=0.5a，0.6a，0.7a，0.75a，0.85a，0.95a，如圖1所示。這種方法可以通過增加圓數(shù)目、改變圓半徑、擺放圓心位置來逼近任意形狀的顆粒。用組合法生成的非圓形顆粒間的接觸完全可以由圓形顆粒間的接觸關(guān)系來確定。這樣做的代價是增大了存儲和計算量。

均勻混合六種顆粒組成單剪體系如圖2所示，上下兩層用顆粒做邊界。下邊界在制備和剪切過程中保持固定，對上邊界施加正應(yīng)力σzz，體系穩(wěn)定后在水平方向施加剪切速度vx，中間為受壓剪切顆粒，水平方向采用周期性邊界條件。

1.2 顆粒間法、切向力

1.3 顆粒間的滾動

顆粒間的滾動阻力有多種形式，文獻[25]用黏彈塑性滯回模型（adhesive-elasto-plastic contact model）計算法、切向力，用球形單剪顆粒體系對五種滾動阻力模型進行了驗證，確定了五種模型中唯一一種符合客觀性和率無關(guān)的滾動阻力模型，即顆粒間的相對滾動速度取vr=-aijn×（ωi-ωj），滾動阻力和滾動阻力矩用類似與切向力的方式確定，只需要將切向滑動速度vt改為滾動速度vr，滑動摩擦系數(shù)改為滾動摩擦系數(shù)μd，滾動阻力矩qrolling=aijn×fr，aij=aiaj/（ai+aj）為兩顆粒的約化半徑。此滾動阻力模型對接觸模型也是普適的，將此模型用于線性和Hertz接觸模型得到的剪切、剪脹曲線都符合率無關(guān)（參見第2節(jié)）。

2 滾動對剪切曲線的影響

用單剪顆粒體系（圖2）驗證公式（1）、（2）所給的滾動阻力模型對線性和Hertz接觸模型的普適性，選用表1中的參數(shù)，計算滾動阻力對剪切曲線的影響。物理參數(shù)的選取大致是按照離散元模擬中常用的值，其實表1中值的大小并不是很重要，因為程序中計算用的是無量綱的數(shù)值，同一個程序采用不同的標準單位可以表示不同的材料，如果時間單位采用10-5 s，質(zhì)量單位用10-2 kg，長度單位取10-2 m，相同的程序表示密度為2.65×104 kg/m3，楊氏模量10 GPa的材料，可以任意選取時間、質(zhì)量和長度單位，來給出想要的材料參數(shù)值[34]。

圖4給出的是缺少滾動阻力的應(yīng)力-應(yīng)變和體應(yīng)變和剪應(yīng)變間的關(guān)系，圖4（a）縱坐標表示剪應(yīng)力和正應(yīng)力之比（也被稱作有效摩擦系數(shù)），橫坐標表示剪切位移和體系長度之比即剪應(yīng)變，圖4（b）縱坐標表示體積變化百分比。缺少滾動阻力不同速度下的剪切曲線是不對的，模擬時體系慣性數(shù)Ι=γ˙d/P/ρ<10-3滿足準靜態(tài)條件，其中d為顆粒體系粒徑尺寸（通常用平均半徑衡量），P為體系正應(yīng)力，ρ為體系平均密度，γ˙為剪切速率和體系高度之比，不同剪切速度的剪切曲線和剪脹曲線應(yīng)該是重合的，而缺少滾動阻力時剪切強度和剪脹卻表現(xiàn)出和速度的相關(guān)性。準靜態(tài)條件下的剪切體系，顆粒間的滾動是必須考慮的，和球形顆粒類似并非所有的滾動阻力模型都能夠滿足率無關(guān)，文獻[18]、[19]以及[11，17，20]中沙漏體系中選用的滾動阻力模型（μr″=0.3，其他參數(shù)數(shù)值不變）用于相同的體系得到不同速度下的剪切曲線如圖5-7所示。

三種滾動阻力模型不同速度的剪切曲線重合得不是太好，而且有效摩擦系數(shù)偏小，將本文所給的滾動阻力模型（公式（1）、（2））用于Hertz接觸模型，計算結(jié)果如圖8所示。

可見本文所給的滾動阻力模型對率無關(guān)的滿足要好于圖5-7所給的三種滾動阻力模型，而且有效摩擦系數(shù)可以和文獻[30-32]中的實驗結(jié)果對比，體應(yīng)變也較大（即剪脹比較明顯）。

將本文所給的滾動阻力模型用于線性的法、切向作用力，也得到了滿足率無關(guān)的剪切和剪脹曲線，結(jié)果如圖9所示。

改變顆粒形狀和剪切體系得到的剪切曲線也和實驗一致[30-32]，圖10是任意形狀顆粒直剪體系，圖11是應(yīng)力-應(yīng)變和體應(yīng)變隨剪應(yīng)變變化的計算結(jié)果：

3 討論和結(jié)論

對于顆粒氣體或者沙漏、斜槽流等稀疏體系，顆粒間的滾動作用相比于法向碰撞較弱，引入不同的模型對結(jié)果影響不大，主要起到數(shù)值穩(wěn)定作用[18]，而且對于較為稀疏的顆粒體系沒有一個類似率無關(guān)類似的必須滿足的條件，離散元模擬結(jié)果和實驗結(jié)果對比時缺少這樣一個檢驗標準，只要實驗關(guān)注的某個物理量和模擬符合較好，那么就認為模型是合理的；對于準靜態(tài)條件下的剪切顆粒體系，滾動阻力是不可以忽略的相互作用，不同的滾動阻力模型單一剪切速度下的曲線和實驗都是符合的，如果沒有率無關(guān)的限制模型的優(yōu)劣是無法區(qū)別的。

準靜態(tài)體系滿足率無關(guān)對滾動阻力模型優(yōu)劣的甄別提供了一個標準，正確的滾動阻力模型要滿足

以下條件：

第一，滾動力引入機制要正確，這一點現(xiàn)有的滾動阻力模型都有自己的機制，都有自己的合理性；

第二，滾動阻力的引入要滿足守恒律，如角動量的守恒等；

第三，滾動速度的引入要滿足客觀性，即顆粒間的滾動速度不隨觀察者和坐標系的改變而改變，換句話說，在不同的參考系里滾動速度有相同的表達式和值。

第四，要滿足準靜態(tài)顆粒體系的率無關(guān)；

第五，要和實驗結(jié)果符合；

第六，對不同的法向、切向力接觸模型要有普適性，對顆粒形狀和體系也要有普適性。

如果用這幾點來判斷滾動阻力模型優(yōu)劣，本文給出的模型滿足，所給的非圓形顆粒滾動速度vr=-ricnrjcn/（ricn+rjcn）n×（ωi-ωj）滿足客觀性，因為前面的系數(shù)是標量不隨坐標系改變，顆粒繞質(zhì)心旋轉(zhuǎn)，并不改變顆粒間的相對角速度，質(zhì)心系中接觸面法線方向單位矢量只是將坐標原點進行一個平移并不改變方向，所以滾動速度在質(zhì)心系和實驗室坐標系中有相同的表達式和值，滿足客觀性；顆粒間滾動阻力矩qrolling=ricnrjcn/（ricn+rjcn）n×fr保證了兩個顆粒受到的力矩大小相等方向相反，滿足角動量守恒；圖8-11計算結(jié)果滿足率無關(guān)和實驗相符[29-32]，對接觸模型和形狀也是普適的，滾動速度率無關(guān)的論證涉及變量和公式較多已在文獻[25]和[34]中論證了，在此不再列出，具體證明可以參考此兩篇文獻。文獻[18]、[19]中的模型qrolling=-urRiFnijωi、qrolling=-u′rFnijωi因為兩個顆粒在實驗室中的角速度方向并不時時刻刻一致，導(dǎo)致計算出來的滾動力矩對兩個顆粒并不相等，角動量不守恒，從這一點而言模型是不恰當?shù)?；qrolling=-u′rFnijωij [11，17，21]滿足角動量守恒和客觀性，但只是加了一個和顆粒間相對滾動方向相反的力矩，顆粒間的滾動阻力除了和相對滾動方向有關(guān)外還應(yīng)該和相對滾動位移大小有關(guān)，本文所給的模型考慮了相對滾動形變的大小和方向，還考慮了法向相互作用力和滾動力之間的庫倫條件。

本文在形式上對圓形顆粒間滾動速度進行推廣，給出了非圓形顆粒間的相對滾動速度和滾動阻力力矩，得到了符合率無關(guān)且和實驗相符的剪切曲線，是一個合適的滾動阻力模型。

參考文獻：

[1]MIdI G.On dense granular flows[J].European Physical Journal E Soft Matter，2004，14（4）：341-365.

[2]ODA M，KAZAMA H.Microstructure of shear bands and its relation to the mechanisms of dilatancy and failure of dense granular soils[J].Geotechnique，1998，48（4）：465-481.

[3]MELOSH H J.Dynamical weakening of faults by acoustic fluidization[J].Nature，1996，379（6566）：601-606.

[4]MAJMUDAR T S，BEHRINGER R P.Contact force measurements and stress-induced anisotropy in granular materials[J].Nature，2005，435（7045）：1079-1082.

[5]ZHANG J，MAJMUDAR T S，TORDESILLAS A，et al.Statistical properties of a 2D granular material subjected to cyclic shear[J].Granular Matter，2009，12（2）：159-172.

[6]MUETH D M，DEBREGEAS G F，KARCZMAR G S，et al.Signatures of granular microstructure in dense shear flows[J].Nature，2000，406（6794）：385-389.

[7]KHIDAS Y， JIA X.Probing the shear-band formation in granular media with sound waves[J].Physical Review E，2012，85（5）：051302.1-051302.6.

[8]CUNDALL P A，STRACK O D L.A discrete numerical model for granular assemblies[J].Géotechnique1979，29（1）：47-65.

[9]IMOLE O I，WOJTKOWSKI M，MAGNANIMO V，et al.Micro-macro correlations and anisotropy in granular assemblies under uniaxial loading and unloading[J].Physical Review E，2014，89（4）：042210.1-042210.22.

[10]GONZLEZ S，WINDOWSYULE C R，LUDING S，et al.Forced axial segregation in axially inhomogeneous rotating systems[J].Physical Review E，2015，92（2）：022202.1-022202.9.

[11]LIU S D，ZHOU Z Y，ZOU R P，et al.Flow characteristics and discharge rate of ellipsoidal particles in a flat bottom hopper[J].Powder Technology，2014，253：70-79.

[12]VU-QUOC L，ZHANG X，WALTON O R.A 3-D discrete-element method for dry granular flows of ellipsoidal particles[J].Computer methods in applied mechanics and engineering，2000，187（3）：483-528.

[13]LUDING S.Cohesive，frictional powders：contact models for tension[J].Granular matter，2008，10（4）：235-246.

[14]ROTHENBURG L，BATHURST R J.Numerical simulation of idealized granular assemblies with plane elliptical particles[J].Computers and Geotechnics，1991，11（4）：315-329.

[15]FAVIER J F，ABBASPOUR-FARD M H，KREMMER M，et al.Shape representation of axi-symmetrical，non-spherical particles in discrete element simulation using multi-element model particles[J].Engineering Computations，1999，16（4）：467-480.

[16]JENSEN R P，BOSSCHER P J，PLESHA M E，et al.DEM simulation of granular media—structure interface：effects of surface roughness and particle shape[J].International Journal for Numerical and Analytical Methods in Geomechanics，1999，23（6）：531-547.

[17]ZHOU Z Y，PINSON D，ZOU R P，et al.Discrete particle simulation of gas fluidization of ellipsoidal particles[J].Chemical Engineering Science，2011，66（23）：6128-6145.

[18]YANG R Y，ZOU R P，YU A B.Computer simulation of the packing of fine particles[J].Physical Review E，2000，62（3）：3900-3908.

[19]ZHOU Y C，WRIGHT B D，YANG R Y，et al.Rolling friction in the dynamic simulation of sandpile formation[J].Physica A：Statistical Mechanics and its Applications，1999，269（2）：536-553.

[20]ZHU H P，YU A B.Averaging method of granular materials[J].Physical Review E，2002，66（2）：021302.1-021302.10.

[21]IWASHITA K，ODA M.Rolling resistance at contacts in simulation of shear band development by DEM[J].Journal of Engineering Mechanics，1998，124（3）：285-292.

[22]IWASHITA K，ODA M.Micro-deformation mechanism of shear banding process based on modified distinct element method[J].Powder Technology，2000，109（1）：192-205.

[23]SALAZAR A，SEZ E，PARDO G.Modeling the direct shear test of a coarse sand using the 3D Discrete Element Method with a rolling friction model[J].Computers Geotechnics，2015，67：83-93.

[24]BAGI.A definition of particle rolling in a granular assembly in terms of particle translations and rotations[J].Journal of Applied Mechanics，2004，71（4）：493-501.

[25]ZHAO C，LI C.Influence of rolling resistance on the shear curve of granular particles[J].Physica A Statistical Mechanics Its Applications，2016，460：44-53.

[26]WANGY，ALONSO-MARROQUIN F，GUO W W， et al.Rolling and sliding in 3-D discrete element models[J].Particuology，2015，23（6）：49-55.

[27]HRTL J，JIN Y O.Numerical investigation of particle shape and particle friction on limiting bulk friction in direct shear tests and comparison with experiments[J].Powder Technology，2011，212（1）：231-239.

[28]ZHAO S，ZHOU X，LIU W.Discrete element simulations of direct shear tests with particle angularity effect[J].Granular Matter，2015，17（6）：793-806.

[29]DEPKEN M，VAN S W，VAN H M.Continuum approach to wide shear zones in quasi-static granular matter[J].Physical Review E，2006，73（3）：031302.1-031302.11.

[30]CHEN Q，WANG Q H，ZHAO C，et al.Mechanical response study of glass-rubber particle mixtures[J].Acta Physica Sinica，2015，64（15）：330-338.

[31]ZHANG Q，LI Y C，LIUR，et al.Acoustic probing of the granular solid system under direct shear[J].Acta Phys Sinica，2012，61（23）：234501.1-234501.8.

[32]ZHANG Q，LI Y，HOU M，et al.Elastic waves in the presence of a granular shear band formed by direct shear[J].Physical Review E，2012，85（3）：031306.1-031306.6.

[33]MARKAUSKAS D，KAIANAUSKAS R，DIUGYS A，et al.Investigation of adequacy of multi-sphere approximation of elliptical particles for DEM simulations[J].Granular Matter，2010，12（1）：107-123.

[34]ZHAO C，LI C，HUL.Rolling and sliding between non-spherical particles[J].Physica A，2018，492：181-191.

（責任編輯：周曉南）