基于MBD和DEM耦合的新型CAE软件

王 博, 于哲舟, 袁 军, 付 宏, 于建群

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

播种机通常由机架、 开沟器、 排种器、 覆土和镇压机构等主要部件组成, 以作物种子为播种对象, 工作时还会与肥料和土壤颗粒发生接触作用. 播种机要通过多个部件的综合作业将种子播到土壤中适当的位置, 使种子在田间合理分布. 由于播种机工作过程复杂, 因此目前主要采用试验方法、 传统牛顿力学方法、 多刚体动力学(MBD)分析和光滑粒子法(SPH)等连续介质力学方法进行播种机工作过程的分析和优化设计.

试验方法[1], 即在满足试验条件的基础上进行多次性能和生产试验, 设计者需投入大量的时间和精力, 周期较长且最终得到的结果通常只适于试验的某种具体情况; 使用传统牛顿力学方法只能把土壤作为一个整体考虑, 而土壤由大量形状和大小不同的颗粒组成, 所以该方法无法分析群体中的每个颗粒与播种机零部件间的接触作用, 以及由于接触产生的对播深、 粒距均匀性的影响; 基于多刚体动力学分析的播种机设计研究[2], 重点关注播种机整个机构的运动过程, 且只分析数粒种子的运动, 与散粒群体的接触运动有较大差别; SPH方法[3]用一个或有限个质点描述种子颗粒或土壤颗粒, 未考虑散粒群体中不同形状颗粒的影响, 有一定的局限性. 因此, 上述几种方法均不能较好地解决播种机工作过程分析和优化设计的问题.

本文以播种机为例, 通过建立其多刚体动力学模型, 先利用系统动力学方程数值求解得到每一时步下播种机各刚体质心点的运动信息(位移、 速度和加速度), 再根据每一时刻颗粒与刚体的接触情况, 采用离散元法(DEM)求解颗粒与刚体接触时的相互作用力以及颗粒的运动轨迹, 并更新多刚体系统所受外力情况, 从而实现MBD与DEM的耦合, 在此基础上自主研发了多刚体动力学与离散元法耦合的分析计算软件AgriCAE(agricultural CAE), 为播种机、 深松机、 推土机、 挖掘机等工作过程分析和优化设计提供参考.

1 模型与方法

1.1 多刚体动力学模型

以播种机为例, 其机构运动示意图如图1所示. 由于播种机各构件(刚体)是在一个或几个相互平行的平面内运动, 因此本文采用平面笛卡尔坐标方法建立该系统的动力学模型[4]. 首先, 建立全局坐标系xoy, 如图2所示, 以播种机中某一构件j为例, 其质心点在全局系中的坐标为Dj(xj,yj), 以Dj为原点建立固结在该刚体上的局部坐标系x′o′y′,φj为全局坐标系x轴与局部坐标系x′轴正方向之间的夹角, 称为刚体j的姿态角, 规定逆时针方向为正方向, 从而刚体j的笛卡尔广义坐标矩阵为

若系统(播种机)中包含N个刚体, 每个刚体都按上述过程建立局部坐标系, 则有n=3N个坐标描述其位置, 系统的笛卡尔坐标列阵为

图1 播种机的主要部件及机构运动示意图

图2 全局坐标系及刚体j上局部坐标系

假设刚体j的质量为mj, 转动惯量为Jj, 则可定义刚体j的广义质量矩阵为

由于约束的存在, 各刚体的位置坐标并不是独立的, 设系统有m个位置约束方程且m

(1)

其中：Φi为第i个约束方程左部；q为笛卡尔广义坐标矩阵;t表示时间. 将式(1)对时间求一阶导数和二阶导数, 分别得速度和加速度约束方程为

(2)

(3)

(4)

其中λ是Lagrange乘子矩阵, 其确定了作用在系统上的约束反力和力矩. 将式(4)与式(3)联立可得[4]

(5)

1.2 离散元法模型

离散元法[5]可将土壤模型化为一定数量的形状和大小不同的颗粒集合, 主要由牛顿第二定律求解散粒群体中每个颗粒的运动, 以球形颗粒i为例, 有

(6)

其中：mi为颗粒i的质量；νi为颗粒i的平动速度；t表示时间； ∑fi为颗粒i所受的合外力, 主要由该颗粒自身重力、 颗粒与相邻颗粒碰撞及颗粒与机械部件(刚体)之间的接触作用力组成；Ii,ωi和∑Ti分别为颗粒i的转动惯量、 角速度及所受的合力矩. 接触两体之间作用力计算在颗粒的局部坐标系下进行, 分为法向和切向, 颗粒与边界(刚体)的接触如图3所示, 其中:X方向为法向;Y(Z)方向为切向;Z方向是垂直于纸面方向.

选用离散元法的接触力学模型时, 应尽量接近颗粒对象实际的材料属性和物理性质, 本文用线性黏弹性力学模型[6]进行计算, 其中法向作用力为

(7)

(8)

1.3 多刚体动力学与离散元法耦合模型

(9)

图3 颗粒与边界接触示意图

图4 刚体任意点P位置求解示意图

(10)

为计算接触点处的速度, 可将式(9)对时间求导

(11)

2 软件设计

本文基于MBD-DEM模型, 在VS2010 MFC平台上开发二维多刚体动力学模块, 并封装为动态链接库dll文件, 实现了与自主研发的离散元法(DEM)分析软件[7]的集成, 从而开发出具有自主知识产权、 用于农机产品设计的新型CAE软件----AgriCAE. 该软件集建模、 求解、 显示与分析等功能为一体, 可极大缩短机械的设计研发周期, 其结构和计算流程如图5所示.

图5 AgriCAE软件结构和流程

AgriCAE软件包括4个功能相对独立的部分, 按执行流程有CAD模型设计、 机械部件多刚体动力学和离散元法建模、 MBD与DEM耦合计算和性能分析子系统, 它们通过数据接口和文件耦合成一个整体. CAD模型设计子系统由CAD软件如Pro/E和UG组成, 用户通过人机交互绘制出播种机的三维CAD设计图. 机械部件多刚体动力学和离散元法建模子系统, 用户可采用两种方法建立播种机的分析模型, 即基于数据库的方法[8]和基于数据文件的方法[9]. MBD-DEM耦合计算子系统, 作为求解器, 负责组织模型数据, 计算散粒群体与播种机部件表面之间的接触作用、 刚体运动情况和颗粒运动轨迹等, 并将计算结果保存到一系列二进制文件中. 而性能分析子系统用于评价所设计的播种机机构的工作性能, 包含可视化模块(计算结果在屏幕端的动画演示)和曲线分析模块(播种机机构运动与力分析、 颗粒群体多项运动指标分析). 软件用户界面如图6所示.

图6 AgriCAE用户界面(多刚体建模)

3 软件测试及实例分析

3.1 多刚体动力学计算测试

以图1的播种机为实验对象, 由其CAD模型建立多刚体动力学分析模型, 检验软件中多刚体动力学单独计算模块的正确性. 分别采用AgriCAE软件和Adams软件对图1播种机横梁的加速度、 速度、 位移进行分析, 结果分别如图7～图9所示. 设置仿真边界运动时间为4 s, 计算时步为10-4s, 容许误差值为0.000 1. 由图7～图9可见, 两者数据曲线变化规律非常接近, 吻合程度良好, 从而验证了本文实现的多刚体动力学计算算法是正确的.

图7 采用AgriCAE软件和Adams软件分别得到的播种机工作时横梁质心的加速度曲线比较

图8 采用AgriCAE软件和Adams软件分别得到的播种机工作时横梁质心的速度曲线比较

图9 采用AgriCAE软件和Adams软件分别得到的播种机工作时横梁质心的位移曲线比较

3.2 多刚体动力学与离散元法耦合的仿真实例

为检验软件MBD-DEM耦合模型计算的可行性和合理有效性, 根据图1播种机的CAD模型(CAD设计图)建立多刚体动力学与离散元法耦合的分析模型, 并在此基础上对该播种机的开沟、 播种、 覆土和镇压过程进行仿真分析, 表1列出了耦合计算所用的仿真参数. 播种机工作时四杆仿形机构可使播种机横梁在地表上下浮动, 且可使播种深度保持一致. 弹簧可控制上下浮动的幅值, 特别是限制下仿形幅值, 开沟器可开出一定深度和宽度的沟形, 排种器把种子播到种沟后, 覆土器可将沟形两侧的土壤覆到种子上, 然后用镇压辊压实土壤, 图10为仿真效果截图. 通过仿真计算结果动态显示的观察可知, 播种机部件和散粒物料颗粒的运动情况与实际情况相近, 因此, 采用AgriCAE软件进行播种机工作过程仿真, 可实现播种机的开沟、 播种、 覆土和镇压的功能设计, 初步证明了本文研制软件的可行性和方法的合理性. 而这是用其他方法, 如传统牛顿力学、 多刚体动力学和光滑粒子法等无法实现的.

表1 离散元法与多刚体动力学耦合计算参数选取

图10 采用AgriCAE软件实现的播种机工作过程仿真分析

图11 采用AgriCAE软件实现的深松机工作过程仿真分析

由于深松机、 推土机、 挖掘机的工作过程与播种机相近, 因此都可采用多刚体动力学与离散元法耦合仿真其工作过程. 本文采用AgriCAE软件对一种深松机工作过程进行仿真分析, 图11为深松效果截图. 深松铲对一定深度的土壤进行疏松. 通过观察计算结果的仿真演示可知, AgriCAE软件可由深松机CAD模型(CAD设计图)实现深松机的工作过程仿真分析, 进一步证明了本文研发软件的可行性和有效性. 这也是用其他方法不能实现的.

综上所述, 针对播种机工作过程分析和优化设计中存在的问题, 本文研制了一种新型CAE软件——AgriCAE, 给出了软件中的多刚体动力学模型和离散元模型等关键技术及实现方法, 并通过播种机动力学分析测试与耦合仿真实例初步证明了该软件的可行性和有效性. 该软件可在播种机设计阶段, 由播种机的CAD模型实现播种机工作过程的模拟仿真和工作性能的分析评价, 对于其他类型或不符合设计需求的播种机, 只需重新建立其CAD模型, 更改结构尺寸方案或各项工作参数[10-11]即可再次执行分析过程, 为播种机、 深松机、 推土机、 挖掘机等工作过程分析和优化设计提供了一种新的参考方法.