基于离散元法的玉米脱粒过程分析

付 宏，吕 游，李艳双，于建群

（1.吉林大学计算机科学与技术学院，长春130012；2.吉林大学生物与农业工程学院，长春130022）

我国年产玉米近一亿五千万吨，其中大部分玉米都是收获晾晒后，再进行脱粒操作。因此，玉米脱粒机的使用非常广泛。但由于脱粒过程的复杂性，到目前为止，国内外对玉米脱粒过程的研究，大都采用试验方法、统计分析方法或传统连续介质力学的分析方法［1－3］。试验方法和统计分析方法费时费力，所得结果一般也不具有普遍意义，且不能洞察玉米脱粒过程的物理机理。传统连续介质力学的分析方法，只能分析脱粒过程中单个玉米穗、玉米芯或玉米籽粒的受力和运动情况，或把玉米穗、玉米芯和玉米籽粒当成一个整体来分析，这与实际脱粒过程中，玉米穗、玉米芯和玉米籽粒群体的相互作用和运动过程差别较大。

本文采用颗粒动力学方法——离散元法研究玉米的脱粒过程，建立了一种分析玉米脱粒过程的新方法和软件。

1 玉米穗和脱粒机的建模方法

玉米穗由玉米芯和生长在芯上的籽粒组成，且通常玉米穗成圆台形，而玉米籽粒截面近似梯形，如图1所示。本文采用球颗粒聚合方法，建立了玉米穗和玉米籽粒的分析模型如图2所示。

图1 实际玉米穗和玉米籽粒形状Fig.1 Real shapes of a corn ear and a corn grain

图2 玉米穗和玉米籽粒分析模型Fig.2 Analysis models of a corn ear and a corn grain

为分析玉米穗的脱粒过程，建立图3所示的全局坐标系x-y-z和局部坐标系X-Y-Z。其中全局坐标系是大地坐标系，局部坐标系的原点与玉米穗的质心重合，X轴与玉米穗轴线重合，且与某一玉米籽粒中心重合，Z轴方向符合右手螺旋法则。为了模拟玉米籽粒与玉米芯之间的连接，采用一个球（称连接球）模拟籽粒果柄，且籽粒果柄与玉米籽粒是一个整体，如图4所示，玉米籽粒脱落后，连接球被删除。玉米芯组成单元之间的连接也采用同样的连接球方法处理。

图3 全局坐标系和局部坐标系Fig.3 Global coordinate system and local coordinate system

在分析玉米脱粒过程时，首先在全局坐标系的玉米穗入口随机产生一点，作为待生成玉米穗的质心，然后以该质心为坐标原点，产生局部坐标系X-Y-Z并生成一个玉米穗，下一时步再随机产生一点，并用同样方法生成另一个玉米穗，如此反复。其中玉米籽粒分析模型的高度、厚度和上、下底宽度尺寸、玉米穗分析模型的长度和大小头直径，均可通过程序界面由输入参数控制。

图4 玉米芯、玉米籽粒和籽粒果柄Fig.4 Corn cob，corn grain and kernel stem

玉米脱粒机的建模方法采用本项目组提出的基于CAD模型的边界建模方法［4］，即把与玉米穗接触的脱粒机中零件表面离散成图元，如平面、球面、柱面、锥面等，并添加运动属性和材料特性，由此建立脱粒机的离散元法分析模型。以滚筒式玉米脱粒机为例，其三维CAD模型和三维离散元法分析模型如图5和图6所示。

图5 滚筒式玉米脱粒机的三维CAD模型Fig.5 3D CAD model of a drum type corn sheller

图6 滚筒式玉米脱粒机的三维离散元法分析模型Fig.6 Analysis model of a drum type corn sheller for 3D DEM

2 脱粒过程的作用力计算方法

玉米脱粒过程中的作用力包括籽粒与籽粒之间、籽粒与玉米芯之间以及籽粒与脱粒部件之间的作用力，还包括玉米芯与玉米芯之间和玉米芯与脱粒部件之间的作用力。在采用颗粒聚合体方法建立玉米穗分析模型时，当籽粒未脱落前，籽粒与其生长的玉米芯之间的作用力为连接力（吸引力），同一玉米芯组成单元之间的作用力也为连接力（吸引力）；而籽粒与其他籽粒之间、籽粒与其生长的玉米芯以外的玉米芯之间、籽粒与脱粒部件之间的作用力为接触斥力；玉米芯与其他玉米芯间、玉米芯与脱粒部件之间的作用力也为接触斥力。由于玉米籽粒和玉米芯的建模方法均采用球颗粒聚合的方法，因此计算籽粒和玉米芯上的作用力，可采用球颗粒的计算方法，这样就避免了求解非线性方程组。下面以某一玉米籽粒P为例（见图4），介绍接触作用力的计算方法。

首先在全局坐标系x-y-z中建立如图7所示的局部坐标系X′-Y′-Z′。其中X′轴为玉米籽粒P组成球i的球心与另一个接触籽粒或玉米芯组成球j的球心连线，Y′轴平行于全局坐标系x-y平面，Z′轴符合右手螺旋法则。在X′-Y′-Z′局部坐标系中，当i球与j球的中心距小于两球的半径之和时，则两球接触，其法向叠合量δn为

式中：Ri、Rj、（xi，yi，zi）和（xj，yj，zj）分别为球i和j的半径及中心点在全局坐标系x-y-z中的坐标。当采用线性黏弹性力学模型计算两球的接触斥力时，两球间的法向接触作用力为

图7 两颗粒组成球接触时的局部坐标系Fig.7 Local coordinate system when two balls contact

当球i与脱粒部件某一表面接触时，如以平面为例，局部坐标系X′轴为球i的球心到该平面的垂线，Y′和Z′轴的确定方法及法向和切向接触作用力的计算方法与上述相同。

采用线性黏弹性力学模型计算球i与球j的连接力的方法与上述相同，但刚度系数和阻尼系数选取不同，计算斥力时选取接触刚度系数和接触阻尼系数，计算连接力时选取连接刚度系数和连接阻尼系数。而且在籽粒没脱落前，不采用库仑－莫尔准则对切向力进行修正。

当采用上述方法计算某籽粒与其生长的玉米芯之间的法向及切向连接力和玉米芯组成单元之间的法向及切向连接力时，还需分别与由试验测定的连接作用力最大值作比较，如果计算的连接力大于试验测定的连接力最大值，则玉米籽粒从玉米芯上脱离——脱粒，或玉米芯组成单元之间分离——破碎。

当采用上述方法计算作用在玉米穗、脱落的玉米籽粒或破碎的玉米芯单元上每个组成球的法向和切向作用力后，还需把作用力分别移到玉米穗、玉米籽粒或玉米芯单元的质心上，并变换到全局坐标系下，然后再求出作用在玉米穗、玉米籽粒或破碎玉米芯单元质心上的合力和合力矩。

3 颗粒运动计算方法

在玉米脱粒过程中，颗粒运动包括玉米穗的运动、脱落籽粒的运动、破碎玉米芯单元的运动和未脱落籽粒相对玉米芯的运动。

通过上面的分析计算后，可求得作用在玉米穗、脱落玉米籽粒和破碎玉米芯单元质心上的在t时刻下沿全局坐标系3个坐标轴方向所受的合力，继而可得t＋Δt时刻颗粒质心在全局坐标系3个坐标轴的新位置 （x（t＋Δt），y（t＋Δt），z（t＋Δt））为

玉米穗、脱落玉米籽粒和破碎玉米芯单元在t时刻的转动角速度满足如下欧拉动力学方程［5］

求解上述方程即可得到颗粒在局部坐标系X-Y-Z中的角位移及新的局部坐标系，限于篇幅求解过程从略。

图8 计算玉米籽粒转动的局部坐标系Fig.8 Local coordinate system when calculated corn grain rotating

计算未脱落籽粒相对于玉米芯的运动，还需建立局部坐标系X″-Y″-Z″，参见图3和图4。其中X″轴与X轴重合，Y″轴为被分析玉米籽粒质心与玉米芯中心连线，当被分析籽粒为P时，Y″轴与Y轴重合，当分析其他籽粒时，Y″轴需变化；Z″轴符合右手螺旋法则。在X″-Y″-Z″局部坐标系中，求出籽粒上作用的合外力，并求出籽粒与玉米芯的连接力，根据合外力与连接力即可求出未脱落籽粒相对于玉米芯的新位置并对其更新。

通过上述计算方法按时步迭代，即可求出玉米穗、脱落玉米籽粒、破碎玉米芯单元和未脱落籽粒相对于玉米芯在每一时刻所受的力、力矩、运动速度和新位置，由此即可分析玉米的脱粒过程。

4 软件设计和实例验证

在上述工作的基础上，本文设计了玉米脱粒过程分析软件，并实现了与CAD软件的集成，从而开发出一种集设计和性能分析评价为一体的玉米脱粒过程分析和脱粒机优化设计软件［6－7］，其结构如图9所示。其流程为由CAD软件进行玉米脱粒机的设计，然后由边界建模模块自动建立脱粒机的离散元法分析模型，接着即可进行脱粒过程的离散元法分析计算，计算结束后可根据计算结果文件进行脱粒过程的性能分析，如脱净率、籽粒损失率和籽粒沿脱粒机轴向分布等，还可以动态显示玉米的脱粒过程。

图9 玉米脱粒过程和脱粒机优化设计新方法及其软件的结构Fig.9 Proposed method and software structure for analysis of corn threshing process and optimization of thresher

图10和图11为采用上述方法及其软件，由CAD模型实现的两种结构滚筒式玉米脱粒机工作过程的二维离散元法仿真分析。图12为采用上述方法及其软件，由图5所示玉米滚筒式脱粒机的三维CAD模型实现的玉米脱粒过程的三维离散元法仿真分析。从图12可以看出，随着仿真计算时间的延续，玉米籽粒逐渐从玉米穗上脱落并从脱粒机凹板上的孔漏出。图13为软件给出的破碎率、脱净率和脱落籽粒沿脱粒机轴线分布的分析结果。由此初步证明了本文方法和软件的可行性。有关仿真和试验的详细对比分析将另文撰述。

图10 由CAD模型实现玉米脱粒过程的二维仿真分析Fig.10 2D DEM simulation and analysis by CAD model

图11 改变CAD模型实现的玉米脱粒过程二维仿真分析Fig.11 2D DEM simulation and analysis by changing CAD model

本文方法和软件的优点是：①在设计阶段通过修改脱粒机的CAD模型、玉米籽粒、玉米芯和玉米穗的分析模型、连接和接触作用的力学模型、离散元法计算参数，能分析不同品种玉米、不同工况、不同原理、不同结构和尺寸的玉米脱粒机的性能，由此实现玉米脱粒机结构方案和尺寸参数的优化；②通过脱粒机的CAD模型，能进行脱粒过程的动态仿真，由此分析玉米脱粒机的工作机理或工作过程，还可以发明新原理和新结构的玉米脱粒机，这是现有玉米脱粒机的研究和设计方法不能做到的。

图12 玉米脱粒过程的三维离散元法仿真分析Fig.12 3D DEM simulation and analysis of corn threshing process

图13 本文软件给出的脱粒过程分析结果Fig.13 Analysis results of threshing process by software

5 结束语

根据玉米脱粒机研究和设计中存在的问题，本文提出采用离散元法研究玉米的脱粒过程，由脱粒机的CAD模型建立了脱粒机的离散元法分析模型，采用颗粒聚合体方法建立了玉米穗的分析模型，采用接触和连接力学模型计算脱粒过程中的作用力，研究了玉米脱粒过程的离散元法计算方法，在此基础上建立了一种分析玉米脱粒过程的新方法和软件。通过对该软件的实例验证，初步证明了新方法的可行性，为玉米脱粒过程的分析和脱粒机的优化提供了一种新思路。

［1］Petkevichius S，Shpokas L，Kutzbach H D.Investigation of the maize ear threshing process［J］.Biosystems Engineering，2008，99：532-539.

［2］Miu P I，Kutzbach H D.Modeling and simulation of grain threshing and separation in threshing units-Part I［J］.Computers and Electronics in Agriculture，2008，60：96-104.

［3］毛欣，衣淑娟.轴流脱粒、分离机理及仿真研究的发展探讨［J］.农机化研究，2007，10：202-205.

Mao Xin，Yi Shu-juan.The developmental discussion of mechanism research and simulation of the axial flow threshing mechanism［J］.Journal of Agricultural Mechanization Research，2007，10：202-205.

［4］于建群，付宏，刘振宇，等.基于CAD模型的离散元法边界建模方法：中国，200510016835.7［P］. 2006-07-28.

［5］Lin X，Ng T T.A three dimensional discrete element model using arrays of ellipsoids［J］.Geotechnique，1997，47（2）：319-329.

［6］李艳双.颗粒聚合体的三维离散元法算法研究与软件设计［D］.长春：吉林大学，2011.

Li Yan-shuang.Research of algorithm for particles agglomerate in 3D DEM and design of the software［D］.Changchun：Jilin University，2011.

［7］吴玄辰.基于二维离散元法的玉米脱粒过程仿真分析［D］.长春：吉林大学生物与农业工程学院，2011.

Wu Xuan-chen.Simulation and analysis of maize threshing based on 2D DEM［D］.Changchun：Jilin University，College of Biological and Agricultural Engineering，2011.