邊界條件對橢球顆粒在超臨界水中受力與傳熱特性影響的數值模擬

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

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

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

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

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

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

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

1 數學模型

本文采用的是三維穩態方程：

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

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

圖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為本文數值模擬的示意圖，其中計算區域為圓柱形，長為60 dp，圓柱截面直徑隨λ的改變而不同.

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

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

2 模型驗證

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

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

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

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

3 結果與討論

3.1 截面溫度和壓力系數分布

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

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

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

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

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

圖5（c）（d）給出了超臨界水中單橢球形顆粒平均Nu數隨λ和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 新關聯式的提出

4.1 對前人關聯式的評估

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

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

4.2 流化顆粒問題

對于流化顆粒問題，通過對現有數值結果的回歸分析，上述式（3）中的新系數取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）中可以看出，這些新系數可以很好地預測Cdf.最大和最小相對偏差11.373%和0.0413%，相關系數為98.99%.

通過對現有數值結果的回歸分析，上述式（4）中的新系數取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）中可以看出，這些新系數可以很好地預測Nuf.最大和最小相對偏差0.2% 和0.002%，相關系數為99.567%，大多數數據相對偏差遠低于0.1%.

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

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

4.3 沉降顆粒問題

對于沉降顆粒問題，對數值模擬結果進行同樣分析，得到Cds方程中的新系數取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%，相關系數為99.601%.

同樣地，Nus方程中的新系數取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%，相關系數為99.333%.

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

5 結 論

（1）在兩種顆粒擾流問題中，Cd和Nu均隨著Ar及λ的增大而減小．由于墻壁效應，在同一Ar和λ條件下，流化顆粒對應的Cd值較大．

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