孫曉霞， 孟文俊， 牛雪梅， 袁媛

(1.太原科技大學 機械工程學院，山西,太原 030024; 2.山西能源學院，山西,太原 030024)

隨著垂直螺旋輸送機的廣泛應用，作為其重要組成部分的取料裝置也有了各種形式，其中最重要的一種形式就是相對旋轉式喂料頭(簡稱喂料頭)，其喂料量和喂料效率關系到輸送機的輸送效率[1-2]. 在國際上，垂直螺旋輸送機及其喂料頭的技術相對先進，提出了如多片導料擋板式反向旋轉喂料頭、雙口螺線式喂料頭、螺旋擋板式喂料頭、傾斜翼板式喂料頭等多種結構的喂料頭[3].對喂料頭的研制是通過試驗，測得不同尺寸的喂料頭模型在不同的喂料頭轉速的情況下的充填率、功率及生產率，由試驗曲線描述取料段性能指標隨喂料頭結構參數的變化規(guī)律，確定較佳的設計變量. 相對于國外先進的技術優(yōu)勢，喂料頭在我國目前尚未有統(tǒng)一的標準，每個廠家都有各自的加工標準和輸送量確定的標準.

喂料頭的喂料效率和喂料量直接影響到垂直螺旋輸送機的輸送效率，而對喂料頭的喂料量產生影響的主要因素就是喂料頭的喂料曲面. 因此，喂料曲面的投影、展開平面以及特征曲線的形狀及方程也是各有特點. 顆粒進入喂料頭后，由于料壁內螺旋的高速旋轉，會產生氣流影響顆粒的運動，所以對于喂料頭中的顆粒，不僅要考慮顆粒與顆粒、喂料頭之間的相互作用，還要考慮顆粒與氣體之間的相互作用. 顆粒在喂料頭中的運動比較復雜，受力也比較復雜，難以在理論上深入探討顆粒在喂料頭中的運動過程[4].

近年來，隨著計算機技術的快速發(fā)展，數值模擬已逐漸成為與試驗手段相輔相成的研究工具. 宋斌等[5-6]采用DEM模型及EDEM軟件仿真對喂料頭進行了研究，分析了垂直螺旋輸送機內顆粒在不同條件下的運動情況，得出影響顆粒運動的因素. 但是喂料頭中需要考慮氣體對顆粒的作用，所以需考慮一種數值模型，不僅需要考慮顆粒與顆粒、喂料頭之間的相互作用，還要考慮顆粒與氣體之間的相互作用[7]. 文中采用E/L法對垂直螺旋輸送機進行模擬，由于顆粒密度較高，顆粒之間的碰撞不能忽略，而離散單元法(DEM)是顆粒碰撞模型中最應用廣泛的一種. 對于模擬喂料頭這種復雜的顆粒流系統(tǒng)，建立氣固兩相DEM模型，采用離散單元法和計算流體力學(CFD)相結合的方式，可以跟蹤每個顆粒的運動，獲得大量的微觀信息. 在DEM模型及其仿真的方法基礎上，結合理論分析與氣固兩相DEM模型，及DEM+CFD耦合仿真分析顆粒的喂料過程，考察顆粒在不同曲面喂料頭中的運動狀態(tài)及輸送量情況.

1 氣固兩相DEM模型及其仿真方法

在喂料頭中，因為顆粒所受外力及摩擦力的復雜性和不確定性，難以用傳統(tǒng)的牛頓第二定律的方法來確定顆粒的速度，而且，顆粒還受到空氣的影響，因此，下面采用氣固兩相DEM模型及其仿真分析方法來觀察顆粒的輸送量情況.

1.1 氣固兩相DEM模型

1.1.1 顆粒固相控制方程

在氣固兩相DEM模型中，顆粒的運動通過求解牛頓第二定律獲得，

ma=Fg+Fp-p+Fp-g+Fge,

(1)

式中:a為固相顆粒對于喂料葉片的加速度;Fg為重力;方向垂直向下;Fp-p為顆粒與顆粒、喂料葉片之間的作用力，方向不確定;Fp-g為氣體與顆粒之間的作用力，方向不確定;Fge為離心力.

1.1.2 氣相控制方程

在氣固兩相DEM模型中，氣體的運動規(guī)律由N-S方程描述，設εg為氣體空隙率，vg為氣體表觀速度，則氣相控制方程為

(2)

(3)

其中，氣相控制方程中動量方程的流固耦合項為

Fp=β(vg-vp),

(4)

式中β為氣固相間曳力系數，其模型采用Gidaspow曳力模型. 曳力是氣固相間最主要的作用力.可以看出，氣體對顆粒的曳力Fd與顆粒對氣體的力Fp是一對相互作用力.

1.2 接觸曳力模型

文中采用DEM模型及支持此模型的仿真軟件EDEM軟件進行仿真，其提供的顆粒之間的接觸模型有Hertz-Mindlin(no slip)模型等[8]. 文中仿真的顆粒都是無黏結性小顆粒，根據實際情況選擇Hertz-Mindlin 無滑動接觸模型，該模型能高效計算兩個顆粒接觸時的接觸作用力，由法向力和切向力以及滑動摩擦力和滾動摩擦力組成. 法向力分量由Hertzian的接觸理論計算，切向力由Middlin-Deresiewicz計算，它們都有阻尼分量，其阻尼系數和恢復系數有關. 庫倫摩擦定律決定滑動摩擦力，接觸獨立定向恒轉矩模型決定滾動摩擦力.

在顆粒的運動過程中，顆粒受到來自氣體的曳力，而氣體又因為顆粒的存在對其流動產生影響，因此兩者之間存在緊密的耦合作用，曳力影響了顆粒的輸送，因此曳力模型的準確性是非常關鍵的.

1.3 氣固兩相耦合

氣固兩相DEM模型模擬程序由空氣流場Fluent計算和顆粒動力學EDEM計算兩部分構成，模擬時，顆粒場與空氣流場的求解交叉進行. EDEM+Fluent耦合求解過程如圖1所示. 在一個時間步長內，Fluent通過UDF和EDEM中的API接口相結合，先在Fluent中進行流場計算，當迭代收斂時，觸動EDEM開始當前時步的顆粒計算，并把Fluent中的流場數據傳遞給EDEM，根據Gidaspow曳力模型計算氣體與顆粒間的相互作用力，再根據Hertz-Mindlin(no slip)接觸模型用以計算顆粒與顆粒、喂料葉片、螺旋葉片之間的相互作用力，然后，根據牛頓動力學方程給顆粒定位. 當以上步驟在EDEM中完成后，再把氣體與顆粒間的作用力反饋回Fluent中，開始新一代的迭代.

圖1 EDEM+Fluent軟件耦合求解過程Fig.1 Coupling and solving process of EDEM with Fluent

2 喂料葉片曲面

2.1 曲面曲線方程

當顆粒進入喂料頭處的窗口時，隨著螺旋徑向尺寸的逐漸減小，迫使顆粒進入螺旋輸送管底部中心附近或灑落在螺旋葉片上，螺旋軸的高速旋轉使顆粒向上輸送，螺旋擋板沿高度的傾斜有效地防止了顆粒離心作用引起的濺出. 圖2是文中研究的喂料頭結構圖.

圖2 喂料頭結構Fig.2 Structure of feeder head

喂料頭的輸送量由進入的顆粒和流出的顆粒之差組成，顆粒進入的越多，流出的越少，輸送量越多，輸送量的多少是由螺旋的轉速、喂料頭的轉速和喂料頭的結構決定的，而喂料頭的結構主要考慮喂料導擋葉的曲面尺寸[8].

利用反向工程，得到某種高效的喂料頭的表面數據，并利用這些數據得到了實物原型的三維模型以及喂料曲面的展開圖，文中的喂料頭的數據都是源于宋斌等[5]、翟曉晨[6]基于反向工程的數據獲得的.

喂料頭單個導擋葉曲面俯視簡圖如圖3所示.

圖3 喂料頭單個導擋葉曲面俯視圖Fig.3 Sketch for top-view of single feeder surface of guide blade

要擬合的曲線由俯視狀態(tài)下的點云數據擬合獲得，通過測量、分析，在曲線上與中心O的連線達到最大值的點處將邊緣曲線分成q1、q2兩部分，曲線q1的r值變化較大，曲線q2的r值變化較小，分別對這兩部分曲線進行擬合.

以夾角θ為自變量，r為相應的因變量，利用Matlab軟件對能夠反映導擋葉喂料曲面的特征曲線進行了方程擬合，得到的曲線為一個阿基米德螺旋線，設該曲線為曲線3，曲線方程為r=a+bθ，其中r為極徑，θ為極角(弧度).

經擬合后q1的擬合曲線為

r=0.358+0.106θ， 0≤θ<π/3,

(5)

q2的擬合曲線為

r=0.478+0.025θ， 0≤θ≤π/2.

(6)

以上是在喂料口半徑為0.358 m時的擬合曲線，其喂料面積和該曲線是成正比的，所以喂料口半徑也是和該曲線是成正比的，設置喂料口半徑為Rf，則

q1的擬合曲線為

r=Rf+0.3Rfθ， 0≤θ<π/3,

(7)

q2的擬合曲線為

r=1.34Rf+0.07Rfθ， 0≤θ≤π/2.

(8)

當顆粒在喂料葉片中運動時，其有效速度是朝向喂料頭中心的速度vn1，它的朝向中心的受力和運動速度由喂料葉片的阿基米德螺旋線決定，而且喂料頭的轉速不變，顆粒朝向中心的有效速度就不會改變. 根據關系式

t=θ/ωf,

(9)

式中:t為喂料頭旋轉的時間;θ為喂料頭轉過的角度(弧度);ωf為喂料頭角速度.

顆粒流向喂料頭中心的速度為

vnl=(r-a)/t.

(10)

將式(9)(10)聯立，并將其簡化得

vnl=bωf,

(11)

式中:a，b為阿基米德螺旋線的系數，a為θ=0時的極徑;b為每旋轉1°時極徑的增加量. 從式(11)可以看出，流向中心的有效速度vn1只與螺旋線的參數b和喂料頭的轉速ωf有關.

下面改變阿基米德螺旋線b的值，觀察不同導擋葉曲線的喂料狀態(tài). 以漸開線在曲線3的基礎上逐漸增大，b分別增大10和20 mm，得到曲線2和曲線1，以漸開線在曲線3的基礎上逐漸減小，b分別減小10和20 mm，得到曲線4和曲線5. 即漸開線從大到小分別為曲線1～5. 圖4分別為各曲線的q1的曲線圖和q2的曲線圖.

圖4 喂料頭導擋葉片各曲線的曲線圖Fig.4 Curve graph of all kinds of guide blade of feeder head

因此，導擋葉俯視圖的各曲線q1、q2的方程分別為

曲線1:

(12)

曲線2:

(13)

曲線3:

(14)

曲線4:

(15)

曲線5:

(16)

喂料頭單個導擋葉左視如圖5所示.

圖5 單個導擋葉的左視圖Fig.5 Left optic of single feeder surface of guide blade

曲線s和曲線m分別為喂料頭導擋葉的上下曲線，根據反向工程得到曲線s和曲線m的數據，并進行擬合.

可得s曲線擬合函數為

s=0.000 47θ+0.159.

(17)

式中:θ為喂料頭轉過的角度;s為曲線s離底端的距離.

可得m曲線擬合函數為

m=-0.000 007 6θ2+0.002 29θ+0.007 35.

(18)

設導擋葉的曲線m和曲線s不變，分別與前面q1、q2的曲線1～5的5個曲線組成不同的5個曲面，下面分析這5個曲面下的導擋葉的喂料頭的顆粒情況.

2.2 模型建立

由于喂料頭里面的顆粒還受到空氣的影響，采用氣固兩相DEM模型及EDEM+Fluent耦合仿真方法分析喂料頭內的顆粒狀態(tài). 建立模型，取螺旋管徑為0.315 m，喂料口直徑為0.34 m，喂料窗口個數3，將上述在Solidwork中建立的模型導入到EDEM軟件中，然后再結合Fluent軟件進行EDEM+Fluent耦合仿真分析，得到在考慮氣體時喂料頭內的顆粒狀態(tài).

不同時刻的喂料情況如圖6所示.

圖6 喂料頭的喂料過程Fig.6 Feeding states of feeder head

從圖6中可以看出，當運行到第2.0 s時，喂料頭全部埋入顆粒堆中，顆粒開始沿著螺旋向上輸送，此后開始穩(wěn)定向上輸送顆粒. 顆粒在喂料頭中的運動軌跡如圖7所示.

由圖7可以看出，顆粒在喂料頭中的運動非常復雜，顆粒的運動速度方向有朝向中心，進入喂料頭的；也有離開中心，從喂料頭中濺出的. 所以喂料頭的輸送量由兩部分組成，即進入喂料頭的顆粒量，濺出喂料頭的顆粒量.

圖7 顆粒運動軌跡圖Fig.7 Movement trace of bulk materials

顆粒與喂料頭之間的靜摩擦因數是指顆粒與喂料頭之間摩擦力與法向正壓力的比值，它不僅與喂料頭的材料有關，而且與顆粒的表面形狀和粗糙度有關，該摩擦因數對喂料頭中顆粒的運動狀態(tài)起很大的作用，所以文中采用不同的靜摩擦因數μ進行分析，各種仿真工況如表1所示.

表1 喂料頭數值模擬參數表Tab.1 Various simulation working condition of feeding head

如表1所示，分析6種不同靜摩擦因數的顆粒的喂料狀態(tài)，喂料頭轉速ω分別為10，20，30，40，50，60，70，80，90 r/min，顆粒半徑r=2 mm，堆積密度ρ=0.7 t/m3，喂料頭轉速ωf為0～90 r/min. 螺旋管壁里面的螺旋轉速根據文獻[12]中的最佳螺旋轉速定義得出.

2.3 最佳曲面分析

利用氣固兩相DEM模型及EDEM+Fluent耦合仿真分析，分別考察以上5種曲面，并比較不同摩擦因數下不同曲面的輸送量的對比情況.

當在不同的顆粒與喂料頭的摩擦因數下，各曲面在不同轉速下的輸送量分別如圖8所示.

由圖8中可以看出，隨著顆粒與喂料頭摩擦因數的增大，推薦使用漸開線更大些的曲面，因為摩擦因數增大后，顆粒受到喂料葉片的摩擦力越大，顆粒從喂料口濺出量減小，所以曲面漸開線大些好，進料量大，但是曲面漸開線不是越大越好，太大了也容易導致顆粒從喂料口濺出，最大漸開線的曲面1的輸送量并不比曲面2大，在任何摩擦因數下都不建議使用曲面1. 同理，隨著顆粒與喂料頭摩擦因數的減小，推薦使用漸開線更小些的曲面，因為摩擦因數減少后，顆粒受到喂料葉片的摩擦力越小，顆粒從喂料口的濺出量增大，所以曲面漸開線小些好，但是曲面漸開線不是越小越好，太小了顆粒進入喂料口的量減小，所以最小漸開線的曲面5的輸送量并不比曲面4大，在任何摩擦因數下都不建議使用曲面5. 不同摩擦因數下推薦使用不同的導擋葉曲面. 當μ≤0.2時,推薦使用曲面4, 當0.2<μ≤0.5時,推薦使用曲面3; 當μ>0.5時,推薦使用曲面2.

圖8 不同摩擦因數下的5種曲面的輸送量Fig.8 Throughput of five kinds of curve surface in different friction coefficient

3 不同曲線過料面積及輸送量計算

3.1 過料面積

由前面顆粒進料的動力學分析，可知

Qt=Qin-Qout,

(19)

式中:Qt為喂料頭的輸送量;Qin為喂料量，即通過喂料窗戶進入喂料頭里面的顆粒量;Qout為濺出量，即隨著喂料頭的轉動，從喂料窗口濺出喂料頭的顆粒輸送量，它與導擋葉的漸開度和顆粒與喂料葉片的摩擦因數有關.

喂料量Qin表示為

Qin=60V,

(20)

式中V為每分鐘喂料的體積.

V=2πRfDfnωf，

(21)

式中:Df為喂料口面積;Rf為喂料口半徑;n為喂料進口個數;ωf為喂料頭轉速.

由此可得，喂料頭的喂料量為

Qin=120πnDfRfωf，

(22)

單個喂料口的面積Df如圖9中的ABCD所示，垂直輸送管的結構中心線為DE，q1曲線的終點A點與中心E的連線為AE，連線AD與喂料頭的m曲線組成一個截面，該截面被定義為喂料頭上的過料面，可近似地認為是一個梯形.

其中，DE為喂料口半徑，AE為角度為65°時q1曲線的值，DE和AE之間的角度相差45°，AD的值可由AE和DE求出.BC線段為曲線m上兩點之間的線段，B點為角度為60°時曲線m的值，即為梯形中AB的值，C點為角度為120°時曲線m的值，即為梯形中CD的值. 前面已知q曲線、m曲線的方程，從而可求出梯形ABCD中各線段的值，并求出梯形ABCD的面積，即單個喂料面的面積.

圖9 過料面積Fig.9 Single material area

由q1、q2曲線方程(12)～(16)、s曲線方程(17)、m曲線方程(18)，可得出不同的5種曲線的過料面積，如表2所示.

表2 各曲線的過料面積Tab.2 Material area of all kinds of curve

3.2 輸送量公式

設Qt=kQin，則分別計算不同摩擦因數下的喂料量Qin，并得出不同摩擦因數下的輸送量仿真值Qt.

由前面分析可得，當摩擦因數為0.1，0.2時，選擇曲面4，并由表3得出喂料面積為0.018 078 m2，當摩擦因數為0.3，0.4，0.5時，選擇曲面3，并得出喂料面積為0.018 257 m2，當摩擦因數為0.6時，選擇曲面2，并得出喂料面積為0.018 648 m2，然后由式(22)，可計算出各喂料頭轉速下的喂料量Qin，并與前面所得的仿真輸送量值進行比較，得到不同摩擦因數下的喂料量與輸送量比較圖，如圖10所示. 由圖10可以看出，隨著摩擦因數的增大，輸送量和喂料量之間的差距越來越小，即顆粒的濺出量隨著顆粒與喂料頭之間的摩擦因數的增大而減小.

圖10 不同摩擦因數下的喂料量與輸送量對比圖Fig.10 Comparison chart between feeding quantity and transport quantity under different friction coefficient

從式(22)中可得出，在同一摩擦因數和同一曲面下，喂料量和喂料頭轉速成線性關系，設Qin=aωf，則根據式(22)，可得a=120πRfDfn. 從圖10中可以看出，輸送量和喂料頭轉速的關系也基本成線性關系，設Qt=bωf，對圖13中的每個仿真輸送量值進行擬合可得出b. 設k=b/a，則不同摩擦因數下的k便可以計算出，計算結果如圖11所示.

擬合函數為

k=-0.075μ2+0.273μ+0.782.

(23)

所以喂料頭的輸送量為

Qt=(-0.075μ2+0.273μ+0.782)×

(120πnDfRfωf).

(24)

已知垂直螺旋輸送機的輸送量時，由式(24)，可得喂料頭的轉速為

ωf=

(25)

圖11 不同摩擦因數下的k值Fig.11 Value of k under different friction coefficient

4 結 論

針對相對旋轉式喂料頭進行工作機理研究和動力學分析，分析了其氣固兩相DEM模型及其仿真方法，以某種高效的喂料頭為研究對象，研究喂料頭導擋葉5種曲面的曲線方程，并利用EDEM+Fluent耦合仿真方法分析喂料頭內的顆粒狀態(tài)，比較不同摩擦因數下不同曲面的輸送量的對比情況，得到不同摩擦因數下推薦使用的喂料頭導擋葉曲面. 進行喂料頭內輸送量和喂料量的分析，計算導擋葉5種曲面的喂料窗口的過料面積，進而推導出喂料頭的輸送量公式，得到喂料頭的轉速公式，對喂料頭的設計有非常重要的指導意義.