基于离散元法玉米秸秆双层粘结模型力学特性研究

刘禹辰，张锋伟，宋学锋，王 锋，张方圆，李相周，曹晓庆

（甘肃农业大学机电工程学院，兰州 730070）

农作物秸秆作为重要生物质能源，全球年产量高达20 亿t。我国秸秆产量居世界之首，水稻、小麦和玉米秸秆资源量占比较大。2018 年全国主要农作物秸秆产量为82 850.24万t，其中玉米秸秆产量为26 746.51 万t，占比达32.28%[1]。我国对于农作物秸秆的利用始终坚持农用优先策略，大力倡导推进秸秆肥料化、饲料化和基料化利用，不断推动秸秆综合利用产业提质增效[2]。建立玉米秸秆离散元模型，开展机具耦合相关研究工作，可为饲草料加工机械关键零部件参数选取提供物料模型支撑，提高秸秆加工效率，提升利用率。

目前国内外学者针对秸秆二次利用已开展相应研究工作。张锋伟等早期研究结果表明，分析玉米秸秆破碎揉丝过程中受力，可表征玉米秸秆力学特性[3]。Chevanan 等通过直剪单元万能试验机对玉米秸秆作切碎试验和抗剪强度试验，结果表明粉碎后的玉米秸秆具有最高角摩擦，碎料之间无屈服强度[4]；张李娴通过万能试验机对玉米秸秆外表皮与内穰作力学特性试验，得出秸秆表皮法向粘结系数、切向粘结系数、法向刚度系数、切向刚度系数及玉米秸秆表皮弹性模量[5]；侯俊铭等采用旋转组合试验方案进行弯曲力学试验，结合有限元法得出不同生长部位纤维导管含量对茎秆抗弯强度的影响规律[6]。

赵吉坤等利用EDEM建立单层水稻秸秆离散元模型并标定相关粘结参数[7]。采用双层粘结模型标定秸秆力学特性的相关研究较少，通常采用各向同性（均质化）的方法建立数值模型，忽略秸秆表皮与内穰的力学特性差异。但玉米秸秆外表皮属于木质部，内穰属于纤维部，两者之间存在显著力学特性差异。通过对玉米秸秆外表皮与内穰的大量力学试验发现，内穰呈现各向同性，外表皮呈现各向异性，两者在弹性模量和剪切模量等相关参数上差异显著。

因此，建立离散元模型过程中区分相关参数差异，可提高建模准确性和可靠性。本文采用离散元Hertz-Mindlin with bonding 颗粒接触方法建立玉米秸秆双层粘结双峰分布数值模型，通过力学试验测定玉米秸秆相关特性参数，结合力学试验与虚拟仿真试验对比校正该模型，旨在为采用离散元法建立秸秆双层粘结模型，提供新建模思路，为青贮玉米秸秆饲料制造研究提供理论依据。

1 建立玉米秸秆双层粘结模型

1.1 bonding模型

离散元法（Discrete element method，DEM）是20 世纪70 年代由Cundall 和Starck 首先提出[8-9]，为分析求解复杂离散系统动力学问题提供新的数值方法[10]，目前已广泛应用于岩土工程、化学工程、农业工程等相关领域。采用离散元法直观反映秸秆模型在外力作用下粘结键破裂变化情况。农作物季节性较强，计算机数值模拟仿真一定程度上可替代繁杂的田间试验，省时、省力[11]，同时实现试验过程可视化，便于参数优化及控制，降低研究周期与成本。因玉米秸秆外表皮与内穰力学特性差异大，不能简单看为各向同性，采用有限元法分析存在一定局限性，所以在玉米秸秆离散元模型建立过程中，采用Hertz-Mindlin with bonding（bonding 模型）颗粒接触模型，该模型使相邻两颗粒在接触点处发生平行粘结，平行键作用效果相当于分布在球形基础颗粒截面上的一组弹簧，通过接触点力F1与力矩M1描述材料力学性能见图1[12]。该接触模型适用于模拟破碎断裂等问题，在模型建立过程中颗粒互相粘结形成粘结键，受外力作用时颗粒间粘结键发生断裂，通过粘结键破裂情况反映模型强度等相关力学特性。

图1 接触模型Fig.1 Contact model

图1 中Fi,b为颗粒A 作用于颗粒B 上合力（N）；、分别为法向、切向力矩（N）；ni、τi分别为法向和切向分量。Lb为颗粒A 与颗粒B 重叠量（mm）；Rb为粘结键半径（mm）。

1.2 理论计算

采用bonding 模型进行秸秆粘结，颗粒之间产生粘结键主要粘结参数有法向刚度系数Kn、切向刚度系数Ks、法向临界应力σ、切向临界应力γ和粘结半径Rj。

式中，va和vb为颗粒泊松比；Ea和Eb为颗粒弹性模量（MPa）；ra和rb为颗粒半径（mm）。

法向临界应力σ与切向临界应力γ计算公式如下：

式中，R为压缩面半径（mm）；F为临界压力（N）；c为秸秆内聚力（MPa）；φ为内摩擦角（°）。

将力学试验测得数据代入公式（1）～（4）中，得出玉米秸秆外表皮颗粒之间法向刚度系数Kn为2.55×108N·m-3，切向刚度系数Ks为1.71×108N·m-3，法向临界应力σ为6.34×107MPa，切向临界应力γ为7.45×106MPa，粘结半径Rj取0.36 mm，为颗粒半径的1/2。玉米秸秆内穰法向刚度系数Kn为11.41×107N·m-3，切向刚度系数Ks为7.61×106N·m-3，法向临界应力σ为4.9×104MPa，切向临界应力γ为4.25×104MPa 粘结半径Rj=0.6 mm。外表皮-内穰颗粒之间法向刚度系数Kn为1.02×108N·m-3，切向刚度系数Ks为6.82×107N·m-3，法向临界应力σ为6×107MPa，切向临界应力γ为6.5×106MPa，粘结半径取颗粒半径的1.2倍，Rj为0.48 mm。

1.3 建立玉米秸秆双层粘结离散元模型

模型几何尺寸依据试验所测样品尺寸确定，外表皮厚度，取值范围为0.82～2.86 mm，内穰取值为21.34～24.28 mm。为表征玉米秸秆外表皮木质部与内穰纤维部力学差异，在EDEM 中建立两种形状一致参数不同的散体颗粒，分别填充秸秆外表皮与秸秆内穰，参照所测外表皮厚度，设置外表皮颗粒粒径为0.3 mm，采用双峰分布的颗粒生成方式填充物料。在颗粒工厂生成颗粒时选择静态生成颗粒[13]，大颗粒生成后根据颗粒间空隙大小，自动生成双峰分布范围之内的小颗粒对空隙处作二次填充，该填充方式减小颗粒间空隙，提高玉米秸秆离散元模型填充密度，增加颗粒间粘结键数量，保证模型具有良好粘结强度。采用双峰分布后玉米秸秆外表皮平均粒径为0.312 mm。为保证两类颗粒粒径大小相差在一定范围之内，选取内穰颗粒粒径为0.5 mm，内穰颗粒平均粒径为0.528 mm，生成玉米秸秆内穰颗粒10 271 个，外表皮颗粒9 941 个。计算域网格大小为最小颗粒粒径的3.5倍，保证颗粒间接触计算可靠性。本研究中物料参数与接触参数总结见表1[14]。通过先填充后粘结方式，将玉米秸秆内穰颗粒-内穰颗粒、玉米秸秆外表皮颗粒-外表皮颗粒、玉米秸秆外表皮颗粒-内穰颗粒粘结，生成粘结键总数为108 305，其中内穰-内穰粘结键58 803、外表皮-外表皮粘结键32 501、外表皮-内穰粘结键17 001。建立双层玉米秸秆粘结接触模型，秸秆底面直径25 mm高25 mm。玉米秸秆外表皮颗粒以蓝色表示，内穰颗粒以黄色表示，外表皮颗粒之间粘结键以橙色表示，外表皮与内穰颗粒粘结键以绿色表示，内穰颗粒之间粘结键以深蓝色表示，见图2。

图2 玉米秸秆双层粘结双峰分布离散元模型Fig.2 Discrete element model of corn straw double-layer bonding and bimodal distribution

表1 离散元物料参数Table 1 Discrete element method material parameter

2 玉米秸秆力学特性试验

不同品种、成熟度和根节玉米秸秆在力学特性存在差异。本研究选取试验样本为大京九23 号蜡熟期青贮玉米秸秆，以离地第二根节开展相关力学试验。试验所用玉米秸秆内穰直径均值为23.16 mm。依据GB/T14699.1—2005《饲料采样》和GB/T6435—2014《饲料中水分和其他发挥性物质含量的测定》含水率测定方法[15-17]，采用HKSF-2型快速水分仪测得青贮玉米秸秆样品含水率为69.53%，青贮玉米秸秆收获要求含水率为65%～75%，样品符合收获要求。玉米秸秆外表皮和内穰力学特性试验在甘肃农业大学机电工程学院力学实验室完成，分别使用新三思实验设备有限公司生产电子万能材料试验台MTS-MODEL E43，试验台最大载荷为5 000 N，仪器精度为0.001 N 和SMS 公司生产TA.XT plus 超技质构仪，质构仪最大载荷为300 N，仪器精度为0.1 g。因质构仪最大载荷较小，无法完成青贮玉米秸秆轴向压缩试验，故使用两种仪器完成单轴力学压缩试验。对青贮玉米秸秆外表皮和内穰作拉伸、压缩和剪切试验，依据试验数据标定相关参数。

2.1 玉米秸秆外表皮与内穰拉伸试验

将青贮玉米秸秆外表皮与内穰均分成矩形条状样品，采用游标卡尺测量秸秆外表皮厚度与宽度。分成5 组，每组5次拉伸试验。试验时为防止出现秸秆与夹具相对滑动，影响测量准确性，将玉米秸秆外表皮与内穰按统一生长方向放入夹具中，调节夹具至最紧状态[5]，见图3（a）。为防止引伸计夹持造成样品表面出现应力集中现象，在不影响试验准确性的基础上对夹持处用一次性医用纱布做简单包裹，见图3（b）。将样品放入夹具后，在试验平台上将力和位移清零，取5次试验平均值作为试验值，青贮玉米秸秆外表皮与内穰拉伸试验见图3，拉伸结果见表2。试验结果表明青贮玉米秸秆外表皮弹性模量为982.52 MPa，青贮玉米秸秆内穰弹性模量为25.67 MPa。青贮玉米秸秆外表皮最大拉伸力为817.29 N，青贮玉米秸秆内穰最大拉伸力为47.66 N。青贮玉米秸秆外表皮抗拉强度为1 105 MPa，青贮玉米秸秆内穰抗拉强度为14.65 MPa。

表2 青贮玉米秸秆外表皮与内穰拉伸试验Table 2 Tensile test of silage corn straw skin and corn stalk pith

图3 青贮玉米秸秆外表皮与内穰拉伸试验Fig.3 Tensile test of silage corn straw skin and corn stalk pith

2.2 玉米秸秆压缩试验

对试验样品作预处理，样品高度为2.5 cm，采用游标卡尺测量样品直径。对切段面进行平整处理，使压缩力始终沿轴向。将青贮玉米秸秆与内穰径向和轴向分别放入试验台，开展单轴压缩试验。设置加载速度为2 mm·s-1将力和位移清零后重复作多组试验，取5次平均值作为试验值，轴向压缩结果见表3，径向压缩结果见表4。

表3 青贮玉米秸秆和内穰轴向压缩试验Table 3 Axial compression test of silage corn straw and corn stalk pith

表4 青贮玉米秸秆和内穰径向压缩试验Table 4 Radial compression test of corn stalk and corn stalk pith

2.3 剪切试验

青贮玉米秸秆外表皮因存在弧度，剪切过程中影响试验结果，因此处理样品时将其均分为更小份数，选取平整样品作剪切试验。通过控制器，校正质构仪力和高度。在试验方案中将剪切刀具位移速度设置为2 mm·min-1。取5次平均值作为试验值，剪切试验过程见图4，试验结果见表5。

表5 青贮玉米秸秆外表皮和内穰剪切结果Table 5 Shear results of silage corn straw skin and corn stalk pith

图4 青贮玉米秸秆外表皮与内穰剪切试验Fig.4 Shear tests of silage corn straw skin and corn stalk pith

随位移增加，剪切力逐渐上升后突然下降，此时剪切顺序依次为维管束、皮层和表皮。随维管束密度不断增加，剪切力也增加，直至秸秆发生断裂达到最大剪切力[19-20]。青贮玉米秸秆外表皮最大载荷为75.05 N，剪切强度为77.23 MPa，剪切模量为387.58 MPa。青贮玉米秸秆内穰的最大载荷为9.64 N，剪切强度为10.71 MPa，剪切模量为9.09 MPa。

玉米秸秆物理特性复杂，通过青贮玉米秸秆外表皮和内穰拉伸、压缩和剪切试验，结果表明玉米秸秆外表皮与内穰拉伸弹性模量、压缩强度和剪切模量均存在显著差异，玉米秸秆外表皮主要有表皮和管束组织，其力学特性呈现各向异性，管束组织结合强度可分为轴向和径向两个方向，径向为圆弧形状，容易折断；轴向结合强度高为其主要力学特性之一[21]。

3 离散元模型校核

采用离散元法进行力学特性仿真可直接观察玉米秸秆粘结键在受力时变化情况。本研究通过粘结键的断裂情况作为玉米秸秆受外力时破裂程度评定指标[22]。秸秆压缩和剪切仿真试验中，所用三维模型参照力学试验所用仪器大小，通过SolidWorks 绘制并将模型以.igs 格式导入到EDEM[23]。

3.1 压缩仿真试验

对秸秆模型分别进行单轴径向与轴向压缩，设置上压缩板以速度为2 mm·min-1匀速向下运动。因接触模型采用bonding 粘结模型，为减少试验误差Rayleigh时间步长设置为12%。颗粒粘结模型受力情况在后处理中导出，青贮玉米秸秆径向和横向压缩试验与仿真试验对比见图5。轴向压缩时，玉米秸秆双层粘结双峰分布离散元模型最大压缩力为1 846.1 N，与力学特性试验中青贮玉米秸秆轴向压缩受力大小一致，EDEM后处理中载荷与位移云图曲线变化趋势与力学特性试验结果一致（见图6）。在压缩过程中，轴向压缩试验粘结键断裂数为21 744 个，径向压缩试验粘结键断裂数为43 744 个，在相同试验条件下径向粘结键断裂数量高于轴向粘结键断裂，从仿真试验角度验证玉米秸秆轴向结合强度高是其主要力学特性之一。轴向外表皮-外表皮粘结键断裂12 993个断裂率66.85%，外表皮-内穰粘结键断裂5 593 个，断裂率50.45%，内穰-内穰粘结键断裂25 058个，断裂率74.21%，从粘结键断裂情况反映玉米秸秆受外力作用时外表皮抗压强度大于内穰抗压强度。玉米秸秆模型外表皮受力破裂情况与力学试验破裂情况一致，见图5（a）、图5（b），内穰颗粒与外表皮颗粒受力分离情况与力学试验也具有高度相似性。图5中粘结键主要断裂部位集中于玉米秸秆中心部位，见图5（c）、（d），断裂情况与力学特性试验破裂情况一致。通过对比径向与轴向力学试验与仿真试验结果，证明所建立玉米秸秆双层粘结双峰分布模型可靠，见图6。

图5 青贮玉米秸秆压缩试验与仿真试验效果Fig.5 Experimental and simulation test effect of silage corn straw compression

图6 青贮玉米秸秆轴向压缩力学试验与仿真试验结果Fig.6 Experimental and simulation results of silage corn straw axial compression

3.2 剪切仿真试验

设置剪切仿真试验刀具加载速度为2 mm·min-1，保持与力学试验运行速度一致。在刀具剪切玉米秸秆外表皮仿真过程中，可清楚观察到粘结键法向力和切向力变化情况，见图7。仿真所得最大剪切力为75.17 N，与力学试验数据一致。青贮玉米秸秆外表皮剪切离散元仿真受力情况见图8。

图7 外表皮粘结键受剪切力变化云图Fig.7 Cloud diagram of bond subjected to shear force

玉米秸秆外表皮离散元模型受剪切力作用时，颗粒间粘结键发生断裂，通过图8刀具剪切面的受力情况可以得出，在剪切力达到75.17 N 时，外表皮颗粒之间粘结键发生断裂，刀具持续向下运动，颗粒间粘结键断裂，但颗粒间粘结键非规律性粘接，各粘结键之间有平行、交叉等情况，在离散元剪切试验结果图中载荷出现先上升后下降再继续上升的情况，其主要原因为粘结键排布不规律。从剪切力折线图上升趋势，对比力学试验所获得相关结果和外表皮破裂情况，充分证明所建立的玉米秸秆外表皮粘结离散元模型与青贮玉米秸秆试样力学特性一致。

图8 青贮玉米秸秆外表皮剪切力学试验与仿真试验结果Fig.8 Shear mechanics test and simulation test results of silage corn straw outer skin

4 结论

对比秸秆离散元模型与秸秆试样在力学特性测试结果表现一致，证明基于离散元法建立双层玉米秸秆双峰分布粘结模型可靠。所建立玉米秸秆双层粘结双峰分布离散元模型可表征秸秆外表皮与内穰力学特性差异，相比现有单层玉米秸秆离散元模型具有更高精度。该模型通过减小颗粒间空隙度提高模型质量。虽然粘结键位置关系对剪切试验过程中载荷的细微变化有干扰，但在宏观上表现出良好力学特性。研究为采用离散元法探究秸秆不同部位力学特性差异提供新方法，并为进一步建立离散元模型，开展机具耦合互作，增速青贮玉米饲料制造的发展等相关研究提供理论指导。