陳桂香，劉文磊，劉超賽，鄭德乾，岳龍飛

(河南工業大學土木工程學院，鄭州 450001)

固定床是一種充滿隨機填充顆粒的固定容器，廣泛應用于農業干燥、化學反應以及核反應堆之中[1-3]。糧食籽粒作為常見的顆粒物質，在糧堆固定床中存在隨機堆積過程，使床層內的孔隙分布十分復雜，對床層內氣體流動、溫度傳遞過程產生影響，因此研究固定床內部流體的流動阻力與傳熱特性，對提高糧食干燥效率與固定床的設計具有指導作用[4, 5]。

國內外學者對糧堆氣流阻力與傳熱特性進行大量的實驗與模擬研究。陳竹筠[6]利用設計的實驗臺，研究了不同流速、孔隙率以及床層高度對玉米堆通風阻力的影響，并對Ergun方程進行了修正。朱英開等[7]采用均勻的多孔介質假定對糧層阻力進行模擬，得到了3種糧食的壓降分布規律。Kashaninejad等[8]對散裝開心果進行氣流阻力測定，研究了填充方式等對結果的影響，對比3種氣流阻力模型，Ergun方程具有最佳的擬合效果。Zare等[9]設計了深床干燥實驗裝置，通過改變入風時氣流溫度與流量，測定了糧食溫度變化規律，給出了稻谷傳熱效率的計算方式。Ranjbaran等[10]采用計算流體力學(CFD)方法對深床干燥器中的稻谷進行數值模擬，利用編寫的自定義函數，對深床稻谷傳熱過程及傳熱性能進行了研究。

多數學者采用宏觀水平的多孔介質假定模型對糧堆阻力進行模擬研究，沒有考慮糧食顆粒在隨機堆積過程中孔隙分布不均勻對氣流特性的影響。黃凱[11]利用EDEM軟件構建玉米顆粒進行隨機堆積，根據孔隙結構參數信息，建立了孔道網絡通風模型。Rusinek等[12]采用離散單元法(DEM)預測了油菜籽在筒倉中的溫度分布，描述了油菜籽的自發熱過程。劉立意等[13]采用CFD-DEM耦合方法，對2個品種稻谷的通風阻力進行數值模擬分析，得到了糧層阻力與糧層深度、表觀風速的關系，但沒有直接將孔隙結構導入到流體計算過程中。離散元(DEM)方法已廣泛地應用于顆粒介質的堆積、壓縮等研究中，基于DEM-CFD耦合的方法，能更加真實的模擬堆積床中氣流流動與傳熱過程。

本研究通過離散單元法(DEM)構建所需孔隙率的球床堆積結構，將球床堆積信息導入前處理軟件中并對顆粒間的接觸進行處理，利用CFD模擬大豆固定床中的氣體流動與傳熱過程，研究不同入口風速對糧堆氣流阻力與傳熱特性的影響，為后續糧堆非均質孔隙結構中的多場耦合研究提供參考。

1 數值模擬模型及方法

1.1 離散元堆積模型

本研究選用PFC3D離散元軟件模擬了顆粒的隨機生成過程，與連續多孔介質模型相比，離散單元法(DEM)將顆粒材料視為不同粒子的集合，球形顆粒由于重力向下運動，每個粒子均遵循基本的物理和力學定律，考慮了粒子物理特性與相互作用特性，當顆粒達到穩定時，模擬堆積過程停止[14]。本研究設定顆粒的運動速度小于10-7m/s時達到穩定，共生成1 418個顆粒，表1為離散元模擬所需的參數值[15]。

表1 大豆顆粒堆積床模擬的離散元參數

在固定床層的實驗與模擬研究中，需要研究壁面效應對結果的影響，Mueller等[16]研究表明，固定床直徑D/dp≤10時，壁面效應對堆積填充的孔隙率影響較大。圖1為本研究的幾何模型，由三部分組成：空氣入口段，顆粒堆積段、空氣出口段。出口與入口段的長度為20 mm，堆積段長度為100 mm，固定床層的直徑為70 mm，顆粒直徑為6.5 mm(D/dp=10.7)。

圖1 幾何模型示意圖

1.2 網格劃分及獨立性驗證

1.2.1 網格劃分

采用自開發的腳本程序，生成上述構建的隨機堆積顆粒與幾何模型，重建出的球形顆粒之間會形成點接觸或狹小的間隙，不利于網格的生成，直接進行網格劃分時，接觸點附近的網格單元會有較大的偏斜率，導致模擬難以收斂，影響計算精度，因此需要對重建顆粒間的接觸進行處理。現有的主要處理方法包括縮徑法、擴徑法、搭橋法和截面法[17]。縮徑與擴徑法改變了球體的半徑，使整體的孔隙率發生變化，截面法與搭橋法只對接觸點局部進行處理，保證了球床幾何模型的完整性。本研究選用搭橋法進行接觸處理，能更準確的模擬動與傳熱的過程，圓柱形橋的直徑為0.1dp，高度為0.05dp，網格劃分與接觸處理如圖2所示。

圖2 網格劃分與接觸處理圖

1.2.2 網格獨立性驗證

為保證模擬結果的準確性，進行網格獨立性研究，并選擇合適的網格尺寸。采用非結構四面體網格對模型進行劃分，設置4種網格尺寸分別為1/5、1/10、1/15dp和1/20dp。在同一模擬設置條件下，網格獨立性計算結果如圖3所示，當網格尺寸為1/15dp時，固定床中定點處流體壓降和溫度計算結果無明顯差別，滿足無關性要求。因此，本次模擬選用最大網格尺寸1/15dp進行研究，網格數量為612萬。

圖3 網格獨立性驗證

1.3 CFD計算模型

本研究采用ANSYS Fluent對流動與傳熱進行計算求解，模型計算區域的進口和出口分別設置為速度入口邊界和壓力出口邊界，側壁與糧食顆粒邊壁設置為固定邊界。Dixon等[18]對固定床層中的流動進行了細致的劃分，研究表明當顆粒雷諾(Reynolds)數Rem大于300時，床層內的流動屬于湍流，在本次模擬中入口風速設置在0.5～2 m/s范圍內，Rem為358～1 432，由于固定床內部通道復雜，故選用RNG-k-ε湍流模型，壁面進行標準函數處理。控制方程采用基于壓力-速度方法進行穩態求解，壓力-速度耦合求解器采用SIMPLE算法，壓力項采用PRESTO!格式離散，其余方程采用二階迎風格式離散，為了得到更精確的結果，所有變量的收斂殘差設置為10-5，表2為模擬初始參數條件。

表2 CFD數值模擬參數

2 計算結果分析

2.1 徑向孔隙率分布

為了得到徑向孔隙率分布，將床層沿徑向劃分為多個同軸圓柱形，如圖4所示。

圖4 球床徑向剖面圖

孔隙率為圓柱面截取的流體面積與圓環面積的比值，計算方法如式(1)所示。

(1)

式中：As為徑向半徑r處圓柱面上固體區域的面積/m2；H為顆粒堆積區域的高度/m。

通過Fluent中的切面工具沿固定床徑向設置45個同心圓柱面，利用Report工具計算出柱面上的孔隙率，離散元(DEM)堆積結構徑向孔隙分布如圖5所示。De Klery[19]提出了一種廣泛用于預測圓柱形填充床內孔隙率分布的經驗公式，如式(2)所示。

(2)

由圖5可知，根據離散元(DEM)堆積結構計算出的徑向孔隙率分布與De Klery[19]提出的經驗公式基本一致，徑向孔隙率向床層中心方向振蕩，離中心距離越近振幅越小，并趨于恒定值。當無量綱距離n=0(r=R)時，徑向孔隙率最大為1，這是因為最外側圓柱面與球顆粒只存在接觸點，隨著無量綱距離增大，孔隙率逐漸變小，最終趨近于球床平均孔隙率。

圖5 球床徑向孔隙率分布

2.2 壓力場分析

糧堆壓力降是固定床中流動特性的重要特征參數，它受顆粒雷諾數與壁面阻力系數的影響。修正后的顆粒雷諾數計算公式為：

(3)

式中：Rem為顆粒雷諾數；ρ為流體密度/kg/m3；u為入口流速/m/s；dp為顆粒直徑/m；ε為球床孔隙率；μ為流體的動力粘度/(Pa·s)。

Ergun方程是用于計算多孔介質中阻力系數與壓力降最為廣泛的方法[20]，Reichelt[21]考慮壁面效應影響，對Ergun方程的相關系數進行修正，圖6為不同入口風速下顆粒雷諾數與阻力系數值。

圖6 不同入口風速下顆粒雷諾數與阻力系數

本研究選用Ergun方程與Reichelt方程驗證DEM-CFD數值模擬的正確性。具體如式(4)～式(8)所示。

(4)

(5)

(6)

(7)

Bw=[1.15(dp/D)2+0.87]2

(8)

圖7為不同入口風速下球床壓力降模擬結果與Ergun、Reichelt方程計算結果對比。DEM-CFD模擬結果與公式計算結果基本一致，糧層壓力降與入口速度近似呈二次函數關系，擬合公式見式(9)。

(9)

在入口速度小于1 m/s時，模擬結果與兩種方程計算結果較為吻合；當入口速度大于1 m/s時，模擬結果與Reichelt方程的計算結果平均誤差為5.41%，由于Reichelt方程考慮了壁面效應對壓力降的影響，與Ergun方程計算結果相差較大。

圖7 不同入口風速下球床壓力降

由圖8可知，4種入口風速下壓力分布趨于一致，沿出口方向壓力逐漸降低，壓力分布沿床層高度呈現明顯分層現象。對比4組模擬結果，隨著入口速度的增加，出口處的壓力下降的越快，與壓力降計算結果一致。在進口與出口段處壓力分布恒定，床層內部由于顆粒的阻擋，使得球床內局部壓力變化明顯。

2.3 速度場分布

由圖9可知，4種入口風速下形成的速度場分布相似，隨著入口速度的增加，球床內部風速也越大。沿球床軸向方向上，入口處風速與填充段前段顆粒接觸時，顆粒會對氣流形成阻擋，使入口風速不均勻分布；在填充段內部，由于顆粒的隨機分布流體通道十分復雜性，使得填充段內速度分布存在明顯差異，局部平均風速增大到入口風速的2.5倍；在出口端，氣流會出現停滯與回流現象。沿球床徑向方向上，近壁面處孔隙率最大，向床層中心方向逐漸減小，可以發現孔隙度較高區域的氣流速度明顯高于其它區域。

2.4 溫度場分布

由圖10可知，隨著入口風速的增加，球床內溫度下降區域越大。由于近壁面處空氣流速高于床層中心區域，流體與近壁面區域的顆粒表面熱交換更為充分，近壁面處流體的平均溫度低于中心區域，在低速條件下，壁面處與中心區域的溫度差更為明顯。球床內部孔道彎曲復雜，流體溫度下降幅度較小，在出口段處，由于流體速度的停滯與回流，出口段流體溫度明顯高于堆積內部溫度。

3 結論

本研究通過離散元法(DEM)在固定床層內生成大豆顆粒的隨機堆積結構，并使用CFD方法直接模擬所創建的顆粒與流體區域。研究了不同入口速度下床層內的流動與傳熱特性，包括床層孔隙率分布、速度場分布、壓力降以及溫度分布等規律。結果表明：

圖8 不同入口風速下球床壓力分布云圖

圖9 不同入口風速下X=0截面處速度分布云圖

圖10 不同入口風速下X=0截面處溫度分布云圖

固定床近壁面處孔隙率最大并趨近于1，徑向孔隙率向床層中心方向呈振蕩趨勢，離床層中心距離越近，孔隙率與振幅逐漸減小，最終趨近于球床平均孔隙率。隨著入口風速增加，壓力降增加，糧堆壓力降與入口速度近似呈二次函數關系，考慮了壁面效應的Reichelt方程較Ergun方程對壓力降的計算更準確。固定床中氣流速度與溫度場分布受孔隙率分布、入口風速影響較大，床層內部速度與溫度場分布不均勻。近壁面區域氣流速度是床層中心速度的2.5倍，相比于中心區域，近壁面處流體的平均溫度降低更快，在出口處會出現氣流的停滯區。