熊 勃，張 浩，郭鵬越，安希忠

（東北大學1.冶金學院；2.教育部多金屬礦山生態(tài)冶金重點實驗室，沈陽 110819）

隨著社會的發(fā)展，對能源的需求越來越大，尋找新能源來替代化石已經(jīng)刻不容緩.生物質(zhì)作為一種新興能源，由于其儲量巨大且在利用過程中對環(huán)境友好，受到人們的廣泛關注［1］.熱化學氣化是目前公認的最具商業(yè)前景的生物質(zhì)利用技術［2］.然而利用此技術處理濕生物質(zhì)時，需要對生物質(zhì)進行加熱干燥，這個過程大量耗能.超臨界水氣化制氫技術克服了這個問題，可以直接處理高含濕率的生物質(zhì)［3］，而且采用超臨界水流化床反應器制氫具有高氣化效率、易搭載催化劑顆粒等諸多優(yōu)勢［4］.目前，該項技術已應用于工業(yè)實踐，它對改善我國能源短缺、過于依賴化石能源的現(xiàn)狀，對治污降霾、改善環(huán)境等具有重大意義.

超臨界水給工業(yè)應用帶來巨大收益的同時，眾多新穎的研究課題也擺在廣大研究者面前，急需攻關解決.反應器中超臨界水與固體顆粒間的受力和傳熱特性就是其中之一.越來越多的研究結(jié)果表明［5－7］，經(jīng)典流化床受力與傳熱理論在超臨界水流化床中不再適用.這些理論知識的缺失使得對反應器的設計和數(shù)學建模無據(jù)可依.本文的研究目的就在于完善超臨界水流化床傳熱和受力理論.

由于超臨界水高溫高壓的特點，導致采用實驗的方法來研究超臨界水與固體顆粒間的受力和傳熱特性變得困難重重.數(shù)值模擬能夠在一定程度上克服極端實驗條件的限制，時下對顆粒－流體兩相流動的研究存在兩類比較流行的模型［8］，一種是雙流體模型（TFM），它把固體顆粒和流體都當作連續(xù)相［9］；另一種是CFD-DEM模型，它把固體顆粒當作離散相，而流體作為連續(xù)相［10－12］.但不管是哪種模型，流固兩個求解器間均需通過作用力和熱通量實現(xiàn)數(shù)據(jù)的交換.因而，曳力系數(shù)（Cd）和平均努塞爾數(shù)（Nu）是必須要考慮的核心參數(shù)，它們可由式（1）計算得出，其作用是封閉上述模型［13］.

式中，fd為曳力，N；ρ為流體密度，kg／m3；dp為等體積當量球直徑，m；he為流體對流換熱系數(shù)，W／（m2·℃）；A為迎風面積，m2；uc為入口處流體速度，m／s；Q為熱通量，W／m2；k為流體的導熱系數(shù)，W／（m·K）；Ts和Tf分別為顆粒和流體的溫度，K.

在數(shù)值模擬過程中，通過關聯(lián)式（經(jīng)驗公式）確定上述兩個參數(shù)是常用的方法.為了保證模擬的準確性，Cd和Nu的關聯(lián)式必須要能描述各因素對能量傳遞的影響，如雷諾數(shù)（Re）［14］、顆粒形狀［15］、顆粒運動狀態(tài)［16－17］、傾角［18］等.除此之外，墻壁因子（反應器直徑D與顆粒當量直徑dp之比，λ）［19－23］和墻壁邊界條件［19］也是必須要考量的因素.相關學者對固體顆粒在超臨界水中的動量、能量傳遞特性做了數(shù)值模擬研究，也構建了部分關聯(lián)式［24－27］.Xiong等［28］研究了壁面效應對球形顆粒在超臨界水中受力與傳熱的影響并構建了Cd和Nu的關聯(lián)式，但沒有考慮顆粒形狀的影響.本文將采用顆粒尺度數(shù)值模擬方法，通過對不同形狀橢球顆粒在不同λ和反應器壁邊界條件下超臨界水中的強制對流模擬，探究流化顆粒（反應器壁保持靜止）和沉降顆粒（反應器壁隨流體移動）兩類問題［19］，分別得到了相應工況下的Cd和Nu，再通過回歸分析的方法得到兩個參數(shù)的關聯(lián)式，該關聯(lián)式可為CFD-DEM等宏觀多相流模擬所采用.

1 數(shù)學模型

本文采用的是三維穩(wěn)態(tài)方程：

式中，ρ為流體密度，kg／m3；u為流體速度，m／s；p為流體壓力，Pa；μ為流體動力黏度，Pa·s；h為焓，kJ／kg；k為流體導熱系數(shù)，W／（m·K）；T為流體溫度，K.

圖1給出的是臨界點附近超臨界水的物性隨溫度的變化情況.由圖可知，4個物性參數(shù)均隨流體溫度發(fā)生著劇烈變化.鑒于此，本文模擬中根據(jù)流體實時溫度分布情況更新其物性參數(shù)信息.

圖1 超臨界水物性隨溫度變化圖Fig.1 Physical properties of SCW in the pseudo-critical zone

本文定義顆粒沿x，y，z方向的極徑分別為a，b，c，顆粒長徑比Ar=a／b=a／c．模擬中固定Re=100，變化不同的墻壁因子（λ=5，7.5，10，20，40）、反應器壁邊界條件（保持靜止或隨流體等速移動）和顆粒長徑比（Ar=0.5，0.75，1，2，2.5），考察不同工況下顆粒在超臨界水中的強制對流特性，這些信息對于揭示超臨界水中顆粒運動和傳熱機理以及完善宏觀耦合模型都起著非常重要的作用．顆粒當量直徑dp固定為0.1 m，顆粒形狀Ar變化時需保證顆粒體積恒為（πd3p）／6.圖2為本文數(shù)值模擬的示意圖，其中計算區(qū)域為圓柱形，長為60 dp，圓柱截面直徑隨λ的改變而不同.

圖2 超臨界水中單顆粒強制對流模擬示意圖Fig.2 Schematic diagram of forced convection of an isolated particle in SCW

圖2 中白色區(qū)域為網(wǎng)格加密區(qū)域，本文對顆粒表面的網(wǎng)格進行加密處理，目的是確保計算的準確性.計算過程中，顆粒中心位于（2 m，2 m，2 m）并保持靜止狀態(tài)，顆粒溫度TH為657 K并保持不變，流體初始溫度TL為647 K，流體初始速度uc由雷諾數(shù)決定.顆粒表面采用第一類溫度邊界.

2 模型驗證

為了驗證模型的準確性，常規(guī)流體（T=300 K，P=0.1 MPa）下，沉降球的壁面效應對Cd和Nu的影響也被加以研究并與文獻中的數(shù)據(jù)進行了對比．表1給出了數(shù)值模擬計算的Cd與文獻［29］結(jié)果對比情況．從表中可以看出，本文模擬計算結(jié)果與文獻報道中的數(shù)據(jù)具有較高的吻合度．因此，用該模型研究墻壁效應對Cd的影響合理.

表1 λ＝5條件下數(shù)值模擬計算的曳力系數(shù)與文獻［29］結(jié)果對比Table 1 Comparison of Cd between the data in reference［29］and current results withλ ＝5

表2給出了數(shù)值模擬計算的Nu與文獻［18］結(jié)果對比情況．從表中可以看出，本文模擬計算結(jié)果與文獻報道中的數(shù)據(jù)高度一致．因此，用本模型研究墻壁效應對Nu的影響可行.

表2 λ＝40條件下數(shù)值模擬計算的平均努賽爾數(shù)與文獻［18］結(jié)果對比Table 2 Comparison of Nu between the data in reference［18］and current results withλ＝40

3 結(jié)果與討論

3.1 截面溫度和壓力系數(shù)分布

圖3展示了單橢球形顆粒在不同λ和Ar下過球心且垂直于流動方向上沉降顆粒與流化顆粒溫度差隨離顆粒表面距離的分布．定義Tf為流化球的截面溫度，而Ts為沉降球的截面溫度.由圖可知，給定一個具體的Ar，對于扁平狀和細長型的橢球顆粒，Ts與Tf的差值在某一具體的位置達到頂峰.這說明在相同條件下，沉降顆粒與流化顆粒的傳熱具有相似性.同時，Ts與Tf的極差值隨著λ的增加而減小，這是因為墻壁效應隨著λ的增加而變?nèi)?

圖3 單橢球形顆粒在不同λ和Ar下過球心且垂直于流動方向上沉降顆粒與流化顆粒溫度差隨離顆粒表面距離分布Fig.3 Distribution of Tf－Ts of an ellipsoid in SCW for different values of Ar andλ（Only one slice perpendicular to the x direction and passing the center of the sphere is displayed）

圖4 描述了超臨界水中單橢球形顆粒不同λ下壓力系數(shù)的分布情況．定義Cd1為流化球問題的壓力系數(shù)，而Cd2為沉降球問題的壓力系數(shù)．從圖中可見，Cd1隨λ的變大而減小，這是因為墻壁效應對Cd的影響隨λ的變大而減弱的緣故．改變λ，Cd2的變化不太明顯，這是因為沉降球問題的墻壁是移動的，墻壁效應對沉降球顆粒在超臨界水中受力影響較大．

圖4 墻壁因子對顆粒表面壓力系數(shù)分布的影響Fig.4 Effects of the wall factor on the distribution of the pressure coefficient on the surface of the particle

3.2 墻壁效應對曳力系數(shù)和平均努賽爾數(shù)的影響

圖5 （a）（b）為超臨界水中單橢球形顆粒Cd隨λ和Ar的變化情況圖．其中，Cdf為流化球的曳力系數(shù)，而Cds為沉降球的曳力系數(shù).由圖可知，Cdf和Cds都隨著Ar的增大而變小，這是因為顆粒的迎風面積隨著Ar的減小而變大，而流體與顆粒之間的曳力與迎風面積呈正相關．除此之外，Cdf和Cds都隨著λ的增大而變小，這是因為隨著λ的減小，壁阻滯效應增強，從而使得Cd增加．但Cdf隨λ的變化趨勢相比Cds更明顯，說明墻壁效應對流化顆粒更強烈.

圖5（c）（d）給出了超臨界水中單橢球形顆粒平均Nu數(shù)隨λ和Ar的變化．其中，Nuf為流化球的平均Nu，而Nus為沉降球的平均Nu．由圖可知，Nuf和Nus都隨著Ar的增大而變小，這是因為顆粒的迎風面積隨著Ar的減小而變大，而流體與顆粒之間的對流傳熱與迎風面積呈正相關．除此之外，Nuf和Nus都隨著λ的增大而變小，這是因為隨著λ的減小，壁阻滯效應增強，從而使得平均Nu增加．但Nuf隨λ的變化趨勢相比Nus更明顯，說明墻壁效應對流化顆粒更強烈.

圖5 單橢球形顆粒在不同λ和Ar下Cd和Nu的分布Fig.5 Distribution of an ellipsoid for different values ofλand Ar

4 新關聯(lián)式的提出

4.1 對前人關聯(lián)式的評估

由于當前的努爾系數(shù)與Kishore等［20］在常規(guī)流體上的數(shù)值結(jié)果具有相似的趨勢，為了使公式適用于超臨界水中，本文對原公式［19－20］進行改進，如式（3）（4）所示.

式中，λ為墻壁因子；Ar為顆粒橫縱比；C1～C7以及a1～a7為常數(shù).為了檢驗Kishore等［19－20］提出的公式是否適用于本文，將目前的數(shù)值模擬結(jié)果與Kishore等［19－20］的公式計算結(jié)果進行比較，如圖5所示.由圖可知，由文獻中的公式計算得到的經(jīng)驗值與實際數(shù)值模擬得到的值存在很大誤差，因此提出新的擬合公式意義重大.

4.2 流化顆粒問題

對于流化顆粒問題，通過對現(xiàn)有數(shù)值結(jié)果的回歸分析，上述式（3）中的新系數(shù)取C1=6.913 48e6，C2=－4.087 82，C3=0.188 42，C4=0.525 11，C5=－1.057 19，C6=－0.171 81，C7=0.641 5，C8=－0.389 94．

從圖7（a）中可以看出，這些新系數(shù)可以很好地預測Cdf.最大和最小相對偏差11.373%和0.0413%，相關系數(shù)為98.99%.

通過對現(xiàn)有數(shù)值結(jié)果的回歸分析，上述式（4）中的新系數(shù)取a1=0.472 91，a2=0.195 68，a3=5.884 13，a4=－0.559 22，a5=46.562 95，a6=0.001 66，a7=－36.619 47.

從圖7（b）中可以看出，這些新系數(shù)可以很好地預測Nuf.最大和最小相對偏差0.2% 和0.002%，相關系數(shù)為99.567%，大多數(shù)數(shù)據(jù)相對偏差遠低于0.1%.

圖6 公式［19－20］預測Cd和Nu與當前數(shù)值模擬結(jié)果對比Fig.6 Predicted by the correlation［19－20］with the present numerical results

圖7 根據(jù)提出的公式預測Cdf和NufFig.7 Predicted Cdf and Nuf by the proposed correlations with numerical results

4.3 沉降顆粒問題

對于沉降顆粒問題，對數(shù)值模擬結(jié)果進行同樣分析，得到Cds方程中的新系數(shù)取C1=4.950 35e8，C2=－4.865 92，C3=0.024 53，C4=2.147 65，C5=－2.591 41，C6=－0.128 43，C7=1.121 21，C8=－0.648 42．

從圖8（a）中可以看出新公式對Cds有很好的擬合效果.最大和最小相對偏差4.581% 和0.105 6%，相關系數(shù)為99.601%.

同樣地，Nus方程中的新系數(shù)取a1=－2.579 33，a2=－2.773 51，a3=3.328 52，a4=－0.787 97，a5=6.098 12，a6=－0.069 3，a7=6.471 26．

從圖8（b）中可以看出，新公式對Nus有很好的擬合效果.最大和最小相對偏差1.212% 和0.013%，相關系數(shù)為99.333%.

圖8 根據(jù)提出的公式預測Cds和NusFig.8 Predicted Cds and Nus by the proposed correlations with numerical results

5 結(jié) 論

（1）在兩種顆粒擾流問題中，Cd和Nu均隨著Ar及λ的增大而減?。捎趬Ρ谛?，在同一Ar和λ條件下，流化顆粒對應的Cd值較大．

（2）超臨界水中的曳力系數(shù)與平均努賽爾數(shù)和前人的經(jīng)驗公式計算值存在很大差距．因此，在0.5≤Ar≤2.5和5≤λ≤40的條件下，本文分別建立了兩種墻壁邊界條件下的關于Cd和Nu的關聯(lián)式，可在CFD-DEM 等宏觀多相流模擬中采用.