基于EDEM的粮食立筒仓卸料动态侧压力研究

刘强, 张茜茜, 刘涛, 丁永刚, 许启铿, 孙启帅

(1.河南工业大学土木工程学院, 郑州 450001; 2.河南省粮油仓储建筑与安全重点实验室, 郑州 450001;3.机械工业第六设计研究院有限公司, 郑州 450001)

立筒仓以其占地面积小、易于机械化、自动化作业、流通费用低、造价低等特点,被广泛应用于粮食储备。但贮料[1]与仓体存在着复杂的相互作用,贮料会受到仓壁的约束作用,在卸料过程中又会对仓壁产生一定的反作用力,发生超压现象从而加大仓体破坏的风险,影响粮食的供给保障。

随着筒仓需求量越来越大,相关的研究也在不断完善。Wang等[2]运用细砂和玻璃微珠两种散体材料,对比研究筒仓卸料时试验和模拟的流动规律及法向壁压分布情况,得出试验和模拟的结果基本一致。Wójcik等[3]通过多次装卸料试验,得出在整体流状态时,最大水平压力和竖向压力位于仓壁底部范围内。Fullard[4]采用PIV对筒仓的重力驱动卸料试验进行录像分析,表明了在一个比例时间内,远离排放孔,流动行为是相同的。原方等[5]对比分析普通筒仓安装减压管时筒仓单侧壁卸料的动态侧压力,得出减压管能够减小卸料时储料对仓壁的冲击,从而提高卸料时筒仓的稳定性。王振清等[6]进行了50%、80%和100%仓容3 种状态下的侧压力试验,发现管状流动的出现位置与初始的储粮仓容相关,结果又与规范值进行对比,可为粮食筒仓设计提供参考。徐志军等[7]分析颗粒在普通筒仓双侧壁卸料和带流槽侧壁卸料过程中的力学行为,验证了溜槽的减压机理。程远浩[8]将曲线漏斗筒仓与锥形漏斗筒仓对比,曲线漏斗改善流态,从而降低了筒仓动态侧压力。石鑫等[9]建立的改进颗粒模型与传统单元模型比较,发现其仓壁侧压力峰值显著减小,力链分叉现象明显。余汉华等[10]采用3种机器学习方法预测筒仓动态侧压力,同时进行随机抽样得到筒仓动态侧压力的概率分布,为筒仓结构的可靠度提供了理论基础。李坤由等[11]通过EDEM(engineering discrete element method)离散元软件,对比分析不同摩擦情况下力链的细观参数,得出颗粒间摩擦系数影响卸料时间,增大颗粒间滚动摩擦会增加拱效应。

综上可知,前人研究中针对卸料过程中形成的不同流动状态没有具体分析其成因;在静态、动态侧压力等方面的研究成果虽有很多,但是对于超压系数没有固定的取值,且最大值范围也并不确定;同时卸料过程中形成的动态拱具有连续形成和倒塌的特性,会对研究带来较大的困难,有关动态拱的形成机理仍缺少相应解释。为此,应用离散元软件EDEM(engineering discrete element method),对比分析粮食立筒仓卸料时静态、动态侧压力和超压系数差异规律,探究筒仓卸料时成拱机理与仓壁侧压力的关系,为立筒仓动态侧压力提供一种新的研究途径。

1 离散元模型建立

1.1 立筒仓模型

以某粮库储藏大豆的柱承式钢筋混凝土立筒仓为研究对象,在本团队已开展的试验研究[12]基础上,通过EDEM(engineering discrete element method)离散元软件建立粮食立筒仓模型,采用有机玻璃材料属性,模型总高度为1.62 m,其中仓壁高度1.28 m,内径0.48 m,外径0.50 m。结构模型物理参数泊松比为0.33,密度为1 180 kg/m3,弹性模量为3.05×103MPa,结构模型整体与生成贮料区域重力加速度为9.8 m/s2,粮食立筒仓模型如图1所示。

图1 立筒仓模型Fig.1 Silo model

1.2 贮料颗粒模型

所采用的贮料为陶粒,其为球形颗粒形态,在EDEM(Engineering Discrete Element Method)软件前处理模块的Bulk Material中创建贮料颗粒,选取Spheres,设置半径为0.001 5 m,将颗粒的粒径分布设置为正态分布,平均值设为1即表示所有颗粒平均半径大小在长度方向与当前所建颗粒模型大小一致,标准差设为0.05,并通过Calculate Properties自动计算出颗粒质量、体积和转动惯量X、Y、Z。

考虑到所用散体贮料为陶粒,颗粒数量多且在卸料过程中颗粒之间、颗粒与筒仓结构之间接触时间较长,可通过加速度定律求解颗粒间的接触碰撞力;同时颗粒之间、颗粒与筒仓结构之间无黏结,故选择Hertz-Mindlin(no-slip)模型以模拟其接触关系。依据材性试验与相关文献[13-14],颗粒接触属性参数与流动参数如表1所示。满仓装粮线为立筒仓顶部下降0.1 m,即贮料计算高度hn=1.18 m(hn为装粮线至仓身下部的高度)。

表1 仿真参数Table 1 Simulation parameters

1.3 测点选定

根据立筒仓模型结构的对称性和试验条件,仿真模型数据测点的选取与试验[12]所布测点位置保持一致,其中P1～P5用来采集贮料颗粒对仓壁产生的侧压力,分别对应测点深度0.2、0.5、0.8、1.1、1.256 m,且由于卸料方式为中心卸料,同一高度处的侧压力数据是相同的,故选取一列数据进行分析。因筒仓为三维对称结构,故截取筒仓的中部垂直薄片进行数据分析,该垂直截面薄片设置为20 mm厚。在与侧压力测点P1～P5对应高度处,设置颗粒网格组(高度30 mm),其中V1～V5测点以获取与侧压力测点高度一致的不同高度处颗粒速度,如图2所示。

1.4 卸料工况

在进行卸料侧压力分析时,由于不同流速状态下所产生的动态侧压力[15]与超压系数会存在明显不同,故将其分为两种工况,流速慢与流速快工况分别命名为卸料a、卸料b,进而分析两种流速卸料工况下粮食立筒仓静态侧压力、动态侧压力与超压系数。

2 模拟结果与分析

2.1 贮料流动状态与速度分布

根据EDEM数值仿真和试验结果,卸料流速不同时,两种工况下的流动状态性质相同,故进行流动状态与速度分布分析时,选取卸料a(流速慢)工况为例,分别采用直接观察和整体流量指数法进行卸料流动状态分析。

2.1.1 直观观察卸料流态

通过截取不同时刻所选取的筒仓垂直截面薄片,直接观察筒仓卸料流态,如图3(a)～图3(g)所示;同时根据筒仓卸料的实际工作情况,选取贮料顶部的5个颗粒进行标记,隐藏其它颗粒,观察贮料颗粒的运动轨迹,如图3(h)所示。

分图题为卸料时间

由图3(a)可知,当粮食达到稳定状态后,筒仓开始进行卸料;由图3(b)、图3(c)可知,贮料为整体流动状态,贮料上部表现为整体下降,下部呈漏斗流,卸料流动下降较缓慢;由图3(d)、图3(e)可知,当贮料下降至0.8 m测点深度处时,开始出现管状流动趋势,中间部位贮料流速明显更快于仓壁两侧的贮料流速,仓壁附近的贮料颗粒向中心移动,并逐渐变为管状流动阶段,贮料下降明显,流动均匀;由图3(f)、图3(g)可知,在卸料的末期,贮料表现出漏斗状,流动加快,以卸料口正上方范围内的颗粒速度最快先行流出,而靠近漏斗壁的颗粒由于阻碍作用形成死料区导致流出较晚,直至卸料完成。

随着贮料颗粒的流动,在筒仓下部的颗粒之间、颗粒与仓壁之间的摩擦和挤压更显著,导致流动阻力增大,并在筒仓下部范围内达到最大,但是此时筒仓上部的贮料颗粒不受此影响,阻力较小且流动较快,使得筒仓下部贮料更加密实,流动更加困难,相对应使得仓壁侧压力增加。

2.1.2 速度分布与整体流量指数

通过建立的V1～V5颗粒网格组,对整体卸料的贮料颗粒速度分布情况进行统计分析,根据模型内径宽度,每个网格组在水平方向可均分为19个小网格,用于速度数据的采集,所提取速度代表为某一时刻颗粒经过该网格组时的竖向平均速度,如图2(b)所示,在不同高度的竖向平均速度如表2所示。

表2 不同高度处竖向平均速度分布Table 2 Vertical average velocity distribution at different heights

通过整体流量指数[16](global flow index, MFI)对卸料流态情况进行量化分析,分析每个颗粒网格组速度分布,计算公式为

(1)

式(1)中:Vwall为筒仓内壁边界处的贮料颗粒运动速度;Vcenter为筒仓内中心处的贮料颗粒运动速度。

当整体流量指数值大于0.3时该贮料流态为整体流,小于0.3时该贮料流态为漏斗流。因此根据表2计算不同高度处MFI值,其中V1和V2高度处的贮料表现为整体流,而从V3开始,贮料改变流态,逐渐呈漏斗流,结合直观观察卸料流态,并与试验流态进行对比,两者卸料时的流动状态大致相同。

由图4可知,在卸料整体过程中,中间贮料流动速度一直大于仓壁贮料流动速度,且筒仓仓壁[图3(h)]竖向范围内由上至下的不同深度位置的速度值逐渐减小,中间部位[图3(h)]竖向范围内由上至下的不同深度位置的速度值逐渐增大,即MFI值在逐渐减小,说明越靠近漏斗部分,贮料卸料流动状态越接近漏斗流;在筒仓漏斗部分(V5测点处),因其接近于卸料口,故一直表现为漏斗流,筒仓在进行卸料作业时,其贮料所表现出的流态与所处于筒仓的位置有关。

括号中数值为MFI

2.2 立筒仓仓体侧压力分析

2.2.1 静态侧压力

图5为立筒仓满仓静态侧压力模拟值、试验值和规范值对比结果。由图5可知,静态侧压力的模拟值同试验值、规范值[17]一致,均随着测点深度的增加均呈增大趋势,最大值出现在仓壁与漏斗连接部位范围内,静态侧压力的模拟值较大于规范值,规范公式计算值具有一定的局限性且偏于保守,本次卸料仿真的模拟值与试验值接近,为后续的动态侧压力与超压系数分析提供可靠的依据。

图5 立筒仓满仓静态侧压力模拟值、试验值和规范值Fig.5 Simulation, test and specification values of static lateral pressure of silo full

2.2.2 动态侧压力与超压系数

图6为立筒仓两种卸料工况下动态侧压力时程曲线。由图6(a)可知,13 s时在0.8 m测点深度处压力突增明显,动态侧压力达8.79 kPa,是静态侧压力的2.14倍;由图6(b)可知,14 s时在0.8 m测点深度处压力突增明显,动态侧压力达8.77 kPa,是静态侧压力的2.27倍。在卸料a与卸料b模拟时间0～150 s和0～48 s,贮料向下滑移,由于贮料颗粒之间、颗粒与仓壁之间的相互挤压摩擦,一直对仓壁产生冲击,致使仓壁上的P1～P4测点均出现了压力突增的现象,并且发生震荡现象;在漏斗处即测点P5处,由于离卸料口处最近,且在此处流态型式一直为漏斗流,会相应产生一部分死料区,故从卸料0 s开始,动态侧压力明显迅速下降至某一段范围内发生震荡现象直至试验结束,且无超压现象发生。

图6 立筒仓两种卸料工况下动态侧压力时程曲线Fig.6 Dynamic lateral pressure time history curve of silo under two discharge conditions

在卸料过程中呈现出的震荡现象,是由于卸料时贮料流动的过程中在1/3hn范围内形成了超压动力拱,使贮料内部反复的结拱和塌拱,拱脚对仓壁不断进行冲击造成。

由表3、图7可知,EDEM仿真模拟所得的静态、动态侧压力数据值与试验做对比,两者相差较小;筒仓仓壁且除漏斗处的各测点动态侧压力值相较于静态侧压力值均有不同程度的增加,存在超压现象,产生这种现象的原因是由于上部贮料重力和下部贮料的持续流出,导致发生成拱和破拱的现象,在0.8 m测点深度处成拱和破拱速度更快。在漏斗处即1.265 m测点深度发生动态侧压力减小的现象。随着物料的流出,仓内贮料不断地重复“结拱起始—结拱完成—拱塌落”的动态拱过程[18],因而仓壁卸料压力呈现“较小—局部较大—较小”的震荡分布。

表3 立筒仓侧压力模拟值与试验值Table 3 Simulation and test values of silo lateral pressure

图7 立筒仓静态侧压力和动态侧压力Fig.7 Silo static lateral pressure and dynamic lateral pressure

由图8可知,超压系数整体呈现先升后降的趋势,其中在0.8 m测点深度处最大,模拟卸料a超压系数(2.14)和模拟卸料b超压系数(2.27)均超过规范中水平侧压力修正系数Ch,其他测点超压系数均在规范值内;在0.8 m测点深度处模拟卸料a、b相较于试验卸料a、b的超压系数误差在3.4%、7.6%,差别较小且在合理误差范围内;且对比在不同测点深度处两个工况超压系数的大小,可发现不同卸料流速下,流速快工况的超压系数均大于流速慢工况的超压系数。

3 结论

通过对粮食立筒仓卸料进行EDEM(Engineering Discrete Element Method)仿真模拟分析,通过直接观察和整体流量指数分析贮料流态,对比研究了不同流速下模拟值与试验结果的贮料静态压力数据、动态侧压力分布规律与超压系数,得出如下结论。

(1)筒仓下部1/3hn范围内贮料的流动状态发生显著变化;贮料所表现出的流态与所处于筒仓的位置有关,越靠近漏斗部分,贮料卸料流动状态越接近漏斗流;卸料的发展进程为整体流动—管状流动—漏斗流。

(2)静态侧压力数值模拟值与试验值较为接近、动态侧压力分布规律一致,EDEM离散元卸料仿真模型动态侧压力研究具有可行性和合理性。

(3)卸料过程中,仓壁出现侧压力震荡增大现象是由于在仓体下部1/3hn范围内出现 “成拱—塌拱—再成拱—再塌拱”的动态拱过程。

(4)超压系数最大值出现于0.8 m测点深度处;相同位置范围内不同流速下的超压系数均高于规范中水平侧压力修正系数Ch;不同测点深度处,流速快工况的超压系数均大于流速慢工况的超压系数。