杨 摄， 郑明军， 吴文江， 赵晨磊

(1 石家庄铁道大学 机械工程学院，河北 石家庄 050043;2 石家庄铁道大学 教务处，河北 石家庄 050043)

0 引言

在研究风沙地区农业机械以及除沙机械等相关设备前，沙土的力学特性及沙土与机械之间相互作用参数的准确性将直接影响机械优化及仿真结果，不同的土壤类型有不同的力学特性。由于风沙地区沙粒物理特性参数的缺失，因而无法建立针对风沙地区沙粒的土壤模型。沙粒同机械及沙粒同橡胶传送带仿真的过程中，参数尚无法确定，常出现与试验结果较大的出入。应用试验及离散元方法研究风沙地区沙粒的物理特性，可以为风沙地区农业机械及除沙机械设计研究奠定基础。目前国内外针对风沙地区沙粒物理参数特性的研究较少。张锐等[1]认为沙土颗粒的形状将会对物理特性参数有很大影响，王宪良等[2]利用砂壤土建立了相关土壤模型。

目前颗粒的物理特性参数，主要是通过试验直接获取以及虚拟参数标定，单一的试验没有考虑到几种参数相互作用时对试验结果的影响[3-9]。沙粒直径较小，不易直接获得颗粒间静摩擦系数、滚动摩擦系数以及恢复系数，目前国内学者主要采用的方法为测量实际堆积角值，采用离散元仿真以及试验数据处理，建立拟合曲线，反推得到[10-18]。

在现有的相关研究基础上，提出一个系统标定沙粒参数的方法。通过大量试验对风沙地区沙粒物理特性参数进行标定，如堆积密度、沙粒直径、质量分数、含水率、空隙率以及堆积角。其次是通过试验确定沙粒与钢板和橡胶的静摩擦系数、滚动摩擦系数以及恢复系数进行了标定。最后是在保证参数正确的情况下，进行相应的微调。运用EDEM软件进行仿真和试验的对比，进行验证，建立针对风沙地区沙粒的土壤模型，为机械除沙以及相关后续研究提供合理可靠的试验及仿真参数。

1 沙粒间参数标定方法与原理

土壤及种粒建模通常采用的方法为离散元法。离散元法是分析和求解复杂离散系统动力学问题的一种新型数值方法，主要通过建立固体颗粒系统的参数化模型，进行颗粒行为模拟和分析[19]。基本思想是将不连续介质离散为刚性颗粒单元，颗粒由于接触，重力位置和速度发生改变，该过程遵循牛顿第二定律[20-22]。目前离散元方法已经广泛应用于矿山、能源化工、散粒体运输和农业等许多领域，郑明军等[23]运用离散元方法分析了沙粒在铁路除沙车内集沙、排沙、抛沙的动态过程和数值规律并进行了相关优化，严战友等[24]通过离散元方法研究了对隧道盾构机施加不同的激励条件下，刀盘推进力和扭矩的变化规律。

图1 离散元接触力学模型

由于风沙地区沙粒的粘结力较小，且沙粒本身具有离散性，又由于在EDEM中Hertz-Mindlin无滑动接触模型为默认计算模型，相比较于其他模型，具有计算精准、速度快等特点，故采用Hertz-Mindlin无滑动接触模型。在此模型中不考虑表面粘结的情况下，法向力采用Hertz理论进行计算，切向力采用Middlin-Deresiewicz 理论进行计算，Hertz理论是颗粒曲面弹性接触问题的理论基础，广泛应用于球体、柱体等曲面体的弹性接触。Hertz与Middlin-Deresiewicz理论均针对颗粒的弹性接触，但是由于变形量和接触力之间是非线性的，法向力与切向力相互耦合很难求解，所以通常采用微量叠加法分别计算法向和切向的位移与接触力[19]，其接触力学模型如图1所示。其中，kn、kt、Cn、Ct、Fn、Ft分别为弹簧法向刚度系数、弹簧切向刚度系数、法向阻尼系数、切向阻尼系数、颗粒间法向力以及颗粒间切向力。

颗粒间的法向力为

(1)

接触时的切向力为

(2)

(3)

2 沙粒本征参数测定

2.1 风沙地区沙粒的含水率

风沙明显区别于河沙及海沙，河沙成分较为复杂，表面具有一定的光滑性，含土量较高。海沙含盐量较高，对金属材料会产生电化学腐蚀。风沙则是由强风将沙漠中粒径较小的沙吹起形成的，其成分主要为二氧化硅，且纯度较高，杂质较少。取风沙试样100 g，取沙地点位于内蒙古乌日根塔拉镇，东经112.767°、北纬42.23°。加热6～8 h至衡重，通过测量加热前后的质量变化，经过多次实验得到加热后沙土质量降低(0.04±0.002) g，其含水率约为0.04%。而陈忠达等[25]，李振纲[26]通过试验认为当沙粒含水率低于1%时，沙粒的干密度较大，此时沙粒的结构十分松散，几乎没有粘结力，仅有内摩擦力。

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

由于沙粒极难溶于水，所以测量沙粒的堆积密度采用排水法。经过30次试验测量，风沙地区的沙粒堆积密度为1 613 kg/m3,其空隙率为38.12%。

2.3 沙粒直径分布与质量分数

由于风沙地区沙粒直径较小，为了给离散元软件颗粒工厂提供准确的颗粒直径，取沙土样本500 g，用GZS-1型高频土壤筛机，将不同直径的沙粒通过标准筛进行颗粒分级，量取不同级别的沙粒质量，从而获得不同直径的沙粒所占的质量分数。试验重复10次，筛出5个级别沙粒试样，如图2所示，试验结果的平均值如表1所示。

图2 筛后沙粒试样

表1 土壤粒径分级及质量分数

查阅《公路土工试验规程》[27]可知, 沙粒范围在2.0～0.5 mm为粗沙;0.5～0.25 mm为中沙;0.25～0.074 mm为细沙;0.074～0.002 mm为粉粒，从试验结果可知试验用风沙主要由中沙、细沙和粉粒组成。

3 沙粒与钢板及橡胶参数测定

3.1 沙粒与钢板及橡胶的恢复系数测定

恢复系数是反映物体碰撞时，物体变形恢复能力的参数，只与材料相关，广泛用于农业领域及化工领域。测定恢复系数主要使用的方法为，碰撞前后的两物体在接触点处的法向相对分离速度与法向相对接近速度之比[28]。根据上述定义，恢复系数的表达式为

(4)

式中,vn′为物料碰撞时法向相对分离速度；vn为物料碰撞时法向相对接近速度。

图3 恢复系数试验原理图

观察颗粒的恢复系数，目前大部分学者采用的方法是将碰撞前后的速度之比转化为高度之比，利用高速摄像机进行捕捉，但是由于沙粒直径过小，因而不易观察到前后高度变化。本文采用如下方法测定风沙地区沙粒与钢板的恢复系数。原理如图3所示。

图3中，H为沙粒自由落体到钢板的高度，H1为沙粒与钢板碰撞后的下落速度，令H=H1，L为沙粒与钢板碰撞点与沙粒在粘板落点之间的水平距离。钢板与水平面夹角为45°。

设沙粒以初速度为0从距离水平面高度为2H处自由下落，经过时间t与钢板接触，沙粒碰撞前速度为v，碰撞后速度的水平分量为v′,碰撞后沙粒下落至粘板的时间为t0。根据运动学原理

(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。通过制作试验平台，经过30次测量取平均值，得出沙粒与钢板的恢复系数为0.478±0.021，沙粒与橡胶的恢复系数为0.175±0.018。

3.2 沙粒与钢板及橡胶的静摩擦系数测定

斜面滑动法是用来测量静摩擦系数常用的方法，在研究土壤及种粒与钢板之间静摩擦系数中应用较为广泛。设质量为m的单个沙粒放在斜面仪的斜面上，α为斜面与水平面的夹角，静摩擦系数μ与夹角α的表达式如下

(11)

图4 静摩擦系数测量装置

试验装置如图4所示，取钢板(Q235)放置于斜面仪上，将沙粒放在钢板之上，逐步增加夹角α，直至沙粒开始下落，当沙粒刚开始下滑时，记录夹角度数，取其正切值，为了便于观察，单个沙粒的粒径选取范围取0.5～1 mm。经过30次测量取平均值，沙粒与钢板的静摩擦系数为0.59±0.01，沙粒同橡胶的静摩擦系数为0.613±0.008。

3.3 沙粒与钢板及橡胶的滚动摩擦系数测定

沙粒与其他材料滚动摩擦系数不易测定，主要采用仿真方法确定。李贝等[9]通过试验方法确定了铁球的滚动摩擦系数，并且通过试验验证，确定了方法的可行性。

图5 滚动摩擦系数测定原理图

将沙粒放置在长为L、倾角为θ、高度为h1的斜面上，让其无初速度滚下至脱离斜面时的速度为v0(水平速度为v1，竖直速度为v2)；自由落体至粘板，落点与斜面终点的水平距离为H，竖直距离为h3。设滚动摩擦系数为k，小球的质量为m，小球自由落体的时间为t。原理示意图如图5所示。

由能量守恒定律可知

(12)

V1=V0cosθ

(13)

V2=V0sinθ

(14)

v1t=H

(15)

(16)

(17)

图6 参数优化流程图

将Q235钢板放在斜面仪上，沙粒放在钢板上，逐步增加夹角θ，直至沙粒全部滑落，记录夹角度数，取其正切值，为了便于观察，取单个沙粒的粒径选取范围是0.5～1 mm。经过30次测量取平均值，沙粒与钢板之间滚动摩擦系数为0.28±0.009，沙粒同橡胶之间滚动摩擦系数为0.45±0.011。

4 沙粒间物理参数虚拟标定

由于沙粒直径较小，沙粒间静摩擦系数、滚动摩擦系数以及恢复系数无法通过试验方法直接测量出来，故采用仿真与试验相结合的方法，通过数值分析及优化，对以上3个参数进行测定。参数优化流程如图6所示。

通过以上试验及相关参考文献确定离散元模型参数，如表2所示。

表2 离散元仿真模型参数

由表1知，颗粒直径在0.1～0.5 mm的比例占总数的97.58%，颗粒建模时以直径为0.1～0.5 mm为基础。如果按照实际尺寸及比例生成颗粒，仿真时间将会延长，效率降低。为了加快仿真速度，提高效率，将沙粒的直径放大0.8～2倍，并利用颗粒工厂对其进行随机生成。沙粒的形状主要为球形、长条形、菱形[1,32]。为了能够更加真实准确模拟堆积角，选取标准球型，对以上3种形状进行填充。

在多次预试验过程中发现，在下落高度一致的情况下，堆积角能否形成取决于颗粒数量的大小，颗粒数量过多将会使仿真时间延长，仿真效率降低。颗粒数量较少，不能形成堆积角，经过大量预试验，确定颗粒工厂共生成颗粒30 000颗，生成速率为15 000颗/s，颗粒生成后，从漏斗中下落。为了确保颗粒能够达到稳定状态，将仿真时间设置为4 s，然后对其堆积角进行测量。

颗粒间相应参数，通常采用的方法为先测量堆积角，通过试验数据及分析反推得到相应的参数。经过大量预试验及查阅文献[30]和文献[33]，确定沙粒间静摩擦系数所在区间为0.2～0.28，沙粒间滚动摩擦系数所在区间为0.05～0.2，沙粒间恢复系数所在区间为0.15～0.35。

因待确定因素较多，故本次试验采用回归正交试验，建立拟合方程，对自变量参数进行求解。仿真试验因素取值范围、设计及结果如表3、表4所示。

表3 仿真试验因素取值范围

表4 仿真试验方案设计与结果

对试验数据进行分析，得到沙粒堆积角与试验变量的二阶回归模型，其回归方程为

θ=-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)

对仿真结果进行方差分析，拟合模型中的P值小于0.000 1，表明利用此模型表示自变量与因变量的关系极为显著，其中沙粒-沙粒静摩擦系数二次项P值均小于0.000 1，相比较于其他因素，沙粒静摩擦系数对于堆积角的影响最为显著。试验模型整体极为显著，表明对指标有显著影响的因素已经考虑到，说明试验合理有效。决定系数R2=0.937 3,校正决定系数Radj2=0.893 9,两者均接近于1，表明方程可靠度较高。

5 试验验证

研究采用沙粒堆积角测试试验装置如图7所示，将500 g沙粒倒入漏斗中，沙粒在漏斗下方形成颗粒堆。为了减少人为测量误差，更加准确地对堆积角进行测量，从4个方向对堆积角进行拍照，使用Matlab对图像进行灰度化、二值化处理、孔洞填充、提取边界，利用最小二乘法对边界曲线进行直线拟合，拟合直线的斜率即为堆积角的正切值。仿真及试验结果如图8及表5所示。

图7 漏斗法试验装置

图8 堆积角图像处理

表5 堆积角试验与仿真数据对比

以实际堆积角27.94°为目标值，对回归方程进行有约束目标的优化求解。将优化结果通过离散元软件进行仿真确认并进行进一步优化确认，得到沙粒间静摩擦系数为0.23、滚动摩擦系数为0.1、沙粒间恢复系数为0.25。由图8及表5中可以看出，仿真的堆积曲线与试验堆积曲线拟合较好，沙粒间静摩擦系数、滚动摩擦系数与恢复系数接近于真实值，分别为0.23、0.1、0.25。

6 结论

(1)以风沙地区沙粒作为研究对象，通过高频土壤筛机、堆积角测量装置及相关试验设备，确定了风沙地区沙粒的颗粒直径、堆积角、堆积密度、含水率，提出一种系统标定沙粒与其他材料相互作用参数的方法，并首次测量得到沙粒与钢板(Q235)及橡胶的静摩擦系数、滚动摩擦系数以及恢复系数。

(2)利用仿真与试验相结合的方法，通过回归正交试验方法，对离散元仿真数据进行分析，建立回归方程，采用漏斗法进行堆积角试验，利用Matlab图像技术测量堆积角，以实际堆积角为目标值，对方程进行求解，依次求得沙粒间静摩擦系数为0.23，沙粒间滚动摩擦系数为0.1，沙粒间恢复系数为0.25。通过试验对优化参数进行验证，对比2条曲线，拟合度较高，堆积角差异值为0.39°，相对误差值为1.39%，误差值较小。说明针对风沙地区沙粒物理特性参数的标定方法是可行的，测量数据真实有效，从而为风沙地区农用机械及除沙机械仿真、设计与优化提供可靠的试验数据。