基于CFD-DEM的飼料調質器物料運動模擬與試驗

彭 飛 方 芳 王紅英 黃志剛 付宗強 高東明

(1.北京工商大學材料與機械工程學院， 北京 100048; 2.鄭州大學化工與能源學院， 鄭州 450001;3.中國農業大學工學院， 北京 100083)

0 引言

飼料作為動物主要的食物來源，是畜禽和水產養殖業的物質基礎，2017年全國飼料總產量達到2.24億t[1]。調質是飼料加工過程中非常關鍵的環節，是影響顆粒飼料質量的重要因素[2]。但由于調質器作業過程的封閉性、濕黏飼料原料物性的復雜性，目前國內外對調質過程的機理和調質器作業的性能研究較少，難以進行精確的定量分析[3]。

離散元法(Discrete element method，DEM)是分析顆粒離散體物料的一種方法，1971年由美國Cundall教授基于分子動力學原理提出，被大量應用在復雜物理場作用下粉體動力學現象、具備較為復雜結構的材料力學特性和多相混合材料介質等研究中，涉及粉體生產加工、研磨和混合攪拌等生產實踐領域[4-6]。計算流體動力學(Computational fluid dynamics，CFD)能夠很好地模擬流體在腔體內的氣壓分布和特定層面的速度分布，目前針對調質器工作過程的仿真模擬還比較少，前人做過類似調質器結構的建模研究[7-10]。DEM能夠構建黏性顆粒模型，還可以與CFD耦合計算。李洪昌等[11]通過CFD-DEM耦合方法模擬了風篩式清選裝置中物料在篩面上的運動，其中物料使用EDEM建模，連續流體使用CFD建模，研究并分析入口氣流速度對該裝置篩分性能的影響。劉佳等[12]基于CFD-DEM耦合方法，模擬了機械氣力組合式精密排種器的工作過程，模型采用EDEM軟件中的bonding和API替換的方法，建立了非球形虛擬玉米顆粒。彭飛[13]初步構建了調質器CFD-DEM耦合模型，但是該模型僅局限于定性描述，缺乏試驗數據和定量分析等，模型的可靠性與合理性尚未驗證，存在一定的局限和不足。以上文獻為深入分析與研究調質器CFD-DEM耦合模型提供了理論基礎和方法參考。

本文應用CFD-DEM耦合計算方法，分析飼料原料在調質器內的運動和黏結情況。基于構建的調質器耦合模型，通過三因素五水平正交回歸模擬試驗，分析槳葉安裝角、調質軸轉速、填充率3個因素對其出料量的影響，并利用試驗裝置進行試驗驗證，為調質器結構設計和作業參數優化提供參考。

1 調質器結構和工作原理

如圖1所示，調質器部件主要由進料單元和調質單元組成，具體包括：進料口、進料螺旋、蒸汽腔、蒸汽進口、調質器外壁、主軸、槳葉、調質腔、出料口及聯軸器等[14-15]。

圖1 飼料調質器結構簡圖Fig.1 Structure diagram of feed conditioner1.進料口 2.進料螺旋 3.蒸汽腔 4.蒸汽進口 5.調質器外壁 6.調質器主軸 7.槳葉 8.調質腔 9.出料口

蒸汽進口位于調質單元物料進入處，并與調質腔體連通。在調質器驅動電動機和變頻器的帶動下，進口螺旋葉片推動原料至調質腔內；蒸汽發生器產生蒸汽，經蒸汽腔、環形分布的蒸汽加工孔，進入到調質腔內；在調質腔內旋轉扇形槳葉的攪拌作用下，飼料原料和蒸汽受到擠壓、剪切、翻滾和拋出等強制混合作用，在劇烈的相對運動中均勻混合并產生水熱反應，由出料口流出，完成調質過程。

2 調質器CFD-DEM數值模擬

2.1 數學模型

利用CFD-DEM耦合模擬時，耦合模型主要分為Lagrangian模型和Eulerian模型[16]。其中，Eulerian模型考慮了顆粒對流場的影響(適用于顆粒相占總體積大于10%的情況)[17]，根據飼料調質器作業實際工況，原料試驗顆粒對流體相的作用明顯，因此本文采用Eulerian模型對飼料原料調質過程進行模擬。考慮到原料對流場的影響，在守恒方程中加入一個額外的體積分數，即

(1)

式中ε——氣體體積分數

ρ——氣體密度，kg/m3

t——時間，s

u——氣體速度，m/s

動量守恒方程為

(2)

其中

式中g——重力加速度，m/s2

μ——氣體動力黏度，Pa·s

S——動量匯，即作用在網格單元內流體阻力的總和，N/m3

Fi——第i顆粒對氣流的阻力，N

V——網格單元的體積，m3

由于飽和蒸汽對飼料原料的影響主要為曳力，故選擇Free-stream阻力模型，該自由流阻力的計算公式為

Fd=0.5CDρA|v|v

(3)

式中CD——單個顆粒單元的曳力系數

A——顆粒投影面積，m2

v——氣體相對顆粒單元的流動速度，m/s

飼料顆粒的曳力系數CD取決于雷諾數Re，其計算公式為

(4)

其中

式中L——原料顆粒直徑，m

α——CFD網格單元的自由體積，m3

針對仿真對象的不同，選擇相應的接觸模型[17]。由于飼料原料進入調質器后處于濕黏狀態，為典型的黏性散粒體，需要引入黏結鍵來表征顆粒彼此接觸時的運動狀態。Hertz-Mindlin計算顆粒間法向和切向黏結力分別定義為Fn、Ft，法向和切向黏結力矩分別為Tn、Tt，單個顆粒法向和切向的速度分別為vn、vt，單個顆粒法向和切向的角速度分別為ωn、ωt。黏結力Fn、Ft和力矩Tn、Tt隨著時間步長的增加，按公式從零開始增加。

(5)

式中RB——黏結半徑，m

δt——時間步長，s

Sn——顆粒法向剛度，N/m

St——顆粒切向剛度，N/m

A1——顆粒間的接觸面積，m2

當法向和切向應力超過某個定義值時，黏結就會被破壞。因此，定義法向應力最大值σmax和切向應力最大值τmax為

(6)

2.2 調質器建模與仿真

構建調質器內腔幾何模型，采用Fluent 12.0對調質器內腔模型進行參數設置，采用Fluent 12.0中的DEM模塊對散粒體參數進行設定，DEM由軟件EDEM 2.2求解，兩者耦合安裝并計算，耦合方法的基本思路：利用CFD方法計算壓力場及流場分布，利用DEM方法計算顆粒系統的受力和運動情況；兩者基于一定的模型進行質量、動量和能量等物理量的作用與交換，實現耦合計算。

2.2.1結構模型構建和參數設置

基于調質器尺寸，構建調質器Pro/E模型，步驟如下：

(1)調質腔模型主要幾何參數：圓柱腔體直徑為100 mm，長度為600 mm；空心主軸位于腔體中心，直徑30 mm，長度與調質腔相同；扇形槳葉共18組，距離調質腔內壁2 mm，按軸向間隔30 mm、徑向間隔90°分布，安裝角度可調；蒸汽進口為10組，直徑為6 mm，呈環形均勻分布于調質腔進料口處的端蓋上；出料口位于調質腔右端40 mm，其結構為Φ24 mm的圓柱體。由于裝配體零件較多，模型導入到Gambit軟件時，導致Gambit模型中線、面、體較多，達到上千個，不便于下一步在Fluent軟件中設置邊界類型和邊界條件以及后續的模擬計算，因此需要對模型進行簡化。使用Pro/E 4.0建立零件1(調質器內腔)、零件2(調質軸與槳葉)這兩個零件(圖2、3)，然后進行裝配。

圖2 調質器內腔幾何模型Fig.2 Geometric model of conditioner inner cavity

圖3 調質軸與槳葉幾何模型Fig.3 Geometric model of conditioner shaft and blade

圖4 Gambit中調質器模型Fig.4 Model of conditioner in software Gambit

圖5 利用Gambit對模型進行邊界條件設置和網格劃分的結果Fig.5 Boundary condition setting and meshing of model by using Gambit

(2)將裝配體保存為igs格式并導入到Gambit 2.4.6中(圖4)，得到由2個獨立構件組成的幾何體，共計226個面；結合實際作業情況和邊界類型設定方法，將這些面定義為4種類型：進氣口(inlet)、出氣口(outlet)、壁面(wall)、旋轉部分(wall)。對圖4所示模型進行體網格劃分，網格尺寸為8 mm，零件1得到219 205個網格，零件2得到21 938個網格，結果如圖5所示；最后將劃分完成后的裝配體以mesh格式保存，作為Fluent的mesh文件。

(3)將mesh文件導入到Fluent，流體模型采用RNGk-ε黏性模型；設置進氣口為壓力入口，出氣口為壓力出口；基于物料所占腔體體積比例，選擇Eulerian法進行耦合，耦合參數匹配設置方法參考文獻[18-19]。調質過程中，原料與飽和水蒸氣發生傳熱傳質反應，原料含水率由10%～12%增加到17%～19%，物料黏度增加；基于容重法國標[20]，測得該含水率范圍條件下，飼料物料密度平均值為443 kg/m3；在EDEM中將物料間離散元參數設置為黏結模型(Hertz-Mindlin with bonding)。EDEM模型中材料參數設置如表1所示。

表1 EDEM模型中材料參數設置Tab.1 Material parameter setting of EDEM model

(4)其它參數設置：設置離散元模型中時間步長為2×105s，為瑞利時間步長的14%，該值在瑞利時間步長比例的建議范圍之內，保存時間間隔為0.01 s。設置顆粒生成量為unlimited number。由于該調質器調質時間為15～25 s，為便于模擬調質器穩定工作前及穩定工作一段時間后的作業狀況，仿真時間應適當大于實際生產達到穩定的調質時間，因此在EDEM中將仿真時間設置為30 s。

2.2.2調質過程仿真分析

仿真預試驗時，試驗因素均設定為零水平，分別是：槳葉安裝角為45°、調質軸轉速為250 r/min、填充率為35%。圖6為調質器仿真模擬30 s過程中迭代收斂曲線，由圖分析可知，仿真收斂性較好，在0～30 s內調質器內部的流場速度趨于平穩。在EDEM中統計并分析調質器內顆粒數量隨時間變化情況，結果如圖7所示，可以看出，大概在25 s以后，整個腔體內的顆粒數目趨向穩定、基本不再增加或減少，即進料口和出料口顆粒數量基本維持平衡，此時調質器模型趨于穩定狀態，可以用于模擬調質器實際生產中達到穩定的過程。

圖6 調質器模擬工作30 s過程中迭代收斂曲線Fig.6 Iterative convergence curves of simulation during 30 s

圖7 0～30 s過程中調質器內顆粒數量變化曲線Fig.7 Changing curve of particle number in conditioner during 0～30 s

2.3 仿真結果與分析

圖9 0～25 s調質器穩定工作前顆粒分布和運動狀態Fig.9 Particle distribution and movement during 0～25 s

圖8a、8b分別為30 s時調質器速度場與壓力場分布圖，由垂直的兩個截面顯示。由圖可以看出，10個蒸汽進口處氣體流速較大；隨著不斷遠離蒸汽進口，氣體流速逐漸減慢；沿著物料運動方向，速度場與壓力場整體變化均勻。在調質器作業過程中，蒸汽與原料顆粒充分接觸并發生水熱反應，這與該調質器結構和調質工作原理基本一致。

圖8 30 s時調質器速度場與壓力場分布圖Fig.8 Velocity and pressure field in conditioner at 30 s

在DEM模型中分別統計0、5、10、15、20、25 s情況下調質器內顆粒的運動和分布情況，結果如圖9所示。由圖可知，在該時間段內，顆粒一直不斷生成，調質器進料速率大于出料速率，因此調質腔內顆粒物料不斷增多，腔體物料填充率不斷提高；同時調質軸旋轉并帶動槳葉攪拌顆粒物料，槳葉有一定的安裝角度，旋轉的槳葉起到攪拌、拋起并推動物料前進的作用，在25 s左右進出調質器腔的顆粒數量基本持平，此時填充率基本不再變化，這與調質器開始作業階段的真實狀態基本一致。圖10為調質器穩定工作過程階段，此階段腔體內顆粒數目和填充度變化不大，基本處于穩定的狀態。由圖10f可知，顆粒黏結在一起呈團簇狀，調質反應后的物料也呈團簇狀被翻起攪拌，與實際生產情況相符合，可知選用的顆粒黏結模型基本合理；由顆粒速度分布規律可知，顆粒受到槳葉的攪拌作用，不斷被拋起并向前推進，在調質腔中部及上部、靠近槳葉末端處速度較快；靠近調質器底部顆粒的運動速度較慢，有可能導致物料黏附在該處內壁上、產生積料現象，該模擬過程基本與調質器作業過程一致。物料速度場與壓力場分布均勻，顆粒前進方向與蒸汽進氣方向一致、兩者接觸較充分，該模型整體作業過程較為合理，可以用于模擬調質器作業時飼料原料的分布和運動情況。

圖10 26～30 s調質器穩定工作后顆粒分布和運動狀態Fig.10 Particle distribution and movement during 26～30 s

3 試驗設計與指標測定

3.1 試驗設計

由預試驗及文獻[21-22]可知，調質器工作性能的主要影響指標為槳葉安裝角、調質軸轉速、填充率，故本文選擇上述3個因素進行仿真試驗研究。根據理論分析、單因素預試驗結果、調質器結構參數與生產實際，確定各試驗因素的取值范圍為: 槳葉安裝角10°～80°、調質軸轉速150～350 r/min、填充率10%～60%。以槳葉安裝角X1、調質軸轉速X2、填充率X3為試驗變量，基于二次正交旋轉組合試驗原理，建立因素編碼表如表2(x1、x2、x3為各變量真實值)所示。

表2 二次回歸正交試驗設計因素編碼Tab.2 Factors and levels of quadratic regression orthogonal rotating experiment design

3.2 指標測定

調質器在調質作業過程中的工作效率可以通過出料口處的出料量來檢測。通過EDEM軟件的數據導出工具Export Results Data，導出數據并計算出仿真過程中調質器的出料量。

4 試驗結果與分析

4.1 回歸模型建立

以各影響因素編碼值為自變量，以仿真結果測得的出料量Y為評價指標，依據二次回歸正交旋轉組合設置不同試驗組的模型參數，基于CFD-DEM耦合仿真試驗，結果如表3所示(X1、X2、X3為各變量編碼值)。

表3 二次回歸正交旋轉組合試驗設計與結果Tab.3 Experimental design and results of quadratic regression orthogonal rotating test

圖11 基于響應面法的參數組尋優結果Fig.11 Optimization of parameter group by using response surface method

變異來源平方和自由度均方FP模型129.13914.3516.80<0.0001X19.3419.3410.930.0057X220.76120.7624.310.0003X367.70167.7079.27<0.0001X1X20.6310.630.730.4070X1X31.6911.691.980.1826X2X301000.9701X2124.95124.9529.210.0001X220.5210.520.600.4513X233.7513.754.390.0563

采用Design-Expert軟件對試驗進行回歸分析，得到槳葉安裝角X1、調質軸轉速X2、填充率X3與調質器出料量Y的回歸方程

(7)

4.2 回歸模型尋優

為直觀地分析試驗因素與評價指標之間的關系，根據各回歸項對評價指標的影響，結合簡化回歸模型，用“降維法”將任意1個因素固定中心水平[23-24]，得到另外2個試驗因素與評價指標之間的降維回歸模型，由Design-Expert軟件繪制響應面圖，分析各因素對出料量的影響，結果如圖11所示。

在試驗范圍內，將填充率X3固定在中心水平上，得到槳葉安裝角X1和調質軸轉速X2的交互作用對調質器出料量的影響。由圖11a可知，槳葉安裝角和調質軸轉速對出料量的影響為上凸型曲面；當調質軸轉速一定時，出料量隨槳葉安裝角的增加呈現先增大后減小的趨勢，這是因為物料在調質器中主要有軸向和徑向兩種運動形式，當槳葉安裝角較小時，槳葉主要帶動物料圍繞調質軸轉動，軸向推動力較小，因此物料軸向向前運動較弱，出料量較小；隨著槳葉安裝角的增大，物料軸向推進作用增強，出料量隨之增加；隨著槳葉安裝角的進一步增大，槳葉對物料軸向推進作用減弱，主要起到切割物料作用，既不能推進也不能拋起物料。X1和X2的交互作用不顯著。

在試驗范圍內，將調質軸轉速X2固定在中心水平上，得到槳葉安裝角X1和填充率X3的交互作用對調質器出料量的影響。由圖11b可知，槳葉安裝角和填充率對出料量的影響為上凸型曲面；當槳葉安裝角一定時，出料量隨填充率的增加呈現增大的趨勢，且變化幅度隨安裝角的增大呈現先逐漸加劇、后逐漸減緩的規律。根據仿真情況及調質器輸送機理分析，這可能是因為原料靠近槳葉部分出現的翻滾現象造成的。當填充率較低時，大量物料集中于遠離調質軸的區域，徑向速度低而軸向速度高；隨著填充率的增大，更多的原料堆積使得靠近調質軸的原料增多，物料軸向速度變慢，因此出料量增加，但是增長速度變小。楊洋[25]基于離散元法和正交試驗研究了糧食輸送器的輸送效率，羅勝等[26]研究了螺旋不連續加料裝置的出料量變化，有相似的規律和結論。X1和X3的交互作用不顯著。

在試驗范圍內，將槳葉安裝角X1固定在中心水平上，得到調質軸轉速X2和填充率X3的交互作用對調質器出料量的影響。由圖11c可知，調質軸轉速和填充率對出料量的影響為上凸型曲面；調質軸轉速一定時，出料量隨填充率的增加呈現增大的趨勢；填充率一定時，出料量隨槳葉安裝角的增加呈現先增大后減小的趨勢。X2和X3的交互作用不顯著。

4.3 最優參數組合的確定及驗證

由于該調質器最大生產率為50 kg/h，取最大生產率的90%～95%為宜，故Y在該范圍內為宜。基于響應面法，運用Design-Expert軟件中Optimization模塊進行參數優化，在-1.682≤Xi≤1.682 (i=1，2，3)范圍內對各參數進行進一步尋優，獲得調質器最優作業參數組合為：X1=-0.31、X2=-1.13、X3=1.56，即x1=38.49°、x2=182.2 r/min、x3=58.4%，此時調質器出料量Y可獲得最優目標值，仿真預測值為13.1 g/s。

5 樣機試驗

取最優參數組合，在北京市密云區昕三峰飼料廠進行調質器部件車間試驗，如圖12所示；加工對象為正常生產的乳豬料原料，物料理化指標良好。調整設備至最優參數組合，待調質器工作穩定后，測得調質器出料量實測值為12.8 g/s，其它指標性能良好，測得效果基本與優化試驗結果一致，滿足調質器實際生產需求。

圖12 調質器生產試驗Fig.12 Prototype working in workshop

6 結論

(1)基于離散單元法建立了飼料原料黏結顆粒模型，運用CFD-DEM耦合計算方法，對調質器原料顆粒進行運動和分布、顆粒黏結情況的耦合仿真。分別統計0～30 s階段時調質器內顆粒的運動和分布情況，由仿真分析可知，原料受槳葉的攪拌作用，不斷被拋起并向前推進，在調質腔中部及上部、靠近槳葉末端處速度較快；靠近調質器底部顆粒的運動速度較慢，有可能導致物料黏附在該處內壁上、產生積料現象。

(2)基于CFD-DEM耦合模型，通過三因素五水平正交組合試驗，得出各因素對其出料量的影響顯著性順序依次為：填充率、調質軸轉速、槳葉安裝角。通過Design-Expert軟件對試驗結果進行回歸分析和響應面分析，得到調質器最優作業參數組合為：槳葉安裝角38.49°、調質軸轉速182.2 r/min、填充率58.4%，此時出料量仿真值為13.1 g/s，試驗值為12.8 g/s。通過對比仿真值和試驗值，驗證了仿真試驗與回歸模型的有效性。