劉怡琳，李鈺，余亞雄，黃哲慶，周強，2，3

（1 西安交通大學化學工程與技術學院，陜西 西安 710049； 2 新能源系統(tǒng)工程與裝備陜西省高校工程研究中心，陜西 西安710049； 3 西安交通大學動力工程多相流國家重點實驗室，陜西 西安 710049）

引 言

氣固流化床反應器由于其良好的混合、傳質和傳熱性能，在能源化工領域中得到了廣泛應用，掌握氣固兩相流系統(tǒng)的流動和傳熱特性對于提高工業(yè)設備的性能具有重要意義。決定反應器動量傳遞及傳熱傳質的關鍵參數(shù)為曳力、傳熱及傳質系數(shù)。因此國內外有很多研究者針對上述系數(shù)開展了廣泛的理論、實驗及數(shù)值模擬研究，以獲得上述系數(shù)的計算關聯(lián)式。然而，現(xiàn)有研究中一些常用的關聯(lián)式大都是針對于顆粒位置隨機均勻分布的氣固兩相流系統(tǒng)[1]。氣固兩相在非線性相間曳力和固相應力的作用下，很容易產生非均勻的介尺度流動結構，如顆粒聚團等[2-3]。研究表明，計算中未考慮介尺度結構的影響會導致高估氣固相間曳力、傳熱和傳質[4-6]。

由于氣固系統(tǒng)固有的復雜性，系統(tǒng)研究介尺度結構中氣固傳熱現(xiàn)象的實驗較少見諸于文獻[7]。在過去的20 年中，計算流體動力學（computational fluid dynamics,CFD）已經成為模擬氣固系統(tǒng)中流動和傳熱的越來越強大的工具[8]。在模擬中為了解析介尺度結構對流動及傳熱傳質的影響，方法之一是采用足夠小的計算網格（通常為數(shù)倍的顆粒直徑）以滿足近似均勻假設[9-10]，然而，對于大型工業(yè)反應器而言，細網格模擬的計算成本是無法承受的。因此在實際應用中，研究人員通常采用粗網格雙流體法（two-fluid model, TFM）或多相質點網格方法(multi-phase particle-in-cell,MP-PIC)模擬大型流化床，但粗網格模擬無法解析介尺度結構的影響，即微尺度的近似均勻模型在有介尺度結構存在的粗網格模擬中不再適用，這就需要為適合工業(yè)尺度的粗網格模擬方法提供可以考慮亞網格非均勻性的介尺度模型，以得到合理的預測結果[11]。

目前，在不考慮傳熱的氣固兩相流研究中，研究者已揭示了構建介尺度曳力、介尺度應力等模型的必要性[12]，并基于能量最小多尺度方法（energyminimization multi-scale，EMMS）[13-15]、過濾方法[16-18]等開發(fā)了氣固曳力介尺度模型。然而，對介尺度氣固相間傳熱的研究還比較少。Dong 等[19]采用EMMS方法研究了介尺度結構對循環(huán)流化床提升管內氣固傳質的影響，他們發(fā)現(xiàn)，介尺度結構會極大地影響傳質。Hou 等[20]使用EMMS 方法研究了快速流化床中非均勻流動結構與傳遞系數(shù)之間的關系。他們發(fā)現(xiàn)，介尺度結構對曳力系數(shù)、傳質系數(shù)和傳熱系數(shù)有很大影響，他們指出，傳統(tǒng)的均勻模型在有非均勻流動結構存在時不再適用。魯波娜等[21]基于EMMS多尺度動量傳遞及傳熱模型對多產異構烷烴催化裂化反應器進行了模擬，研究表明，相比于傳統(tǒng)的模型，采用多尺度模型可以更為準確預測反應器內的流動結構及溫度分布。Shu 等[22]考慮了顆粒聚團結構對傳熱的影響，對傳統(tǒng)傳熱系數(shù)模型進行了改進，采用改進的氣固曳力模型及氣固傳熱系數(shù)模型對下行床反應器進行了模擬，研究表明采用優(yōu)化模型可以定性捕捉下行床的關鍵傳熱特性。Guo等[23]基于計算流體力學-離散單元法（computational fluid dynamics-discrete element method, CFD-DEM）研究了聚團在氣固傳熱中的作用，他們定義了一個傳熱修正量來量化由聚團引起的局部不均勻性的影響，并提出了一個基于假定形狀概率密度分布函數(shù)的統(tǒng)計模型。Lane 等[24-25]通過過濾雙流體模型（TFM）的計算結果，開發(fā)了一個用于浸沒水平圓柱體的氣固兩相流傳熱的亞網格模型。在他們的模型中，采用固含率、固相速度、圓柱體幾何形狀（直徑和間距）和Peclet 數(shù)來封閉Nusselt 數(shù)。Rauchenzauner 等[26]提出了一個基于漂移溫度的過濾相間傳熱模型。

Agrawal 等[27]構建了基于細網格TFM 模擬的過濾相間傳熱模型。他們指出，過濾后的介尺度相間傳熱系數(shù)比微觀相間傳熱系數(shù)小1～2 個數(shù)量級，他們將微觀相間傳熱系數(shù)的修正量（記為Q）建模為過濾固含率和過濾器尺寸的函數(shù)，在其模型中，過濾尺度小于等于16.7dp（dp為顆粒直徑）時不再對傳熱系數(shù)進行修正（Q=0）。Huang 等[28]通過過濾細網格TFM模擬數(shù)據(jù)，將Q構建為過濾固含率、過濾器尺寸以及兩相之間的無量綱過濾溫差的函數(shù)，在其模型中，過濾尺度小于等于8.3dp時不再對傳熱系數(shù)進行修正（Q= 0）。Li 等[29]對Huang 等的模型進行了修正，取消了Q在（0,1）之間的限制。Lei等[30]通過過濾CFD-DEM 數(shù)據(jù)，構建了介尺度傳熱系數(shù)修正量與過濾固含率、過濾器尺寸及氣固相間溫差的函數(shù)關系。但Huang等[28]、Li等[29]模型中的無量綱過濾溫差需要迭代獲得，在實際的粗網格應用中存在困難，Lei 等[30]的模型對過濾溫差的無量綱方法（除以0.0001 K 使得大部分的氣固溫差處于0～10 之間）缺乏普適性，同樣導致模型應用困難。

另外，在現(xiàn)有的傳熱介尺度模型研究中，如何保持整個系統(tǒng)的熱平衡是一個非常復雜的問題。為了保證氣固兩相間存在持續(xù)的相間溫差，現(xiàn)有文獻中使用了不同的方法。Rauchenzauer 等[26]在每個時間步讓氣相溫度保持線性增加狀態(tài)并在固相中設置熱匯。Lane等[24-25]和Lei等[30]在固相中添加了熱源項，熱源項分別以1 K/s和0.1 K/s的速率加熱氣固系統(tǒng)。他們指出，加熱速率的值是任意的，其大小不會影響模型結果。Agrawal 等[27]和Huang 等[28]在固相和氣相中分別添加了熱源和熱匯項，他們也指出，所構建模型對所施加的熱源熱匯項的大小不敏感。以上方法可以實現(xiàn)氣固相間溫差的維持，然而如何維持氣固相間溫差以獲得更符合真實物理過程（氣固兩相自由傳熱）的介尺度傳熱模型仍是一個問題。

針對以上問題，本文提出了一種重置溫度的方法以維持氣固溫差，使得氣固兩相可自由換熱，并基于細網格CFD-DEM 數(shù)據(jù)過濾獲得了粗網格方便使用的介尺度氣固相間傳熱系數(shù)修正模型。

1 細網格CFD-DEM 傳熱計算

CFD-DEM 方法是基于歐拉-拉格朗日坐標系的一種離散模擬方法，在歐拉坐標系下求解氣相平均方程，在拉格朗日坐標系下求解顆粒方程，分別追蹤單個顆粒的受力、運動及傳熱。有研究表明，CFD-DEM 采用小至1.75～3倍顆粒直徑（1.75dp～3dp）的網格系統(tǒng)可以獲得與解析到顆粒表面的直接數(shù)值模擬（particle resolved-direct numerical simulation,PR-DNS）近似的計算結果[31-32]。因此細網格（1.75dp～3dp）CFD-DEM 計算可以用來為粗網格計算提供非均勻介尺度封閉模型。本文的CFD-DEM 模擬通過對開源軟件MFIX（版本19.2.2，https://mfix.netl.doe.gov/)進行二次開發(fā)實現(xiàn)，顆粒相和氣相的動量、質量控制方程及求解過程可參考Garg 等[33]的研究。顆粒為典型的Geldart A 類顆粒，擬二維周期邊界計算域垂直方向長度Ly為960dp，橫向Lx為240dp，展向Lz為6dp，計算網格數(shù)量為384×96×2，對應x,y方向的網格尺寸為2.5dp。為了平衡周期域上的重力，在垂直方向上設置了恒定氣體流量邊界條件。整體固含率為0.05，初始時刻顆粒在計算域中隨機均勻分布。為了保證氣固兩相間存在持續(xù)的相間溫差，本文中選取了兩種方法：方法一，參考文獻中的方法給氣相能量方程添加熱源項，該熱源項以0.1 K/s 的速度加熱氣體，該設定值的大小不影響計算結果[27-28]；方法二，每隔固定時間重置氣體溫度，所間隔時間為1為顆粒弛豫時間，單位為s，=18μg，其中，ρp為顆粒密度；dp為顆粒直徑；μg為氣相黏度），氣固兩相在1τstp內自由換熱。在計算初始時，氣相溫度設置為1000 K，固相溫度設置為0 K，氣固兩相間傳熱由能量方程控制，氣相的能量方程如式(1)所示，具體可參考Lei 等[30]和Hou 等[34]的研究。

計算持續(xù)50τ，每0.2輸出一次結果，計算后全場平均滑移速度、平均氣固溫差等統(tǒng)計數(shù)據(jù)達到穩(wěn)定狀態(tài)，即系統(tǒng)達到統(tǒng)計穩(wěn)態(tài)[36-37]，后30τstp的數(shù)據(jù)用于模型構建。模型通過對計算結果進行過濾、擬合獲得，所采用的空間過濾方法與Agrawal等[27]相同，具體可參考文獻[9,27]，研究中所采用的過濾尺度為5dp、10dp、17.5dp、20dp、25dp和40dp。

表1 主要計算參數(shù)及設置Table 1 Parameters and settings

2 結果與討論

與文獻[27-29]中相同，本文針對介尺度傳熱系數(shù)修正量Q進行建模。Q定義為：

2.1 兩種維持氣固溫差方法的討論

為了便于表述，將方法一稱為熱源項方法，方法二稱為重置溫度方法。圖1為兩種溫度設置方法下的固含率分布，即特定固含率對應的網格數(shù)量占總計算網格數(shù)量的百分比?？梢钥吹?，在本計算中，固含率在（0.025,0.03）范圍內的網格數(shù)量較多，存在峰值，兩種方法下固含率占比分布基本一致，流動結構基本一致，在此前提下開展兩種維持氣固溫差的方法的對比和討論是有效的。

圖1 兩種方法下的固含率分布Fig.1 The distribution of solid volume fraction under two methods

圖2為兩種方法下計算域中平均氣固溫差隨時間的變化，由圖可知，兩種方法均可以使得系統(tǒng)中存在一定的氣固溫差。方法一可以使得系統(tǒng)中氣固溫差在計算達到穩(wěn)定后處于較為穩(wěn)定的狀態(tài)，波動幅度不大；方法二則是在每次重置氣體溫度后，氣固溫差逐漸下降，直至下一次重置氣體溫度。值得注意的是，兩種方法得到的平均氣固溫差的絕對值不同，但在氣固兩相傳熱穩(wěn)定后該溫差的大小不影響Q的統(tǒng)計結果[27-28]。

圖2 兩種方法下計算域中平均氣固溫差隨時間的變化Fig.2 The average gas-solid temperature difference in the calculation domain versus time under the two methods

為了對比兩種方法對計算域內氣固兩相傳熱的影響，本文統(tǒng)計了不同局部固含率的傳熱量占計算域總傳熱量的百分比，以及局部單位體積傳熱量（某網格內的傳熱量除以網格體積）與總單位體積傳熱量（總傳熱量除以總體積）之比，如圖3、圖4 所示[為了凸顯兩種方法在大固含率下的區(qū)別，對縱坐標取對數(shù)，如圖3(b)和圖4(b)所示]。由圖3 可知，兩種方法下傳熱量占比均在0.025≤?s≤0.03 時存在峰值，即兩種方法下占總網格數(shù)11.0%（圖1）的固含率處于0.025≤?s≤0.03 的網格分別貢獻了25.5%（重置溫度方法）和18.65%（熱源項方法）的傳熱量。由圖4 可知，局部單位體積傳熱量與總單位體積傳熱量之比同樣在0.025≤?s≤0.03時存在峰值，重置溫度方法中該固含率范圍內局部單位體積傳熱量是總單位體積傳熱量的2.32 倍，熱源項方法中該固含率范圍內局部單位體積傳熱量是總單位體積傳熱量的1.71 倍。重置溫度方法在固含率較小時（?s< 0.03）局部單位體積傳熱量與總單位體積傳熱量之比大于熱源項方法，而在固含率?s> 0.03 的范圍內該比值小于熱源項方法，二者差異在固含率大于0.1 時尤為明顯。

圖3 計算域內總傳熱量在不同局部固含率下的分布Fig.3 The distribution of total heat transfer in the calculation domain under different local solid volume fraction

圖4 局部單位體積傳熱量與總單位體積傳熱量之比Fig.4 The ratio of local gas-solid heat transfer per unit volume to the total gas-solid heat transfer per unit volume

為了分析0.025≤?s≤0.03 的特殊性，圖5 展示了某時刻的顆粒位置圖和固含率云圖，可以看到0.025≤?s≤0.03（圖5右側圖片綠色部分）的網格大多處于顆粒聚團邊界位置，由此，下文將固含率范圍為0.025≤?s≤0.03、?s<0.025和?s>0.03的區(qū)域分別稱為界面、稀相和濃相區(qū)域。事實上，界面所對應固含率范圍受整體固含率影響較大，本文整體固含率為0.05，所以界面所對應的固含率值較小。結合圖3～圖5 可以發(fā)現(xiàn)，兩種方法均表明聚團界面位置的局部單位體積氣固傳熱量與總單位體積傳熱量之比最大，重置溫度方法在稀相和界面位置的局部單位體積傳熱量與總單位體積傳熱量之比大于熱源項方法，而在濃相位置的局部單位體積傳熱量與總單位體積傳熱量之比小于熱源項方法。為了更直觀地了解本文所研究的系統(tǒng)中傳熱的分布情況，圖6 統(tǒng)計了稀相、界面和濃相的網格數(shù)量占比（除去固含率為0的網格）及傳熱量占比，可以看到兩種方法的固含率占比一致，而傳熱量占比卻有所區(qū)別。兩種方法的結果中，濃相、界面和稀相的網格占比分別為55.3%、11%和33.7%，而其對應的傳熱量占比卻差別顯著。在重置溫度方法下，濃相、界面和稀相分別貢獻了36.7%、25.5%和37.8%的傳熱量；熱源項方法下，濃相、界面和稀相分別貢獻了55.4%、18.6%和26.0%的傳熱量?？偟貋碚f，重置溫度方法下稀相和界面的傳熱量占比高于熱源項方法，而濃相的傳熱量占比則低于熱源項方法。究其原因，熱源項方法在氣相中添加的熱源（Sg）由式（3）計算：

圖5 某時刻顆粒位置及固含率分布Fig.5 Particle position and distribution of solid volume fraction at a certain time

圖6 稀相、界面和濃相的網格數(shù)量占比及傳熱量占比Fig.6 The percentage of grid number and heat transfer of dilute phase,interface and dense phase

式中，Π?為氣相的能量輸入速率，該設定值的大小不影響計算結果；?sd為計算域整體固含率；?gd為計算域整體氣含率。由式（3）可知，氣相所加熱源項的大小與所在網格的氣含率（?g）有關，不論網格內真實傳熱過程如何，均有熱源加入促使其傳熱。當體系中存在化學反應時，熱源項方法可能與實際較為吻合；當體系中不存在化學反應時，在濃相區(qū)域中氣固相間傳熱一般會較快完成致使氣固相間溫差較小進而傳熱量減少，但熱源項方法會強迫濃相區(qū)域保持一定的溫差進而促進濃相區(qū)域的氣固相間傳熱。重置溫度方法保持了氣固相間自由傳熱特性，僅通過反復重置溫度獲得足夠的統(tǒng)計數(shù)據(jù)來完成建模。在無化學反應的氣固傳熱系統(tǒng)中，氣固兩相間由于溫差的存在會在能量方程控制下自發(fā)傳熱，該過程與本文所提出的重置溫度方法所描述的物理過程類似，因此采用重置溫度方法所構建的介尺度傳熱模型可能會具有較高的準確度。而在有化學反應的氣固傳熱系統(tǒng)中，由于反應和流動、傳熱之間存在復雜的耦合關系，熱源項方法所構建模型是否能較好地近似系統(tǒng)內的傳熱過程尚需更多的模擬或實驗來進行分析和驗證。

2.2 介尺度氣固傳熱模型構建

圖8 為兩種方法在不同尺度下過濾得到的Q的平均值隨過濾固含率的變化，兩種方法下獲得的Q有所不同，在固含率較大時，二者差異更為明顯，整體上重置溫度方法的Q大于熱源項方法獲得的Q。

基于圖8 的數(shù)據(jù)，可將Q構建為過濾固含率和過濾尺度的函數(shù)，重置溫度方法獲得的模型擬合效果如圖9所示，模型為：

圖7 重置溫度后不同時刻的過濾Q隨過濾固含率的變化Fig.7 The change of filtered Q with the filtered solid volume fraction at different time after resetting the temperature

圖8 兩種方法在不同尺度下獲得的Q的平均值隨過濾固含率的變化Fig.8 The average value of Q obtained by the two methods varied with the filtered solid volume fraction at different filter sizes

圖9 重置溫度方法下介尺度模型擬合Q和真實Q隨過濾固含率的變化Fig.9 The variation of fitting Q obtained by the mesoscale model and real Q with the filtered solid volume fraction under the resetting temperature method

值得注意的是，由圖9 可知，當過濾尺度較大時，如為25dp及40dp時，過濾固含率的最大值小于小過濾尺度的最大值。造成此現(xiàn)象的原因為計算區(qū)域有限（240dp× 960dp），過濾尺度為40dp時已占據(jù)了橫向尺寸的1/6，很難出現(xiàn)粗網格上的高固含率。加之全場整體固含率為0.05，隨著過濾尺度的增大，粗網格內的固含率會逐漸接近整體固含率。

2.3 介尺度氣固傳熱模型先驗分析

本節(jié)將對上述模型進行先驗分析，評價指標采用Pearson相關系數(shù)[36]，其計算式為：

式中，x、y分別為模型預測值（Qpredict）的數(shù)據(jù)集和粗網格過濾值（Qexact）的數(shù)據(jù)集；xˉ、yˉ分別為數(shù)據(jù)集的平均值。采用本文所構建的模型[式(4)]對不同過濾尺度和固含率下的Q進行預測，預測值（Qpredict）與基于細網格CFD-DEM 模擬的粗網格過濾值（Qexact）的Pearson 相關系數(shù)如圖10 所示，本文所構建模型的表現(xiàn)明顯優(yōu)于Agrawal 等[27]的模型，而且隨著過濾尺度的增大模型表現(xiàn)在持續(xù)變好。這是因為本文的模型是基于細網格CFD-DEM（網格尺寸為2.5dp）的模擬數(shù)據(jù)，而Agrawal 等[27]的細網格模擬采用的網格尺寸為16.67dp，所構建的模型是基于至少2 倍的基礎網格即33.34dp的數(shù)據(jù)，因此理論上模型適用范圍為過濾尺度大于等于33.34dp，因此導致其模型在小過濾尺度時表現(xiàn)不如本文所構建的模型。更重要的是，采用網格大小為16.67dp的模擬數(shù)據(jù)作為建模數(shù)據(jù)，本身已完全忽略了比16.67dp尺度更小的非均勻結構對更大尺度上流動和傳熱的影響，因此所構建模型準確度較低。事實上，大量研究表明，對氣固兩相流動進行完全解析需要在3dp左右的網格尺寸[38-39]。另外，與文獻中所采用的熱源項方法不同，本文所構建的模型基于重置溫度方法，數(shù)據(jù)源的不同也是導致模型表現(xiàn)不同的原因之一。最后，相較于TFM，CFDDEM 可以追蹤顆粒的運動，計算得到的介尺度結構更為準確。圖11 為所構建模型在不同過濾尺度下模型預測值（Qpredict）與粗網格過濾值（Qexact）相對誤差的概率密度分布?？梢钥吹诫S著過濾尺度的增大，相對誤差逐漸向0 集中，即在模型應用中，在本文所研究的網格尺寸范圍內，所采用的粗網格尺寸越大，模型的表現(xiàn)越好。本文所構建的模型可適用于過濾尺度小于等于40dp時的粗網格氣固傳熱計算。另外需要指出的是，本文所構建的模型是基于周期性邊界的模擬數(shù)據(jù)，未考慮壁面效應，在實際應用過程中，若壁面剪切效應較強，則可能會存在誤差。

圖10 不同過濾尺度下模型預測值（Qpredict）與粗網格過濾值（Qexact）的Pearson相關系數(shù)Fig.10 Pearson correlation coefficients between Qpredict and Qexact under different filter sizes

圖11 不同過濾尺度下模型預測值（Qpredict）與粗網格過濾值（Qexact）相對誤差的概率密度分布Fig.11 The probability density distribution（PDF）of relative errors between Qpredict and Qexact at different filter sizes

3 結 論

本文采用CFD-DEM 研究了氣固兩相流相間傳熱問題，比較了兩種維持氣固相間傳熱溫差的方法的優(yōu)缺點。方法一中給氣相能量方程添加了熱源項；方法二中每間隔一段時間進行氣相溫度重置，重置后氣固兩相進行自由傳熱?；谥刂脺囟确椒ǖ倪^濾數(shù)據(jù)構建了適用于粗網格計算的氣固相間介尺度傳熱模型，先驗分析表明所構建模型具有優(yōu)越性。具體結論如下。

（1）聚團界面位置的局部單位體積氣固傳熱量占比最大，重置溫度方法在稀相和界面位置的局部單位體積傳熱量與總單位體積傳熱量之比大于熱源項方法，而在濃相位置該比值小于熱源項方法。分析表明重置溫度方法更適合計算氣固相間自由傳熱的情況。

（2）基于過濾固含率和過濾尺度構建了介尺度氣固相間傳熱系數(shù)修正因子的封閉模型，所構建模型的Pearson相關系數(shù)隨過濾尺度的增大而增大，所構建模型預測得到的Q與真實Q的Pearson 相關系數(shù)高于文獻中已有的雙參數(shù)模型，模型適用于粗網格尺度小于等于40 倍顆粒直徑的氣固兩相流傳熱計算。

由于介尺度傳熱模型的后驗十分復雜，涉及介尺度曳力模型的選擇與驗證，目前文獻中尚未有介尺度傳熱模型的后驗報道。后續(xù)將繼續(xù)探索介尺度傳熱模型的后驗方法，開展相關的后驗研究，以進一步評估所構建介尺度傳熱模型在實際應用中的表現(xiàn)。