高政國 董朋昆 張雅俊 孫卉竹 迪 亞

* (北京航空航天大學交通科學與工程學院,北京 100191)

? (中國礦業(yè)大學(北京)力學與建筑工程學院,北京 100083)

引言

滾動阻力問題涉及車輛工程[1],土木工程[2]和農(nóng)業(yè)[3-4]等許多領域.與滑動阻力相比滾動阻力通常較小[5],如硬質(zhì)圓形或球型顆粒相對滾動時滾動摩擦系數(shù)一般為10?5～ 10?3,滑動摩擦系數(shù)通常為0.1～ 1.從數(shù)值來看,滾動摩擦系數(shù)比滑動摩擦系數(shù)要小得多,但是滾動摩擦對顆粒的宏觀力學特性有著非常重要的影響.Bardet 等[6]在離散元模型中最早考慮了滾動約束的作用,他們發(fā)現(xiàn)不考慮滾動約束作用,離散元模擬結(jié)果得到的力學參數(shù)在理論值范圍之外.他們認為接觸力偶矩可能起著重要作用,他們指出當約束顆粒滾動自由度后顆粒體系的內(nèi)摩擦角增大.Morgan 等[7]直接將轉(zhuǎn)動阻尼引入離散元模型進行斷層泥的模擬,他們得到了與實驗室評估接近的結(jié)果.Sakaguchi 等[8]將“滾動摩擦”的概念引入離散元模型,進行谷倉清空過程顆粒流阻塞試驗與數(shù)值模擬研究.他們在離散元計算程序中引入了一個滾動摩擦力矩,發(fā)現(xiàn)常規(guī)離散元數(shù)值模擬中圓盤顆粒成拱不穩(wěn)定,易于破壞,考慮滾動摩擦后能有效的模擬出物理試驗中得到的阻塞成拱現(xiàn)象.

目前在滾動阻力理論及工程應用方面許多學者開展了相關研究[9-16].然而,顆粒體系穩(wěn)定過程中滾動阻力作用機制仍不十分清楚,與材料[17-18]及滾動阻力有關的許多力學模型和方法仍有待研究.

滾動阻力通常被學者們稱為“滾動摩擦”.根據(jù)摩擦學理論,滾動阻力的主要來源是接觸表面上的微小滑移、塑性變形、材料的黏滯性、表面附著力和形狀效應[19].

1875 年,Reynolds[20]發(fā)現(xiàn),當金屬圓柱體在橡膠表面上滾動時,在豎向壓力作用下接觸面的切向位移會發(fā)生微小差異,稱為微滑移.類似地,Tabor 等[21]研究了彈性范圍內(nèi),硬質(zhì)球體和硬質(zhì)圓柱體在橡膠軟基上滾動時的摩擦作用.他們發(fā)現(xiàn)接觸面上的切向力始終小于在滾動過程中存在潤滑劑的情況下產(chǎn)生的界面滑移值.因此,他們認為微滑移并非彈性范圍內(nèi)滾動摩擦的主要原因.為了解釋這種情況下的滾動特性,他們提出了彈性滯后是滾動摩擦的主要原因.

Flom 和Bueche[22]采用具有黏滯特性的Voigt模型對彈性范圍內(nèi)材料的彈性滯后進行了研究.根據(jù)該模型,當堅硬球體或圓柱體在軟基上以一定速度滾動時,接觸面的后緣將與軟基脫離形成接觸面壓力的不均勻分布.而滾動阻力力矩正是源自這種不對稱分布的接觸壓力.他們認為,相對堅硬的球體在較軟的基材上的滾動摩擦力會隨滾動速度的變化而改變.對于黏彈性材料的滾動顆粒,滾動阻力主要是由接觸面上體積變形產(chǎn)生的能量耗散引起的,而與表面附著力的關系較小[23].在該研究中,彈性滯后模型被當作具有黏彈特性的力學模型,是速度相關型的.當滾動速度非常小時,滾動阻力接近于零.同樣的,對于在硬質(zhì)基體材料上滾動的軟黏彈性球,滾動摩擦阻力取決于滾動速度,當滾動速度為零時,滾動阻力也將為零[24].但是,當采用這種速度相關型黏彈性模型研究顆粒堆積穩(wěn)定時,若顆粒處于穩(wěn)定狀態(tài),滾動阻力將會消失.顯然,這種穩(wěn)定機制是有問題的.當滾動速度為零時,由于阻力的消失顆粒體系穩(wěn)定性會降低,較顆粒滾動速度不為零時更容易崩塌是不合理的.

Greenwood 等[23]對橡膠薄壁管進行了扭轉(zhuǎn)和拉伸共同作用試驗,來研究材料的彈性滯后現(xiàn)象.試驗結(jié)果顯示,橡膠具有與加載速度無關的彈性滯后.當在彈性范圍內(nèi)進行加載和卸載時,由于應變變化落后于應力,使應力應變曲線的加載和卸載過程不一致,形成了閉環(huán),進而產(chǎn)生了能量耗散.因此,對于彈性材料,彈性滯后引起的滾動阻力包括速度相關型的和速度無關型的兩部分.

1979年,Cundall 和Strack[25]提出了離散元方法(DEM)并迅速引起關注[26-27],它可以方便地再現(xiàn)物理試驗過程細節(jié),并且能有效降低試驗成本,目前已成為一種被廣泛接受的數(shù)值計算方法.離散元模型中滾動阻力力學模型的建立在顆粒體系力學行為模擬中至關重要.

為了研究在剪切帶試驗中觀察到的顆粒間的巨大空隙和高旋轉(zhuǎn)梯度,Iwashita 和Oda[28]在DEM 模型中考慮了滾動阻力,建立了改進的離散元模型(MDEM).該模型滾動阻力由轉(zhuǎn)動方向的彈簧、黏壺、非承拉節(jié)點和摩擦型滑阻器元件表達,其中由彈簧和黏壺組成的Voigt 模型,是一種速度相關型的阻力元件.而速度無關型的滾動阻力通常由摩擦型滑阻器元件表達.研究結(jié)果表明該模型可以有效預測剪切帶的剪脹行為.MDEM 是目前應用比較成功的離散元模型,在此基礎上,許多學者開展了相關的計算應用與模型改進工作[9,29-30].在本文稱MDEM 模型為常規(guī)DEM 模型.

目前對于滾動阻力與滾動速度之間的相關性仍沒有統(tǒng)一的認識,在離散元模擬中通常采用兩種典型的滾動阻力公式[31-33]:第一種滾動力矩與滾動速度無關,力矩大小正比于法向接觸力,與滾動方向相反.當顆粒間接觸力恒定時,滾動力矩是一個恒定值;第二種滾動力矩與相對角速度成正比,表現(xiàn)為一種黏滯力特征.離散元模擬結(jié)果顯示按照速度相關型黏滯滾動阻力公式不能很好地模擬顆粒的穩(wěn)定堆積.而單獨按照速度無關型恒定的滾動阻力公式雖然模擬顆粒堆積穩(wěn)定性有所提高,但顆粒臨近穩(wěn)定靜止時會在平衡位置往復振動,此時滾動阻力大小不變,方向正負變化,微觀上無法達到靜止狀態(tài)[19].因此,在離散元算法中,摩擦型滑阻器元件表達的速度無關型的滾動阻力當顆粒臨近靜止時要退出工作,只有速度相關型黏滯力元件工作.所以,在顆粒堆積問題離散元模擬中,彈性滯后引起的速度無關型滾動阻力是建立顆粒臨近靜止狀態(tài)顆粒接觸力學模型的一個重要思路.

直接采用滾動阻力公式表達的力學元件建立離散元模型應用方便,但由于機理上缺乏深入認識,一些滾動阻力參數(shù)確定往往通過經(jīng)驗或試算得到.滾動阻力通常較小,通過試驗方法直接識別滾動阻力參數(shù)的難度較大[34].

筆者提出了通過圓形顆粒在剛性平面上滾動停止前的往復擺動現(xiàn)象測量滾動阻力參數(shù)的方法,研發(fā)了一個顆粒微動力光學試驗檢測系統(tǒng)[35],可通過測量顆粒的往復擺動曲線識別Voigt 模型的轉(zhuǎn)動剛度與黏滯阻尼參數(shù).試驗研究發(fā)現(xiàn)通過識別參數(shù)的常規(guī)DEM 模型計算得到的顆粒擺動位移曲線與試驗曲線在臨近靜止時刻吻合程度劇烈下降,計算得到的顆粒穩(wěn)定時間較試驗長,且試驗曲線在臨近靜止時刻出現(xiàn)擺動頻率的改變,常規(guī)DEM 模型不能從機理上解釋這一現(xiàn)象.

為此,本文基于彈性滯后理論研究建立了一種滯彈簧力學元件,將與速度無關的材料彈性滯后特性引入,提出一種滯彈性滾動阻力模型,以此建立對試驗中顆粒臨近靜止狀態(tài)滯彈性耗能機理的解釋.與傳統(tǒng)DEM 模型相比,改進后的滯彈性滾動阻力模型計算結(jié)果與試驗結(jié)果更為符合,驗證了滯彈簧滾動阻力模型的有效性.

1 彈性滯后

彈性范圍內(nèi)材料加載卸載時,應變往往落后于應力,使得應力?應變加載線與卸載線不重合圍成一封閉回線,形成彈性滯后現(xiàn)象.與速度相關的彈性滯后現(xiàn)象通常表達為彈性材料的黏性行為,而與速度無關的彈性滯后源于材料加載與卸載過程應力應變不能同步,這一現(xiàn)象的主要原因在于材料分子間的相互作用和弛豫時間,與加載速度無關.通過薄壁橡膠管的純拉伸試驗,Tabor[21]獲得應力?應變滯回曲線.當外力沒有達到極限應變卸載時,應力?應變曲線上將形成一個轉(zhuǎn)折點(εrev1,σrev1),如圖1 所示.當卸載不完全,再次加載時,應變將在最后一次卸載的終點(εrev2,σrev2)繼續(xù)加載.定義彈性滯后上升曲線為加載曲線,下降曲線定義為卸載曲線,并以(εe,σe)表示彈性極限點.因此,可以通過如下冪指數(shù)表達式定義應力和應變之間的關系

圖1 彈性滯后示意圖Fig.1 Schematic diagram of elastic hysteresis

加載

當應變ε∈ (0,εe),表示材料處在彈性范圍.定義參數(shù)β∈ (0,1),表示應變滯后于應力的程度.β值越接近0,加載曲線和卸載曲線所包圍的面積越大,彈性滯后所引起的能量耗散也就越大.

為定義在復雜應力下的彈性滯后的能量耗散過程,可根據(jù)加載和卸載構(gòu)造相應的能量耗散過程.以單個完整加載或單個完整卸載定義為一個子過程,將整個荷載過程分為多個子過程組合,可表達為

Δε是不足一個完整子過程的多余應變.子過程的能量密度表達式可寫為

加載e表示每個荷載過程的能量密度,表示不能構(gòu)成一個子過程的多余應變的能量密度.

2 滯彈簧與HDEM 模型

根據(jù)接觸方向,常規(guī)DEM 接觸模型[28]可分為法向接觸模型,切向接觸模型和滾動阻力模型,如圖2所示.滾動方向阻力模型由滾動彈簧、滾動黏壺、摩擦器和非承拉節(jié)點組成.滾動模型所提供的滾動阻力可表示為

圖2 常規(guī)DEM 模型Fig.2 DEM model

式中Kr是彈簧剛度系數(shù),cr是滾動阻尼系數(shù),μr是滾動摩擦系數(shù).

從上式可看出,顆粒發(fā)生持續(xù)滾動時,滾動力矩總是等于最大摩擦力矩μrFn,當顆粒滾動不能持續(xù)時,滾動阻力小于最大摩擦力矩μrFn,常規(guī)DEM 模型所提供的阻力值與滾動角速度有關.滾動角速度越大,滾動模型所提供的阻力值越大.滾動彈簧是線彈性的,不能耗散能量,動能的耗散僅依靠黏壺黏滯力做功實現(xiàn).為表征滾動摩擦中速度無關的摩擦力,根據(jù)上述建立彈性滯后的應力?應變關系,提出圖3所示的滯彈簧元件.

滯彈簧的角位移與彈性力不同于常規(guī)DEM 滾動模型中滾動彈簧的線彈性關系.參照彈性滯后應力應變表達,滯彈簧的加載和卸載遵循以下關系

加載

式中Δ表示滯彈簧的位移變形,Δe是滯彈簧彈性變形量的極限值,F表示滯彈簧的恢復力值,Fe是彈簧彈性力的極限值.

滯彈簧耗散的能量可以表示為

Esum表示整個滯彈簧運動過程中由于彈性滯后效應而耗散的能量,Δrev表示加載和卸載過程中滯彈簧轉(zhuǎn)折點的位移值,Frev表示加載和卸載過程中滯彈簧轉(zhuǎn)折點的力值.

當滾動顆粒在平衡位置往復擺動時,滾動過程可分為正向加載、正向卸載、負向加載、負向卸載四個階段.根據(jù)滯彈簧變形與滾動角速度,式(11)和式(12)計算卸載與卸載時滯彈簧元件的彈性恢復力,負向時彈性恢復力取負.

將滯彈簧與黏壺、摩擦器及非承拉節(jié)點等元件進行組合,提出一種新的滯彈性滾動阻力離散元模型(HDEM),如圖4 所示.滯彈簧可以表征顆粒堆積穩(wěn)定過程中與運動速度無關的滾動阻力耗能.

圖4 HDEM 模型Fig.4 HDEM model

3 驗證

3.1 滯彈簧有效性驗證

當一個在剛性平面上的滾動圓柱試件速度減小到一定程度后,會在一個平衡位置往復擺動.由于動能耗散,擺動幅度逐漸減小直至為零,如圖5 所示.常規(guī)DEM 滾動阻力模型參數(shù)可以通過測量試件的擺動來識別[36].

圖5 顆粒自由滾動示意圖Fig.5 Particle free-rolling on a flat surface

基于常規(guī)DEM 模型,自由滾動的圓柱試件的運動平衡方程為

其中,θ為滾動角位移,Ks和Kr分別表示切向和滾動彈簧剛度系數(shù),cs和cr分別表示切向和滾動方向阻尼系數(shù),Jz是對接觸點的轉(zhuǎn)動慣量,Jc是對顆粒形心處的轉(zhuǎn)動慣量,Fx是水平慣性力,Mr是慣性力矩,R為圓柱試件半徑.

方程(18)為有阻尼振動方程,可得到滾動角位移表達式為

其中ω為振動圓頻率,ξ為阻尼比,A為最大幅值,α為相位角.而阻尼系數(shù)和滾動剛度可表達為

因此,只需測得顆粒的往復擺動曲線,按照式(21)～式(23)可識別出振動圓頻率ω和阻尼比ξ,進而識別出阻尼系數(shù)cr和滾動剛度Kr.

采用激光位移傳感器微振動位移檢測試驗裝置,可試驗測得圓柱體試件往復擺動曲線,試驗裝置如圖6.

圖6 試驗裝置圖Fig.6 Experimental device diagram

通過檢測試驗,我們可以得到圓柱滾動的時間歷程曲線,如圖7 所示.

圖7 圓柱滾動角位移時程曲線Fig.7 Angular displacement variation versus time for a rolling cylinder

為驗證滯彈簧元件有效性,分別使用常規(guī)DEM 滾動阻力模型與HDEM 模型對顆粒的純滾動過程進行了離散元模擬,如圖8 所示.圖9 是圓柱體在平衡位置擺動過程的數(shù)值模擬結(jié)果.圖中常規(guī)DEM 模型中彈簧角位移與彈性力的線性變化關系,而HDEM 模型中滯彈簧的角位移與彈性力形成滯回環(huán),與所建立的位移?荷載公式一致.

圖8 顆粒轉(zhuǎn)動模型示意圖Fig.8 Particle rotation model

圖9 滾動過程Fig.9 Rolling process

圖10 中黑色曲線是采用激光位移傳感試驗測得的滾動彈簧剛度值,并使用常規(guī)DEM 模型進行離散元數(shù)值模擬得到的動能衰減包絡線.紫色曲線是將常規(guī)DEM 模型中阻尼系數(shù)放大1.5 倍后的動能衰減包絡線,藍色曲線是HDEM 模型模擬出的動能衰減包絡線.在擺動起始時刻將黏壺的阻尼系數(shù)設置為0,可以看出將阻尼系數(shù)放大后,動能衰減曲線近似向下平移.而HDEM 模擬的動能衰減曲線在起始時刻7 s 時與藍色曲線交匯,交匯前位于藍色曲線上方,交匯后位于下方,說明在臨近靜止過程中滯彈簧的耗能大于黏壺.

圖10 HDEM 與DEM 模擬動能變化Fig.10 Kinetic energy evolution with HDEM and DEM simulation

參數(shù)β能反映滯彈簧的耗能能力,其值越小耗能能力越強.通過與試驗數(shù)據(jù)進行對比,經(jīng)過試算擬合可得β值.對10 個不同材質(zhì)圓柱形試件測量結(jié)果擬合得到的β值如圖11 所示.聚氨酯圓柱β平均值為0.844,鋁圓柱β平均值為0.874.聚氨酯的平均值要小于鋁的平均值,分析原因為聚氨酯材料質(zhì)地較軟,在彈性滯后過程中會產(chǎn)生更多的能耗.

圖11 橡膠圓柱與鋁圓柱的β 值Fig.11 β values for the rubber and aluminum cylinder

雖然材料的彈性滯后耗能同樣存在法向與切向.但由于離散元接觸模型中法向與切向剛度比轉(zhuǎn)動方向剛度大得多,法向與切向振動頻率高,速度相關型的黏壺元件耗能要比滯彈簧耗能大得多,因此HDEM 模型中慮法向、切向滯彈簧與轉(zhuǎn)動方向滯彈簧相比作用不顯著,在模型中可不考慮使用.

3.2 耗能分析

圖12 為橡膠圓柱自由滾動激光位移傳感器試驗、常規(guī)DEM 模型和HDEM 模型離散元數(shù)值模擬得到的時間?相對位移圖.從圖中可以看出,常規(guī)DEM 和HDEM 在振蕩早期的動能衰減差異不大.但從整體上,常規(guī)DEM 模型計算得到的滾動到靜止時間較試驗結(jié)果要長,顯示常規(guī)DEM 模型不易達到靜止穩(wěn)定;HDEM 模型與試驗結(jié)果更為接近.因此可以得出HDEM 模型在接近靜止時具有更強的能量耗散能力,與試驗結(jié)果吻合更好.試驗曲線不光滑是由于儀器采集誤差造成的.

圖12 DEM 滾動模型與HDEM 模型模擬結(jié)果Fig.12 Rolling angle versus time for DEM and HDEM model

HDEM 模型中有滯彈簧和黏壺兩個耗能元件.為分析速度無關滾動阻力與速度相關滾動阻力關系,提取圖12 中數(shù)值模擬結(jié)果,得到體系總動能變化如圖13(a)與黏壺作用耗能變化如圖13(b).可以看出,體系動能不斷衰減,黏壺能量耗散逐漸減小.

圖13 體系動能與黏壺耗能Fig.13 Kinetic energy and energy dissipated by the damper

圖14(a)是HDEM 模型中滯彈簧耗能與黏壺耗能變化.圖14(b)是他們的比值隨時間變化趨勢.可以看出隨著滾動速度降低,HDEM 模型中滯彈簧元件所代表的與速度無關的能量耗散比例越來越大.

圖14 滯彈簧與黏壺耗能Fig.14 Energy consumption ratio of hysteresis spring and damper

3.3 頻率擬合結(jié)果

選取4 組橡膠圓柱和鋁圓柱試件的試驗檢測數(shù)據(jù)和常規(guī)DEM 模型、HDEM 模型計算結(jié)果數(shù)據(jù)進行對比(如圖15 和圖16).由于常規(guī)DEM 模型表達的振動具有頻率不變特性,常規(guī)DEM 模型與觀測試驗結(jié)果對比可以看出,擺動速度接近0 時,試驗中圓柱體擺動頻率有增大現(xiàn)象;HDEM 模型的數(shù)值模擬結(jié)果與試驗結(jié)果一致,在擺動速度接近0 時同樣表現(xiàn)出頻率增大的現(xiàn)象.HDEM 模型能夠模擬滾動試件臨近靜止時刻擺動頻率變高的現(xiàn)象.

圖15 橡膠圓柱數(shù)值模擬與試驗對比結(jié)果Fig.15 Comparison of the numerical simulation and experimental results for the rubber cylinder

圖16 鋁圓柱數(shù)值模擬與試驗對比結(jié)果Fig.16 Comparison of the numerical simulation and experimental results for the aluminum cylinder

從模型計算結(jié)果與試驗結(jié)果對比可以看出,本文建立的滯彈簧滾動阻力模型能夠很好地模擬顆粒滾動速度接近于零狀態(tài)時的能量耗散過程,能夠?qū)︻w粒滾動阻力現(xiàn)象進行合理的滾動阻力機理解釋.建立的滯彈性滾動阻力可為顆粒材料堆積穩(wěn)定問題研究提供方法.

4 結(jié)論

本文研究建立了滾動阻力滯彈性表達的HDEM模型,與常規(guī)DEM 模型的數(shù)值模擬結(jié)果進行了比較.通過圓柱試件在平臺上的自由滾動試驗驗證了該模型的有效性,并得出以下結(jié)論:

(1) HDEM 滾動阻力模型模擬結(jié)果與試驗現(xiàn)象吻合,能較好地解釋顆粒在臨近靜止階段的能量耗散特性;

(2) 彈性滯后引起的滾動阻力包括速度相關型的和速度無關型的兩部分,滾動試件臨近靜止時刻,與速度無關的阻力成分占比越來越大;

(3) HDEM 滾動阻力模型能較好地擬合橡膠材料與鋁材料圓柱試件的擺動頻率,且能很好地模擬試驗中臨近靜止時刻頻率變高的現(xiàn)象.