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

(1. 吉林大學(xué) 計算機(jī)科學(xué)與技術(shù)學(xué)院, 長春 130012； 2. 吉林大學(xué) 生物與農(nóng)業(yè)工程學(xué)院, 長春 130022)

播種機(jī)通常由機(jī)架、 開溝器、 排種器、 覆土和鎮(zhèn)壓機(jī)構(gòu)等主要部件組成, 以作物種子為播種對象, 工作時還會與肥料和土壤顆粒發(fā)生接觸作用. 播種機(jī)要通過多個部件的綜合作業(yè)將種子播到土壤中適當(dāng)?shù)奈恢? 使種子在田間合理分布. 由于播種機(jī)工作過程復(fù)雜, 因此目前主要采用試驗方法、 傳統(tǒng)牛頓力學(xué)方法、 多剛體動力學(xué)(MBD)分析和光滑粒子法(SPH)等連續(xù)介質(zhì)力學(xué)方法進(jìn)行播種機(jī)工作過程的分析和優(yōu)化設(shè)計.

試驗方法[1], 即在滿足試驗條件的基礎(chǔ)上進(jìn)行多次性能和生產(chǎn)試驗, 設(shè)計者需投入大量的時間和精力, 周期較長且最終得到的結(jié)果通常只適于試驗的某種具體情況; 使用傳統(tǒng)牛頓力學(xué)方法只能把土壤作為一個整體考慮, 而土壤由大量形狀和大小不同的顆粒組成, 所以該方法無法分析群體中的每個顆粒與播種機(jī)零部件間的接觸作用, 以及由于接觸產(chǎn)生的對播深、 粒距均勻性的影響; 基于多剛體動力學(xué)分析的播種機(jī)設(shè)計研究[2], 重點關(guān)注播種機(jī)整個機(jī)構(gòu)的運動過程, 且只分析數(shù)粒種子的運動, 與散粒群體的接觸運動有較大差別; SPH方法[3]用一個或有限個質(zhì)點描述種子顆粒或土壤顆粒, 未考慮散粒群體中不同形狀顆粒的影響, 有一定的局限性. 因此, 上述幾種方法均不能較好地解決播種機(jī)工作過程分析和優(yōu)化設(shè)計的問題.

本文以播種機(jī)為例, 通過建立其多剛體動力學(xué)模型, 先利用系統(tǒng)動力學(xué)方程數(shù)值求解得到每一時步下播種機(jī)各剛體質(zhì)心點的運動信息(位移、 速度和加速度), 再根據(jù)每一時刻顆粒與剛體的接觸情況, 采用離散元法(DEM)求解顆粒與剛體接觸時的相互作用力以及顆粒的運動軌跡, 并更新多剛體系統(tǒng)所受外力情況, 從而實現(xiàn)MBD與DEM的耦合, 在此基礎(chǔ)上自主研發(fā)了多剛體動力學(xué)與離散元法耦合的分析計算軟件AgriCAE(agricultural CAE), 為播種機(jī)、 深松機(jī)、 推土機(jī)、 挖掘機(jī)等工作過程分析和優(yōu)化設(shè)計提供參考.

1 模型與方法

1.1 多剛體動力學(xué)模型

以播種機(jī)為例, 其機(jī)構(gòu)運動示意圖如圖1所示. 由于播種機(jī)各構(gòu)件(剛體)是在一個或幾個相互平行的平面內(nèi)運動, 因此本文采用平面笛卡爾坐標(biāo)方法建立該系統(tǒng)的動力學(xué)模型[4]. 首先, 建立全局坐標(biāo)系xoy, 如圖2所示, 以播種機(jī)中某一構(gòu)件j為例, 其質(zhì)心點在全局系中的坐標(biāo)為Dj(xj,yj), 以Dj為原點建立固結(jié)在該剛體上的局部坐標(biāo)系x′o′y′,φj為全局坐標(biāo)系x軸與局部坐標(biāo)系x′軸正方向之間的夾角, 稱為剛體j的姿態(tài)角, 規(guī)定逆時針方向為正方向, 從而剛體j的笛卡爾廣義坐標(biāo)矩陣為

若系統(tǒng)(播種機(jī))中包含N個剛體, 每個剛體都按上述過程建立局部坐標(biāo)系, 則有n=3N個坐標(biāo)描述其位置, 系統(tǒng)的笛卡爾坐標(biāo)列陣為

圖1 播種機(jī)的主要部件及機(jī)構(gòu)運動示意圖

圖2 全局坐標(biāo)系及剛體j上局部坐標(biāo)系

假設(shè)剛體j的質(zhì)量為mj, 轉(zhuǎn)動慣量為Jj, 則可定義剛體j的廣義質(zhì)量矩陣為

由于約束的存在, 各剛體的位置坐標(biāo)并不是獨立的, 設(shè)系統(tǒng)有m個位置約束方程且m

(1)

其中：Φi為第i個約束方程左部；q為笛卡爾廣義坐標(biāo)矩陣;t表示時間. 將式(1)對時間求一階導(dǎo)數(shù)和二階導(dǎo)數(shù), 分別得速度和加速度約束方程為

(2)

(3)

(4)

其中λ是Lagrange乘子矩陣, 其確定了作用在系統(tǒng)上的約束反力和力矩. 將式(4)與式(3)聯(lián)立可得[4]

(5)

1.2 離散元法模型

離散元法[5]可將土壤模型化為一定數(shù)量的形狀和大小不同的顆粒集合, 主要由牛頓第二定律求解散粒群體中每個顆粒的運動, 以球形顆粒i為例, 有

(6)

其中：mi為顆粒i的質(zhì)量；νi為顆粒i的平動速度；t表示時間； ∑fi為顆粒i所受的合外力, 主要由該顆粒自身重力、 顆粒與相鄰顆粒碰撞及顆粒與機(jī)械部件(剛體)之間的接觸作用力組成；Ii,ωi和∑Ti分別為顆粒i的轉(zhuǎn)動慣量、 角速度及所受的合力矩. 接觸兩體之間作用力計算在顆粒的局部坐標(biāo)系下進(jìn)行, 分為法向和切向, 顆粒與邊界(剛體)的接觸如圖3所示, 其中:X方向為法向;Y(Z)方向為切向;Z方向是垂直于紙面方向.

選用離散元法的接觸力學(xué)模型時, 應(yīng)盡量接近顆粒對象實際的材料屬性和物理性質(zhì), 本文用線性黏彈性力學(xué)模型[6]進(jìn)行計算, 其中法向作用力為

(7)

(8)

1.3 多剛體動力學(xué)與離散元法耦合模型

(9)

圖3 顆粒與邊界接觸示意圖

圖4 剛體任意點P位置求解示意圖

(10)

為計算接觸點處的速度, 可將式(9)對時間求導(dǎo)

(11)

2 軟件設(shè)計

本文基于MBD-DEM模型, 在VS2010 MFC平臺上開發(fā)二維多剛體動力學(xué)模塊, 并封裝為動態(tài)鏈接庫dll文件, 實現(xiàn)了與自主研發(fā)的離散元法(DEM)分析軟件[7]的集成, 從而開發(fā)出具有自主知識產(chǎn)權(quán)、 用于農(nóng)機(jī)產(chǎn)品設(shè)計的新型CAE軟件----AgriCAE. 該軟件集建模、 求解、 顯示與分析等功能為一體, 可極大縮短機(jī)械的設(shè)計研發(fā)周期, 其結(jié)構(gòu)和計算流程如圖5所示.

圖5 AgriCAE軟件結(jié)構(gòu)和流程

AgriCAE軟件包括4個功能相對獨立的部分, 按執(zhí)行流程有CAD模型設(shè)計、 機(jī)械部件多剛體動力學(xué)和離散元法建模、 MBD與DEM耦合計算和性能分析子系統(tǒng), 它們通過數(shù)據(jù)接口和文件耦合成一個整體. CAD模型設(shè)計子系統(tǒng)由CAD軟件如Pro/E和UG組成, 用戶通過人機(jī)交互繪制出播種機(jī)的三維CAD設(shè)計圖. 機(jī)械部件多剛體動力學(xué)和離散元法建模子系統(tǒng), 用戶可采用兩種方法建立播種機(jī)的分析模型, 即基于數(shù)據(jù)庫的方法[8]和基于數(shù)據(jù)文件的方法[9]. MBD-DEM耦合計算子系統(tǒng), 作為求解器, 負(fù)責(zé)組織模型數(shù)據(jù), 計算散粒群體與播種機(jī)部件表面之間的接觸作用、 剛體運動情況和顆粒運動軌跡等, 并將計算結(jié)果保存到一系列二進(jìn)制文件中. 而性能分析子系統(tǒng)用于評價所設(shè)計的播種機(jī)機(jī)構(gòu)的工作性能, 包含可視化模塊(計算結(jié)果在屏幕端的動畫演示)和曲線分析模塊(播種機(jī)機(jī)構(gòu)運動與力分析、 顆粒群體多項運動指標(biāo)分析). 軟件用戶界面如圖6所示.

圖6 AgriCAE用戶界面(多剛體建模)

3 軟件測試及實例分析

3.1 多剛體動力學(xué)計算測試

以圖1的播種機(jī)為實驗對象, 由其CAD模型建立多剛體動力學(xué)分析模型, 檢驗軟件中多剛體動力學(xué)單獨計算模塊的正確性. 分別采用AgriCAE軟件和Adams軟件對圖1播種機(jī)橫梁的加速度、 速度、 位移進(jìn)行分析, 結(jié)果分別如圖7～圖9所示. 設(shè)置仿真邊界運動時間為4 s, 計算時步為10-4s, 容許誤差值為0.000 1. 由圖7～圖9可見, 兩者數(shù)據(jù)曲線變化規(guī)律非常接近, 吻合程度良好, 從而驗證了本文實現(xiàn)的多剛體動力學(xué)計算算法是正確的.

圖7 采用AgriCAE軟件和Adams軟件分別得到的播種機(jī)工作時橫梁質(zhì)心的加速度曲線比較

圖8 采用AgriCAE軟件和Adams軟件分別得到的播種機(jī)工作時橫梁質(zhì)心的速度曲線比較

圖9 采用AgriCAE軟件和Adams軟件分別得到的播種機(jī)工作時橫梁質(zhì)心的位移曲線比較

3.2 多剛體動力學(xué)與離散元法耦合的仿真實例

為檢驗軟件MBD-DEM耦合模型計算的可行性和合理有效性, 根據(jù)圖1播種機(jī)的CAD模型(CAD設(shè)計圖)建立多剛體動力學(xué)與離散元法耦合的分析模型, 并在此基礎(chǔ)上對該播種機(jī)的開溝、 播種、 覆土和鎮(zhèn)壓過程進(jìn)行仿真分析, 表1列出了耦合計算所用的仿真參數(shù). 播種機(jī)工作時四桿仿形機(jī)構(gòu)可使播種機(jī)橫梁在地表上下浮動, 且可使播種深度保持一致. 彈簧可控制上下浮動的幅值, 特別是限制下仿形幅值, 開溝器可開出一定深度和寬度的溝形, 排種器把種子播到種溝后, 覆土器可將溝形兩側(cè)的土壤覆到種子上, 然后用鎮(zhèn)壓輥壓實土壤, 圖10為仿真效果截圖. 通過仿真計算結(jié)果動態(tài)顯示的觀察可知, 播種機(jī)部件和散粒物料顆粒的運動情況與實際情況相近, 因此, 采用AgriCAE軟件進(jìn)行播種機(jī)工作過程仿真, 可實現(xiàn)播種機(jī)的開溝、 播種、 覆土和鎮(zhèn)壓的功能設(shè)計, 初步證明了本文研制軟件的可行性和方法的合理性. 而這是用其他方法, 如傳統(tǒng)牛頓力學(xué)、 多剛體動力學(xué)和光滑粒子法等無法實現(xiàn)的.

表1 離散元法與多剛體動力學(xué)耦合計算參數(shù)選取

圖10 采用AgriCAE軟件實現(xiàn)的播種機(jī)工作過程仿真分析

圖11 采用AgriCAE軟件實現(xiàn)的深松機(jī)工作過程仿真分析

由于深松機(jī)、 推土機(jī)、 挖掘機(jī)的工作過程與播種機(jī)相近, 因此都可采用多剛體動力學(xué)與離散元法耦合仿真其工作過程. 本文采用AgriCAE軟件對一種深松機(jī)工作過程進(jìn)行仿真分析, 圖11為深松效果截圖. 深松鏟對一定深度的土壤進(jìn)行疏松. 通過觀察計算結(jié)果的仿真演示可知, AgriCAE軟件可由深松機(jī)CAD模型(CAD設(shè)計圖)實現(xiàn)深松機(jī)的工作過程仿真分析, 進(jìn)一步證明了本文研發(fā)軟件的可行性和有效性. 這也是用其他方法不能實現(xiàn)的.

綜上所述, 針對播種機(jī)工作過程分析和優(yōu)化設(shè)計中存在的問題, 本文研制了一種新型CAE軟件——AgriCAE, 給出了軟件中的多剛體動力學(xué)模型和離散元模型等關(guān)鍵技術(shù)及實現(xiàn)方法, 并通過播種機(jī)動力學(xué)分析測試與耦合仿真實例初步證明了該軟件的可行性和有效性. 該軟件可在播種機(jī)設(shè)計階段, 由播種機(jī)的CAD模型實現(xiàn)播種機(jī)工作過程的模擬仿真和工作性能的分析評價, 對于其他類型或不符合設(shè)計需求的播種機(jī), 只需重新建立其CAD模型, 更改結(jié)構(gòu)尺寸方案或各項工作參數(shù)[10-11]即可再次執(zhí)行分析過程, 為播種機(jī)、 深松機(jī)、 推土機(jī)、 挖掘機(jī)等工作過程分析和優(yōu)化設(shè)計提供了一種新的參考方法.