基于離散元的風(fēng)沙沙粒物理特性參數(shù)標(biāo)定

楊 攝， 鄭明軍， 吳文江， 趙晨磊

(1 石家莊鐵道大學(xué) 機械工程學(xué)院，河北 石家莊 050043;2 石家莊鐵道大學(xué) 教務(wù)處，河北 石家莊 050043)

0 引言

在研究風(fēng)沙地區(qū)農(nóng)業(yè)機械以及除沙機械等相關(guān)設(shè)備前，沙土的力學(xué)特性及沙土與機械之間相互作用參數(shù)的準(zhǔn)確性將直接影響機械優(yōu)化及仿真結(jié)果，不同的土壤類型有不同的力學(xué)特性。由于風(fēng)沙地區(qū)沙粒物理特性參數(shù)的缺失，因而無法建立針對風(fēng)沙地區(qū)沙粒的土壤模型。沙粒同機械及沙粒同橡膠傳送帶仿真的過程中，參數(shù)尚無法確定，常出現(xiàn)與試驗結(jié)果較大的出入。應(yīng)用試驗及離散元方法研究風(fēng)沙地區(qū)沙粒的物理特性，可以為風(fēng)沙地區(qū)農(nóng)業(yè)機械及除沙機械設(shè)計研究奠定基礎(chǔ)。目前國內(nèi)外針對風(fēng)沙地區(qū)沙粒物理參數(shù)特性的研究較少。張銳等[1]認為沙土顆粒的形狀將會對物理特性參數(shù)有很大影響，王憲良等[2]利用砂壤土建立了相關(guān)土壤模型。

目前顆粒的物理特性參數(shù)，主要是通過試驗直接獲取以及虛擬參數(shù)標(biāo)定，單一的試驗沒有考慮到幾種參數(shù)相互作用時對試驗結(jié)果的影響[3-9]。沙粒直徑較小，不易直接獲得顆粒間靜摩擦系數(shù)、滾動摩擦系數(shù)以及恢復(fù)系數(shù)，目前國內(nèi)學(xué)者主要采用的方法為測量實際堆積角值，采用離散元仿真以及試驗數(shù)據(jù)處理，建立擬合曲線，反推得到[10-18]。

在現(xiàn)有的相關(guān)研究基礎(chǔ)上，提出一個系統(tǒng)標(biāo)定沙粒參數(shù)的方法。通過大量試驗對風(fēng)沙地區(qū)沙粒物理特性參數(shù)進行標(biāo)定，如堆積密度、沙粒直徑、質(zhì)量分?jǐn)?shù)、含水率、空隙率以及堆積角。其次是通過試驗確定沙粒與鋼板和橡膠的靜摩擦系數(shù)、滾動摩擦系數(shù)以及恢復(fù)系數(shù)進行了標(biāo)定。最后是在保證參數(shù)正確的情況下，進行相應(yīng)的微調(diào)。運用EDEM軟件進行仿真和試驗的對比，進行驗證，建立針對風(fēng)沙地區(qū)沙粒的土壤模型，為機械除沙以及相關(guān)后續(xù)研究提供合理可靠的試驗及仿真參數(shù)。

1 沙粒間參數(shù)標(biāo)定方法與原理

土壤及種粒建模通常采用的方法為離散元法。離散元法是分析和求解復(fù)雜離散系統(tǒng)動力學(xué)問題的一種新型數(shù)值方法，主要通過建立固體顆粒系統(tǒng)的參數(shù)化模型，進行顆粒行為模擬和分析[19]?；舅枷胧菍⒉贿B續(xù)介質(zhì)離散為剛性顆粒單元，顆粒由于接觸，重力位置和速度發(fā)生改變，該過程遵循牛頓第二定律[20-22]。目前離散元方法已經(jīng)廣泛應(yīng)用于礦山、能源化工、散粒體運輸和農(nóng)業(yè)等許多領(lǐng)域，鄭明軍等[23]運用離散元方法分析了沙粒在鐵路除沙車內(nèi)集沙、排沙、拋沙的動態(tài)過程和數(shù)值規(guī)律并進行了相關(guān)優(yōu)化，嚴(yán)戰(zhàn)友等[24]通過離散元方法研究了對隧道盾構(gòu)機施加不同的激勵條件下，刀盤推進力和扭矩的變化規(guī)律。

圖1 離散元接觸力學(xué)模型

由于風(fēng)沙地區(qū)沙粒的粘結(jié)力較小，且沙粒本身具有離散性，又由于在EDEM中Hertz-Mindlin無滑動接觸模型為默認計算模型，相比較于其他模型，具有計算精準(zhǔn)、速度快等特點，故采用Hertz-Mindlin無滑動接觸模型。在此模型中不考慮表面粘結(jié)的情況下，法向力采用Hertz理論進行計算，切向力采用Middlin-Deresiewicz 理論進行計算，Hertz理論是顆粒曲面彈性接觸問題的理論基礎(chǔ)，廣泛應(yīng)用于球體、柱體等曲面體的彈性接觸。Hertz與Middlin-Deresiewicz理論均針對顆粒的彈性接觸，但是由于變形量和接觸力之間是非線性的，法向力與切向力相互耦合很難求解，所以通常采用微量疊加法分別計算法向和切向的位移與接觸力[19]，其接觸力學(xué)模型如圖1所示。其中，kn、kt、Cn、Ct、Fn、Ft分別為彈簧法向剛度系數(shù)、彈簧切向剛度系數(shù)、法向阻尼系數(shù)、切向阻尼系數(shù)、顆粒間法向力以及顆粒間切向力。

顆粒間的法向力為

(1)

接觸時的切向力為

(2)

(3)

2 沙粒本征參數(shù)測定

2.1 風(fēng)沙地區(qū)沙粒的含水率

風(fēng)沙明顯區(qū)別于河沙及海沙，河沙成分較為復(fù)雜，表面具有一定的光滑性，含土量較高。海沙含鹽量較高，對金屬材料會產(chǎn)生電化學(xué)腐蝕。風(fēng)沙則是由強風(fēng)將沙漠中粒徑較小的沙吹起形成的，其成分主要為二氧化硅，且純度較高，雜質(zhì)較少。取風(fēng)沙試樣100 g，取沙地點位于內(nèi)蒙古烏日根塔拉鎮(zhèn)，東經(jīng)112.767°、北緯42.23°。加熱6～8 h至衡重，通過測量加熱前后的質(zhì)量變化，經(jīng)過多次實驗得到加熱后沙土質(zhì)量降低(0.04±0.002) g，其含水率約為0.04%。而陳忠達等[25]，李振綱[26]通過試驗認為當(dāng)沙粒含水率低于1%時，沙粒的干密度較大，此時沙粒的結(jié)構(gòu)十分松散，幾乎沒有粘結(jié)力，僅有內(nèi)摩擦力。

2.2 沙粒的堆積密度及空隙率

由于沙粒極難溶于水，所以測量沙粒的堆積密度采用排水法。經(jīng)過30次試驗測量，風(fēng)沙地區(qū)的沙粒堆積密度為1 613 kg/m3,其空隙率為38.12%。

2.3 沙粒直徑分布與質(zhì)量分?jǐn)?shù)

由于風(fēng)沙地區(qū)沙粒直徑較小，為了給離散元軟件顆粒工廠提供準(zhǔn)確的顆粒直徑，取沙土樣本500 g，用GZS-1型高頻土壤篩機，將不同直徑的沙粒通過標(biāo)準(zhǔn)篩進行顆粒分級，量取不同級別的沙粒質(zhì)量，從而獲得不同直徑的沙粒所占的質(zhì)量分?jǐn)?shù)。試驗重復(fù)10次，篩出5個級別沙粒試樣，如圖2所示，試驗結(jié)果的平均值如表1所示。

圖2 篩后沙粒試樣

表1 土壤粒徑分級及質(zhì)量分?jǐn)?shù)

查閱《公路土工試驗規(guī)程》[27]可知, 沙粒范圍在2.0～0.5 mm為粗沙;0.5～0.25 mm為中沙;0.25～0.074 mm為細沙;0.074～0.002 mm為粉粒，從試驗結(jié)果可知試驗用風(fēng)沙主要由中沙、細沙和粉粒組成。

3 沙粒與鋼板及橡膠參數(shù)測定

3.1 沙粒與鋼板及橡膠的恢復(fù)系數(shù)測定

恢復(fù)系數(shù)是反映物體碰撞時，物體變形恢復(fù)能力的參數(shù)，只與材料相關(guān)，廣泛用于農(nóng)業(yè)領(lǐng)域及化工領(lǐng)域。測定恢復(fù)系數(shù)主要使用的方法為，碰撞前后的兩物體在接觸點處的法向相對分離速度與法向相對接近速度之比[28]。根據(jù)上述定義，恢復(fù)系數(shù)的表達式為

(4)

式中,vn′為物料碰撞時法向相對分離速度；vn為物料碰撞時法向相對接近速度。

圖3 恢復(fù)系數(shù)試驗原理圖

觀察顆粒的恢復(fù)系數(shù)，目前大部分學(xué)者采用的方法是將碰撞前后的速度之比轉(zhuǎn)化為高度之比，利用高速攝像機進行捕捉，但是由于沙粒直徑過小，因而不易觀察到前后高度變化。本文采用如下方法測定風(fēng)沙地區(qū)沙粒與鋼板的恢復(fù)系數(shù)。原理如圖3所示。

圖3中，H為沙粒自由落體到鋼板的高度，H1為沙粒與鋼板碰撞后的下落速度，令H=H1，L為沙粒與鋼板碰撞點與沙粒在粘板落點之間的水平距離。鋼板與水平面夾角為45°。

設(shè)沙粒以初速度為0從距離水平面高度為2H處自由下落，經(jīng)過時間t與鋼板接觸，沙粒碰撞前速度為v，碰撞后速度的水平分量為v′,碰撞后沙粒下落至粘板的時間為t0。根據(jù)運動學(xué)原理

(5)

不計空氣阻力及摩擦，沙粒碰撞后其垂直速度為0，且H=H1，可得

(6)

(7)

vn′=v′sin45°

(8)

vn=vsin45°

(9)

e=(vn′)/vn=(v′sin45°)/(vsin45°)=L/2H

(10)

為了便于觀察，單個沙粒的粒徑選取范圍是0.5～1 mm。通過制作試驗平臺，經(jīng)過30次測量取平均值，得出沙粒與鋼板的恢復(fù)系數(shù)為0.478±0.021，沙粒與橡膠的恢復(fù)系數(shù)為0.175±0.018。

3.2 沙粒與鋼板及橡膠的靜摩擦系數(shù)測定

斜面滑動法是用來測量靜摩擦系數(shù)常用的方法，在研究土壤及種粒與鋼板之間靜摩擦系數(shù)中應(yīng)用較為廣泛。設(shè)質(zhì)量為m的單個沙粒放在斜面儀的斜面上，α為斜面與水平面的夾角，靜摩擦系數(shù)μ與夾角α的表達式如下

(11)

圖4 靜摩擦系數(shù)測量裝置

試驗裝置如圖4所示，取鋼板(Q235)放置于斜面儀上，將沙粒放在鋼板之上，逐步增加夾角α，直至沙粒開始下落，當(dāng)沙粒剛開始下滑時，記錄夾角度數(shù)，取其正切值，為了便于觀察，單個沙粒的粒徑選取范圍取0.5～1 mm。經(jīng)過30次測量取平均值，沙粒與鋼板的靜摩擦系數(shù)為0.59±0.01，沙粒同橡膠的靜摩擦系數(shù)為0.613±0.008。

3.3 沙粒與鋼板及橡膠的滾動摩擦系數(shù)測定

沙粒與其他材料滾動摩擦系數(shù)不易測定，主要采用仿真方法確定。李貝等[9]通過試驗方法確定了鐵球的滾動摩擦系數(shù)，并且通過試驗驗證，確定了方法的可行性。

圖5 滾動摩擦系數(shù)測定原理圖

將沙粒放置在長為L、傾角為θ、高度為h1的斜面上，讓其無初速度滾下至脫離斜面時的速度為v0(水平速度為v1，豎直速度為v2)；自由落體至粘板，落點與斜面終點的水平距離為H，豎直距離為h3。設(shè)滾動摩擦系數(shù)為k，小球的質(zhì)量為m，小球自由落體的時間為t。原理示意圖如圖5所示。

由能量守恒定律可知

(12)

V1=V0cosθ

(13)

V2=V0sinθ

(14)

v1t=H

(15)

(16)

(17)

圖6 參數(shù)優(yōu)化流程圖

將Q235鋼板放在斜面儀上，沙粒放在鋼板上，逐步增加夾角θ，直至沙粒全部滑落，記錄夾角度數(shù)，取其正切值，為了便于觀察，取單個沙粒的粒徑選取范圍是0.5～1 mm。經(jīng)過30次測量取平均值，沙粒與鋼板之間滾動摩擦系數(shù)為0.28±0.009，沙粒同橡膠之間滾動摩擦系數(shù)為0.45±0.011。

4 沙粒間物理參數(shù)虛擬標(biāo)定

由于沙粒直徑較小，沙粒間靜摩擦系數(shù)、滾動摩擦系數(shù)以及恢復(fù)系數(shù)無法通過試驗方法直接測量出來，故采用仿真與試驗相結(jié)合的方法，通過數(shù)值分析及優(yōu)化，對以上3個參數(shù)進行測定。參數(shù)優(yōu)化流程如圖6所示。

通過以上試驗及相關(guān)參考文獻確定離散元模型參數(shù)，如表2所示。

表2 離散元仿真模型參數(shù)

由表1知，顆粒直徑在0.1～0.5 mm的比例占總數(shù)的97.58%，顆粒建模時以直徑為0.1～0.5 mm為基礎(chǔ)。如果按照實際尺寸及比例生成顆粒，仿真時間將會延長，效率降低。為了加快仿真速度，提高效率，將沙粒的直徑放大0.8～2倍，并利用顆粒工廠對其進行隨機生成。沙粒的形狀主要為球形、長條形、菱形[1,32]。為了能夠更加真實準(zhǔn)確模擬堆積角，選取標(biāo)準(zhǔn)球型，對以上3種形狀進行填充。

在多次預(yù)試驗過程中發(fā)現(xiàn)，在下落高度一致的情況下，堆積角能否形成取決于顆粒數(shù)量的大小，顆粒數(shù)量過多將會使仿真時間延長，仿真效率降低。顆粒數(shù)量較少，不能形成堆積角，經(jīng)過大量預(yù)試驗，確定顆粒工廠共生成顆粒30 000顆，生成速率為15 000顆/s，顆粒生成后，從漏斗中下落。為了確保顆粒能夠達到穩(wěn)定狀態(tài)，將仿真時間設(shè)置為4 s，然后對其堆積角進行測量。

顆粒間相應(yīng)參數(shù)，通常采用的方法為先測量堆積角，通過試驗數(shù)據(jù)及分析反推得到相應(yīng)的參數(shù)。經(jīng)過大量預(yù)試驗及查閱文獻[30]和文獻[33]，確定沙粒間靜摩擦系數(shù)所在區(qū)間為0.2～0.28，沙粒間滾動摩擦系數(shù)所在區(qū)間為0.05～0.2，沙粒間恢復(fù)系數(shù)所在區(qū)間為0.15～0.35。

因待確定因素較多，故本次試驗采用回歸正交試驗，建立擬合方程，對自變量參數(shù)進行求解。仿真試驗因素取值范圍、設(shè)計及結(jié)果如表3、表4所示。

表3 仿真試驗因素取值范圍

表4 仿真試驗方案設(shè)計與結(jié)果

對試驗數(shù)據(jù)進行分析，得到沙粒堆積角與試驗變量的二階回歸模型，其回歸方程為

θ=-579.78+4 260.17A+262.05B+340.04C+1 032.43AB-

1 037.93AC+673.83BC-7 866.39A2-2 385.69B2-380.40C2

(18)

對仿真結(jié)果進行方差分析，擬合模型中的P值小于0.000 1，表明利用此模型表示自變量與因變量的關(guān)系極為顯著，其中沙粒-沙粒靜摩擦系數(shù)二次項P值均小于0.000 1，相比較于其他因素，沙粒靜摩擦系數(shù)對于堆積角的影響最為顯著。試驗?zāi)Ｐ驼w極為顯著，表明對指標(biāo)有顯著影響的因素已經(jīng)考慮到，說明試驗合理有效。決定系數(shù)R2=0.937 3,校正決定系數(shù)Radj2=0.893 9,兩者均接近于1，表明方程可靠度較高。

5 試驗驗證

研究采用沙粒堆積角測試試驗裝置如圖7所示，將500 g沙粒倒入漏斗中，沙粒在漏斗下方形成顆粒堆。為了減少人為測量誤差，更加準(zhǔn)確地對堆積角進行測量，從4個方向?qū)Χ逊e角進行拍照，使用Matlab對圖像進行灰度化、二值化處理、孔洞填充、提取邊界，利用最小二乘法對邊界曲線進行直線擬合，擬合直線的斜率即為堆積角的正切值。仿真及試驗結(jié)果如圖8及表5所示。

圖7 漏斗法試驗裝置

圖8 堆積角圖像處理

表5 堆積角試驗與仿真數(shù)據(jù)對比

以實際堆積角27.94°為目標(biāo)值，對回歸方程進行有約束目標(biāo)的優(yōu)化求解。將優(yōu)化結(jié)果通過離散元軟件進行仿真確認并進行進一步優(yōu)化確認，得到沙粒間靜摩擦系數(shù)為0.23、滾動摩擦系數(shù)為0.1、沙粒間恢復(fù)系數(shù)為0.25。由圖8及表5中可以看出，仿真的堆積曲線與試驗堆積曲線擬合較好，沙粒間靜摩擦系數(shù)、滾動摩擦系數(shù)與恢復(fù)系數(shù)接近于真實值，分別為0.23、0.1、0.25。

6 結(jié)論

(1)以風(fēng)沙地區(qū)沙粒作為研究對象，通過高頻土壤篩機、堆積角測量裝置及相關(guān)試驗設(shè)備，確定了風(fēng)沙地區(qū)沙粒的顆粒直徑、堆積角、堆積密度、含水率，提出一種系統(tǒng)標(biāo)定沙粒與其他材料相互作用參數(shù)的方法，并首次測量得到沙粒與鋼板(Q235)及橡膠的靜摩擦系數(shù)、滾動摩擦系數(shù)以及恢復(fù)系數(shù)。

(2)利用仿真與試驗相結(jié)合的方法，通過回歸正交試驗方法，對離散元仿真數(shù)據(jù)進行分析，建立回歸方程，采用漏斗法進行堆積角試驗，利用Matlab圖像技術(shù)測量堆積角，以實際堆積角為目標(biāo)值，對方程進行求解，依次求得沙粒間靜摩擦系數(shù)為0.23，沙粒間滾動摩擦系數(shù)為0.1，沙粒間恢復(fù)系數(shù)為0.25。通過試驗對優(yōu)化參數(shù)進行驗證，對比2條曲線，擬合度較高，堆積角差異值為0.39°，相對誤差值為1.39%，誤差值較小。說明針對風(fēng)沙地區(qū)沙粒物理特性參數(shù)的標(biāo)定方法是可行的，測量數(shù)據(jù)真實有效，從而為風(fēng)沙地區(qū)農(nóng)用機械及除沙機械仿真、設(shè)計與優(yōu)化提供可靠的試驗數(shù)據(jù)。