基于DEM-CFD耦合的辣椒清選仿真研究

張德俊， 張小明， 吳荻， 林蜀云， 張太華， 徐衛平*

（1.貴州師范大學機械與電氣工程學院，貴陽 550001； 2.貴州省山地農業機械研究所，貴陽 550007）

辣椒產業是貴州最具代表性的傳統經濟優勢產業之一，是貴州現代農業的重要經濟支柱，同時，也是農民脫貧、創業增收的重要途徑[1]。2020 年，貴州省辣椒種植面積36.3 hm2，產量724 萬t，產值242億元,生產、加工、銷售規模居全國首位，貴州正從“辣椒大省”向“辣椒強省”邁進[2]。由于貴州山區辣椒種植面積大、田間地塊小及人工收獲成本高等問題，設計并制造了山地自走式履帶辣椒收獲機，通過機收降低了人工收獲成本，但機收的辣椒含雜率較高，果實清選問題亟待解決。目前，谷物、花生等農產品已有成熟的風選設備并投入應用，辣椒收獲機上雖然設有風選裝置，但參數的選擇主要是憑經驗，辣椒風選效果不好。

清選方式一般有兩種：一種是氣流清選，主要是利用混合物的不同懸浮性質；另一種是利用氣流和振動篩的聯合作用來實現清選。氣流清選是以空氣為主要作用介質，在氣流的作用下，根據顆粒密度和顆粒大小對各種顆粒材料進行空氣分選的過程，基本原理為通過氣流將較輕的顆粒水平的吹向遠處，而較重的顆粒吹到近處[3]。

近年來，隨著計算機技術的飛速發展，離散元法（discrete element method，DEM）和計算流體力學（computational fluid dynamics，CFD）在風選設備的分析中應用廣泛。敬志臣等[4]通過Fluent 軟件對小麥進行研究，采用旋風分離筒對分離過程進行仿真分析；李興凱等[5]利用EDEM（extended DEM）軟件對小區小麥聯合收割機的脫粒機構進行仿真分析。使用CFD 技術可以模擬氣流場，使用DEM 可以模擬粒子碰撞和接觸過程。但是，它們只能以單一方式模擬清選過程，不能準確、全面地反映風選室中作物在氣流、容器壁和顆粒碰撞接觸作用下的復雜運動，DEM-CFD 耦合模擬技術將氣流場、風選室和顆粒作為整體系統進行研究，可以準確分析出風選過程。張鳳奎等[6]利用CFDDEM 耦合方法研究了紅棗的負壓氣吸拾撿仿真并進行了分析；郭柄江等[7]通過EDEM 與Fluent 耦合方法研究了谷物風選，采用氣吹與風篩組合的風選方式，不同風速參數下稻米清選裝置工作性能的指標通過風選后稻米損失率和含雜率來進行衡量；丁賀賀等[8]通過DEM-CFD 耦合方法研究了谷物在垂直氣流吹氣下的仿真分離。因此，針對辣椒收獲機收獲的辣椒含雜率比較高的問題，設計了一種辣椒清選分離裝置，采用DEM-CFD耦合方法，對朝天椒風選（吹氣式）分離過程進行模擬仿真，旨在為設計及優化辣椒收獲機的辣椒風選分離裝置提供參考。

1 材料與方法

1.1 風選裝置、原理及阻風效應

1.1.1風選裝置 辣椒風選分離裝置由進風口、進料口、清選室、出口Ⅰ、出口Ⅱ、出口Ⅲ組成，模型結構如圖1 所示。使用SolidWorks 三維設計軟件建模，將風選裝置模型輸出為XT 格式并導入ICEM CFD 網格劃分軟件中，添加進風口、物料入口、辣椒出口Ⅰ、辣椒莖稈及辣椒葉出口Ⅲ，設置完成后劃分網格如圖2所示。

圖1 風選裝置Fig. 1 Wind-sorting device

圖2 風選裝置網格Fig. 2 Grid of wind sorting device

1.1.2風選原理 在水平清選氣流中，4 種力作用于物料顆粒，即水平風力FD、重力FG、阻力FR和浮力F′（圖3）?？諝饷芏群艿停锪项w粒的體積也很小，因此相應的浮力F′很小，對氣流中其他物料顆粒的徑向運動速度影響不大，所以浮力F′可以忽略，因此，垂直方向上合力FGr與水平風力FD的合力等于F，合力F與風力FD的夾角為β。此時，物料顆粒沿合力F方向運動，軌跡大致呈拋物線。物料顆粒的最大運動距離與夾角β的值有關：夾角β越小，物料顆粒運動得越遠，夾角β越大，顆粒運動距離越近。其中，Va為平均風速，V0為顆粒沉降速度，V為顆粒在合力作用下整體的速度。

圖3 水平氣流中的顆粒受力分析Fig. 3 Force analysis of particles in horizontal airflow

顆粒沉降速度V0可表示如下[9]。

式中，d為顆粒的粒子直徑，mm；ρs為顆粒密度，kg·m-3；ρ為風選的氣流密度，kg·m-3；φ為阻力系數；g為顆粒重力加速度，m·s-2。

用沉降速度V0和風速Va的比值表示β的正切值。

從式(2)可以看出，在顆粒水平風速保持不變、顆粒直徑大小相近的情況下，低密度顆粒的tanβ值較低，顆?？梢赃\動得更遠，而高密度顆粒的tanβ值較高，顆粒運動得就會更近，不同的物料顆粒可通過氣流按其密度、形狀和大小進行分離。

1.1.3阻風效應 將多顆粒連續運動模型簡化為顆粒疊加的粒子組模型（圖4）。從圖4 可明顯看出，當氣流中前端的流動力作用于前端粒子組時，氣流受到前端粒子組的阻礙而形成縱向流，沿橫向流過前端粒子組上、下表面。前端粒子組對后端粒子組的阻力屏障形成側流，氣流中的流動力將繞過后端粒子組。前端粒子群抑制氣流方向的向前運動，風力不作用于后粒子群，產生阻礙效應（即阻力效應）影響后端粒子組的整體運動、受力和運動軌跡的改變。

圖4 風選顆粒阻風效應簡化原理Fig. 4 Simplified principle of wind-selected particle windblocking effect

1.2 DEM-CFD耦合模型構建

1.2.1Fluent 流體模型 采用歐拉耦合模型對辣椒及其雜質的風選分離過程進行模擬，流體體積分數項和流體運動微分方程表示如下[10]。

式中，t為時間；ε為流體體積分數項；u為流速；g為重力加速度；ρ為流體密度；S為動量源項；μ為粘滯系數；P為流體微元上的壓強；?為哈密頓微分算子。

動量源項S為氣流作用在單元網格內部氣流阻力的總和[11]。

式中，V為單元網格體積；Fi為第i個顆粒對氣流的阻力。

1.2.2EDEM 固體顆粒接觸碰撞模型 在EDEM軟件中對顆粒力學模型間的接觸力學行為和相互作用力進行仿真分析。由于顆粒間的速度隨接觸力而發生變化，本文采用Hertz-Minderlin 接觸理論和軟球干接觸模型[12]，根據牛頓第二定律獲得第i個顆粒的運動方程[13]。

式中，g為重力加速度；P為顆粒在氣流中運動時受的相互作用力；Ii和mi分別為顆粒i的轉動慣量和質量；Vi和ωi分別為顆粒i的速度和角速度；Ft,ij為切向分力；Fn,ij為法向分力；Tr,ij為滾動摩擦力矩；Tt,ij為切向力矩。

根據彈性阻尼的線性三方程函數模型，每個力都以彈簧、阻尼和滑塊的組合形式將其簡化出來表示[14]。

式中，V為顆粒的運動速度；δ為顆粒間的位移變形；a為扭轉變形；k為剛度系數；f為摩擦系數；η為阻尼系數；L 為重疊量。ij表示顆粒i與顆粒j之間，t表示切向，n表示法向，s表示滑動，r表示周向或滾動。

顆粒在氣流中受的相互作用力P表示如下。

式中，ρ為氣流密度；kg為阻力系數，與物料結構形狀、表面特性和雷諾數等有關；A為物料迎風面積；vw為物料速度；v為物料在氣流中的相對速度；vq為氣流速度。

1.3 DEM-CFD 耦合流程介紹及辣椒風選效果評價指標

本文用CFD-DEM 耦合方法，即Fluent 與EDEM的耦合，其流程如圖5所示。

圖5 DEM-CFD耦合流程圖Fig. 5 Flow chart of DEM-CFD coupling

為深入探討風速大小對辣椒風選效果的影響，采用控制變量法保證其他的各項條件參數相同，模擬辣椒風選的風速數值大小分別為15、20和25 m·s-1時的風選過程。

辣椒清選裝置工作性能的主要衡量指標是風選后辣椒損失率（El）和含雜率（Ed），因此以辣椒損失率El、含雜率Ed為辣椒風選裝置的性能評價指標[15-16]。

式中，W1為仿真過程中通過出口Ⅰ的辣椒質量；W2為仿真過程中通過出口Ⅱ的辣椒質量；M為仿真過程中通過出口Ⅰ的物料總質量。

1.4 模擬參數設置

通過WDW3100 微機控制電子萬能試驗機對朝天椒、辣椒葉及辣椒莖稈進行壓破試驗和拉伸試驗求取泊松比與剪切模量[17-18]，密度用天平和量筒法測量；辣椒、辣椒莖稈和辣椒葉與鋼板的摩擦接觸參數采用斜板試驗（分別為滑動法和滾動法）測量；辣椒、辣椒莖稈、辣椒葉與鋼板之間的恢復系數采用碰撞法測量[19-20]。EDEM軟件中材料的力學性能參數和材料顆粒的接觸系數見表1和表2，辣椒清選分離裝置材質為鋼材。模擬中辣椒、辣椒葉、辣椒莖稈的生成比例設置為6∶1∶0.25，辣椒的生成速度設置為1 200個·s-1，辣椒葉的生成速度設置為200 個/s，辣椒莖稈的生成速度設置為50 個·s-1。EDEM中的仿真時間步長設置為29.6784%，即0.000 007 s，總仿真時間為6 s，Fluent 軟件中的仿真模型采用標準的k-ε湍流模型，時間步長設置為EDEM 的100 倍，即時間步長設置為0.000 7 s，總步數為8 571步。

表1 材料的物理參數Table 1 Physical parameters of the material

表2 相關顆粒的材料接觸系數Table 2 Material contact coefficient of related particles

2 結果與分析

2.1 辣椒顆粒運動軌跡分析

以風速大小為20 m·s-1的耦合仿真最終結果為例，辣椒顆粒的運動軌跡如圖6 所示?？梢钥闯觯苯贰⒗苯啡~和辣椒莖稈由顆粒工廠出料生成，然后自由運動下落并進入風選裝置，當生成的顆粒通過氣流通道流體區域時，顆粒質量較輕的辣椒葉顆粒和辣椒莖稈顆粒能隨氣流運動，被吹到出料口Ⅱ、Ⅲ從而下落，較重的辣椒在重力作用下繼續自由下落，最終從出料口Ⅰ中掉出，從而實現辣椒、辣椒葉和辣椒莖稈的分離。從圖6 中顆粒軌跡圖可以直觀地觀察風選過程中顆粒軌跡，辣椒的紅色軌跡線是一條略帶斜度的直線，表明辣椒受風影響較小，水平位移較小，辣椒可以從出口Ⅰ自由落下；辣椒葉的綠色軌跡線是一條轉彎的斜線，表明辣椒葉受風影響較大，辣椒葉隨風揚動，從出口Ⅲ出來；代表辣椒莖稈的藍線軌跡線是斜線，表示辣椒莖稈受風影響較明顯，辣椒莖稈隨風運動，從出口Ⅱ被吹出。

圖6 辣椒顆粒的運動軌跡Fig. 6 Schematic diagram of the winnowing of pepper particles

2.2 風速20 m·s-1時辣椒風選效果分析

如圖7 所示，在EDEM 軟件中設置區域Ⅰ為辣椒含雜率的統計區域，統計0～6 s 內通過該區域的物料總質量M和辣椒質量W1；區域Ⅱ設置為辣椒損失率的統計區域，統計辣椒在0～6 s 內通過該區域的辣椒質量W2。

圖7 物料顆粒位置瞬態圖（t=1 s）Fig. 7 Transient diagram of material particle position （t=1 s）

從辣椒風選含雜率曲線（圖8）可以看出，0～0.5 s 間的含雜率曲線在X軸上，這是因為該時間段為風選顆粒產生及開始下降的風選過程，經過0.5 s 的緩沖，顆粒經過風力風選位置并進行風選，風選后的辣椒顆粒經過辣椒含雜率統計區域并從出口Ⅰ掉落。通過風選模擬仿真分析可知，經過風力風選之后，辣椒顆粒的含雜率降明顯，通過風選裝置辣椒出口Ⅰ的辣椒平均含雜率為4.624%，通過風選裝置其他出口的辣椒平均損失率（圖9）為1.798%。

圖8 風速20 m·s-1的辣椒含雜率Fig. 8 Impurity rate in pepper with wind speed of 20 m·s-1

圖9 風速20 m·s-1的辣椒損失率Fig. 9 Pepper loss rate at wind speed 20 m·s-1

2.3 不同風速對風選效果的影響

從圖10 可以看出，在較低風速（15 m·s-1）下，在辣椒出口處發現了較多雜質。辣椒顆粒分離度低；在20 m·s-1的較高風速下，混入的辣椒雜質很少，辣椒顆粒的分離度較高；而當風速為25 m·s-1時，幾乎沒有雜質混入辣椒出口處，然而，在莖稈出口處大量的辣椒被吹出，辣椒損失率過高。

圖10 不同風速下的辣椒顆粒分布Fig. 10 Distribution of pepper particles at different wind speeds

不同風速下辣椒含雜率曲線及辣椒損失率曲線如圖11、12所示，通過對辣椒在不同風速下的仿真結果分析，可知當風速大小為15 m·s-1時，辣椒含雜率峰值出現在4~5 s 之間，達到22.234%，風選后的辣椒平均含雜率為11.774%，辣椒平均損失率為0.959%；當風速大小為20 m·s-1時，辣椒含雜率峰值出現在5~6 s 之間，達到8.713%，辣椒平均含雜率為4.624%，辣椒平均損失率為1.798%；當風速大小為25 m·s-1時，辣椒含雜率峰值出現在2~3 s之間，達到3.554%，辣椒平均含雜率為1.676%，辣椒平均損失率為6.476%。根據模擬數據分析可知，在不考慮損失的情況下增加風速強度可以增加風選的效果。

圖11 不同風速下的辣椒含雜率Fig. 11 Impurity rate of pepper at different wind speeds

圖12 不同風速下的辣椒損失率Fig. 12 Loss rate of pepper at different wind speeds

不同風速下的辣椒仿真結果如表3 所示。分析發現，當風速從15 m·s-1增大到20 m·s-1時，風選裝置出口Ⅰ處辣椒的平均含雜率下降了7.150%，而平均損失率增加了0.839%；當模擬的風速再次從20 m·s-1增大到25 m·s-1時，辣椒的平均含雜率下降了2.948%，而平均損失率上升到4.678%。這是由于氣流速度增加，即氣流水平力FD增大，故合力F和風力FD之間的角度β及物料顆粒在水平方向上的位移變大，從而導致辣椒損失率變大，因此增加氣流速度可以降低辣椒的雜質含量，同時也會導致辣椒的損失率增加。

表3 不同風速的仿真結果Table 3 Simulation results at different wind speeds

由以上比分析可知，風速的增加可以增強辣椒的風選效果，但是隨著風速的增加，會導致辣椒從辣椒葉出口吹出的風險。綜合考慮辣椒損失率及含雜率，在風選時應選擇分離率較高及損失率較低的風選速度，即在20 m·s-1風速時分離效果最佳，此時，辣椒的含雜率為4.624%、損失率為1.798%。

3 討論

鄒道鐘等[21]設計了一種適用于板椒收獲機的清選系統，系統總體分為上、下清選機構，復脫機構以及風機風選機構，由于復脫機構中定、動彈齒的高速沖擊，更容易損傷板椒果皮，增大了采收時的損傷；趙飛龍等[22]使用液壓動力代替傳統的機械動力來為風扇提供動力，采用PID（proportion integration differentiation）反饋控制速度回路可以快速響應動態負載變化，減少液體對風機回路的影響，可提高風選裝置的使用壽命及提高風選裝置的輸入動力，能有效提高風選效果。

前人對辣椒的風選研究還停留在憑經驗的機械設計及動力分析上，對風選裝置中物料風選過程的運動特性、流場分布特性之間的關系研究較少。近年來，EDEM 與Fluent 軟件耦合法逐漸成為研究流固兩相間流動性的重要研究手段和方法[23-24]，國內外學者使用CFD-DEM 耦合方法研究農作物的風選，并取得了一些成果：Yuan 等[25]通過DEM-CFD 耦合方法研究了谷物和莖稈的運動行為和篩分特性；蔣恩臣等[26]模擬了聯合收割機分離室中谷物的運動并進行了試驗，表明DEMCFD 耦合方法適用于模擬谷物清選過程；王立軍等[27]使用DEM-CFD 耦合方法研究了玉米清選并使玉米粒的損失率得到了降低。為了進一步研究辣椒的清選，本研究使用DEM-CFD耦合方法對辣椒清選裝置進行仿真。通過非結構化的ICEMCFD 網格模型將仿真模型劃分網格并導入EDEM和Fluent 軟件中，并在2 個軟件之間建立耦合關系以進行模擬仿真。Fluent 軟件中流體的耦合模型采用歐拉模型，歐拉耦合模型使用多顆粒相流耦合框架方程來進行求解，有考慮體積能量分數，該方程考慮不同顆粒相間對流場的相互影響[28]，結合通過進入辣椒的表層進行循環的單向氣流，辣椒層的顆粒對循環氣流產生阻礙力的作用明顯。EDEM 軟件中的辣椒混合物顆粒在建模時，考慮到EDEM 軟件本身自動模擬的一些不足，選擇辣椒清選除雜過程中成分含量高的辣椒、辣椒葉和辣椒莖稈作為主要研究對象進行物料粒子建模。由于EDEM 軟件中的所有粒子模型都用球體表示，因此采用“多球體類聚法”對3 種不同的研究粒子進行簡化、重疊和填充，以代替實際風選的粒子模型[29]。本文采用CFD-DEM 耦合法模擬辣椒在不同風速下的清選分離過程，結合空氣動力學分析辣椒、辣椒葉、辣椒莖稈在流場中的受力情況和運動趨勢。上述研究表明，當風選速度不變時，物料顆粒本身的重力越大，其運動的方向角越大，則物料顆粒的水平位移就越大；此外，也驗證了DEM-CFD 耦合方法對于研究辣椒清選分離是有效的，得到的清選分離參數可以為辣椒清選分離裝置的設計和優化提供指導。