基于DEM-CFD的糧堆氣流阻力與傳熱特性模擬研究

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

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

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

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

多數(shù)學(xué)者采用宏觀水平的多孔介質(zhì)假定模型對糧堆阻力進行模擬研究，沒有考慮糧食顆粒在隨機堆積過程中孔隙分布不均勻?qū)饬魈匦缘挠绊憽｜S凱[11]利用EDEM軟件構(gòu)建玉米顆粒進行隨機堆積，根據(jù)孔隙結(jié)構(gòu)參數(shù)信息，建立了孔道網(wǎng)絡(luò)通風(fēng)模型。Rusinek等[12]采用離散單元法(DEM)預(yù)測了油菜籽在筒倉中的溫度分布，描述了油菜籽的自發(fā)熱過程。劉立意等[13]采用CFD-DEM耦合方法，對2個品種稻谷的通風(fēng)阻力進行數(shù)值模擬分析，得到了糧層阻力與糧層深度、表觀風(fēng)速的關(guān)系，但沒有直接將孔隙結(jié)構(gòu)導(dǎo)入到流體計算過程中。離散元(DEM)方法已廣泛地應(yīng)用于顆粒介質(zhì)的堆積、壓縮等研究中，基于DEM-CFD耦合的方法，能更加真實的模擬堆積床中氣流流動與傳熱過程。

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

1 數(shù)值模擬模型及方法

1.1 離散元堆積模型

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

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

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

圖1 幾何模型示意圖

1.2 網(wǎng)格劃分及獨立性驗證

1.2.1 網(wǎng)格劃分

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

圖2 網(wǎng)格劃分與接觸處理圖

1.2.2 網(wǎng)格獨立性驗證

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

圖3 網(wǎng)格獨立性驗證

1.3 CFD計算模型

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

表2 CFD數(shù)值模擬參數(shù)

2 計算結(jié)果分析

2.1 徑向孔隙率分布

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

圖4 球床徑向剖面圖

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

(1)

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

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

(2)

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

圖5 球床徑向孔隙率分布

2.2 壓力場分析

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

(3)

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

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

圖6 不同入口風(fēng)速下顆粒雷諾數(shù)與阻力系數(shù)

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

(4)

(5)

(6)

(7)

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

(8)

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

(9)

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

圖7 不同入口風(fēng)速下球床壓力降

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

2.3 速度場分布

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

2.4 溫度場分布

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

3 結(jié)論

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

圖8 不同入口風(fēng)速下球床壓力分布云圖

圖9 不同入口風(fēng)速下X=0截面處速度分布云圖

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

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