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

（東北大學(xué)1.冶金學(xué)院；2.教育部多金屬礦山生態(tài)冶金重點(diǎn)實(shí)驗(yàn)室，沈陽(yáng) 110819）

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

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

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

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

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

1 數(shù)學(xué)模型

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

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

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

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

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

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

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

2 模型驗(yàn)證

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

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

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

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

3 結(jié)果與討論

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

圖3展示了單橢球形顆粒在不同λ和Ar下過(guò)球心且垂直于流動(dòng)方向上沉降顆粒與流化顆粒溫度差隨離顆粒表面距離的分布．定義Tf為流化球的截面溫度，而Ts為沉降球的截面溫度.由圖可知，給定一個(gè)具體的Ar，對(duì)于扁平狀和細(xì)長(zhǎng)型的橢球顆粒，Ts與Tf的差值在某一具體的位置達(dá)到頂峰.這說(shuō)明在相同條件下，沉降顆粒與流化顆粒的傳熱具有相似性.同時(shí)，Ts與Tf的極差值隨著λ的增加而減小，這是因?yàn)閴Ρ谛?yīng)隨著λ的增加而變?nèi)?

圖3 單橢球形顆粒在不同λ和Ar下過(guò)球心且垂直于流動(dòng)方向上沉降顆粒與流化顆粒溫度差隨離顆粒表面距離分布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為流化球問(wèn)題的壓力系數(shù)，而Cd2為沉降球問(wèn)題的壓力系數(shù)．從圖中可見(jiàn)，Cd1隨λ的變大而減小，這是因?yàn)閴Ρ谛?yīng)對(duì)Cd的影響隨λ的變大而減弱的緣故．改變?chǔ)耍珻d2的變化不太明顯，這是因?yàn)槌两登騿?wèn)題的墻壁是移動(dòng)的，墻壁效應(yīng)對(duì)沉降球顆粒在超臨界水中受力影響較大．

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

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

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

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

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

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

4.1 對(duì)前人關(guān)聯(lián)式的評(píng)估

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

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

4.2 流化顆粒問(wèn)題

對(duì)于流化顆粒問(wèn)題，通過(guò)對(duì)現(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ù)可以很好地預(yù)測(cè)Cdf.最大和最小相對(duì)偏差11.373%和0.0413%，相關(guān)系數(shù)為98.99%.

通過(guò)對(duì)現(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ù)可以很好地預(yù)測(cè)Nuf.最大和最小相對(duì)偏差0.2% 和0.002%，相關(guān)系數(shù)為99.567%，大多數(shù)數(shù)據(jù)相對(duì)偏差遠(yuǎn)低于0.1%.

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

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

4.3 沉降顆粒問(wèn)題

對(duì)于沉降顆粒問(wèn)題，對(duì)數(shù)值模擬結(jié)果進(jìn)行同樣分析，得到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）中可以看出新公式對(duì)Cds有很好的擬合效果.最大和最小相對(duì)偏差4.581% 和0.105 6%，相關(guān)系數(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）中可以看出，新公式對(duì)Nus有很好的擬合效果.最大和最小相對(duì)偏差1.212% 和0.013%，相關(guān)系數(shù)為99.333%.

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

5 結(jié) 論

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

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