太陽能熱化學反應器內流動和反應的數值模擬

翟占虎，李 巍，劉國棟

(1.哈爾濱工業大學能源科學與工程學院，黑龍江 哈爾濱 150001；2.海洋石油工程股份有限公司，天津 300452)

太陽能因其儲量豐富、清潔可再生等特性，成為當前研究的熱點問題，但是由于其熱流密度低、具有周期性且不易儲存，目前還只是作為輔助能源.近年來太陽能聚焦技術取得長足發展，將光熱轉換和高溫熱化學反應相耦合的光化學利用技術引起重視，該技術以氧載體材料為依托，在高溫下將 CO2、H2O 等熱分解產生H2和CO等碳氫燃料.

圍繞著太陽能反應器，國內外學者開展大量的研究工作.瑞士PSI太陽能技術實驗室Steinfeld[1]教授領導的團隊設計并試驗的太陽能熱化學反應器，能夠將CH4-ZnO在高溫下分解為Zn和合成氣.Chueh等[2]設計的太陽能反應器，參考植物將太陽能轉化為燃料的過程，研究了CeO2在熱化學循環中生產燃料的熱力學和動力學性質以及反應的穩定性，在不進行廢熱回收的前提下，該反應的熱化學學循環效率理論上可達到16%～19%.Meier等[3]設計并實驗的太陽能間接加熱回轉窯，可以在1 200 K～1 400 K溫度區間內實現高溫煅燒反應，并獲得高純度石灰.吳璽[4]針對ZnO與CH4-Ar的重整制氫反應，數值模擬研究了吸收腔內溫度場分布情況，并分析了腔體溫度、反應物料進口速度和CH4-Ar體積比對重整制氫反應效率的影響.王富強等[5]采用數值模擬對比研究了輻射模型多孔介質太陽能接收器傳熱性能的影響.黃興等[6]研究了催化劑NiFe2O4的熱解機理并應用于太陽能熱化學反應器，采用Ansys Fluent研究了反應器內速度場和溫度場，分析了顆粒直徑和氣體流速對反應器內部的溫度分布影響.Lougou等[7]數值模擬研究了太陽能熱化學反應器作用下的甲烷重整反應，分析了NiFe2O4催化作用下，混合氣體流速、操作壓力、催化劑裝載量和輻射溫度分布等參數對反應轉化率的影響.張浩等[8]采用數值模擬結合反應動力學，研究了太陽能反應器中的甲烷重整反應和積碳特性，得到多孔介質的最佳空隙率和顆粒直徑.Bader等[9]設計了太陽能高溫熱化學反應器，以二氧化鈰(CeO2)為固定床催化劑，利用該金屬氧化還原作用的分解CO2，該反應器中催化材料與吸收腔不直接接觸，CeO2相對于Fe3O4和ZnO等催化材料，具有儲放氧性能高、比表面積大、傳熱傳質速率高和結構穩定等優點[10]，避免了催化劑的污染和分離回收，同時獨立的反應器提升了系統穩定性和操作性，適用于存在一定溫度波動的高溫反應器.

太陽能熱化學反應器反應效率受熱流密度分布和金屬氧化物催化劑的氧化還原過程的影響，實驗直接獲得相關的參數比較困難，因此本文數值模擬研究Bader所設計的太陽能熱化學反應器，反應過程中氧空位濃度δ的變化通過求解編寫的UDS(用戶自定義標量方程)實現，質量和能量變化通過UDF(用戶自定義函數)加載到多孔介質區域，最后分析氧化還原反應階段的溫度、氧空位以及產物濃度分布情況.

1 反應器結構與數學模型

1.1 反應器結構及物理模型

反應器整體結構如圖1所示，由一石英窗、光圈、保溫層、圓筒形吸收腔及腔內環形排列著的七個反應器構成，吸收腔內部長為347 mm，半徑為152.5 mm，每個反應器由同軸圓管構成，反應器外管的外/內半徑分別為34.9 mm/31.8 mm，內管的外/內半徑分別為22.2 mm/19.1 mm，圓筒壁面厚為3.1 mm，兩圓筒夾層厚為9.5 mm，夾層中放置718 g固定的金屬氧化物催化劑二氧化鈰(CeO2)，該催化劑為顆粒直徑為5 mm，催化劑顆粒的孔隙率為70%，固定床層的空隙率為40%.

圖1 反應器幾何結構圖

反應器工作流程如下，經過聚焦后的高通量太陽輻射能(3 kW)，透過光圈射在吸收腔壁面和七個圓管反應器的表面，吸收了太陽輻射能的氧化鋁管壁將熱量傳入固定床層的二氧化鈰，催化材料二氧化鈰(CeO2)其表面及體相的晶格氧原子能夠直接參與反應并被消耗[11]，同時形成氧空位(δ)，產出化學燃料.在反應器工作過程中，周期性循環地向反應管內通入N2和CO2，通入N2時長為100 s，質量流量為1.1 g·s-1，該階段發生還原反應，出口產物有N2和O2；通入CO2時長為100 s，質量流量為1.7 g·s-1，該階段發生氧化反應，出口產物有CO2和CO.反應方程式如下所示：

還原反應(Reduction)

(1)

氧化反應(Oxidation)

CeO2-δrd+(δrd-δox)CO2=CeO2-δox+(δrd-δox)CO

(2)

二氧化鈰顆粒及其所構成的固定床，在模擬中被簡化為多孔介質，根據相關公式計算出的輸運性質如表1所示.

表1 二氧化鈰和固定床層輸運性質

為簡化計算，取模型的七分之一進行模擬，首先模擬得到吸收腔內熱流密度分布，然后將熱流密度加載到單個反應器，再進行反應器內流動、傳熱和化學反應模擬.

1.2 控制方程

反應過程中滿足連續性方程、動量方程、組分輸運方程、能量守恒方程.

連續性方程

(1)

公式中：ρ、u、Ri分別為密度、速度和質量源項；下標f、s、i分別為流體、固體和組分i.

動量守恒方程

(2)

公式中：φ、p、μ、K、CE分別為空隙率、壓力、動力粘度、滲透率和慣性阻力系數.

組分輸運方程

(3)

公式中：φeff、Yi分別為有效空隙率和組分i的質量分數.

能量守恒方程

(4)

(5)

(6)

DO輻射方程：

(7)

由于反應器吸收腔內溫度較高，需要考慮其內部的輻射特性，在模擬吸收腔內部熱流密度分布時采用DO輻射模型.

1.3 化學反應動力學方程

目前的文獻中，沒有準確描述二氧化鈰催化分解CO2的反應動力學表達式，本研究中采用的動力學表達式主要參考Bulfie和Rohini Bala的文章[14-15].在二氧化鈰參與氧化還原反應的過程中，其非化學計量數δ決定反應進行速率的大小，以下為δ隨時間變化的微分方程，以及各組分質量源表達式.

δ與t的關系：

(8)

(9)

在氧化與還原過程中，CeO2既能釋放O2，也能吸收O2，因此氧氣質量源項的表達式如下：

(10)

研究表明[17]，高溫下CO2的直接分解作用也很明顯，且為可逆反應，各組分質量源項的表達式如下：

(11)

公式中：i=O2，CO，CO2，相應地vi=0.5，1，-1.

O2的凈產量為CO2分解所釋放的量與CeO2氧化消耗的量之和，如下所示：

RO2，net=RO2，diss+RO2，ceria

(12)

預計熱分解產物CO會在短時間內和O2重新結合成CO2，故本研究認為氧化階段產生的O2，一部分被CeO2吸收來填補氧空位，其余的全部與CO重新結合.

本研究以上述的非化學計量數δ的微分方程、各組分的質量源項表達式為理論基礎，來實現上述化學反應.

2 模擬結果分析

Bader[9]文中給出的太陽能熱化學反應器為概念設計，并未進行實驗測量，反應器在設計運行參數下，工作溫度約為1 800 K，CO產率在預估為3.9×10-4mol·s-1，模擬得到的溫度分布如圖4和圖11所示，以及出口CO摩爾產率如圖17所示，較為接近，可認為模擬所采用模型是合理的.

2.1 吸收腔內熱流密度分布

吸收腔內部熱流密度分布情況，如圖2所示.可以看出圓管反應器壁面的熱流量明顯高于吸收腔壁面，經計算圓管反應器壁面的平均熱流密度為4 600 W·m-2，總熱流量為318.5 W，吸收腔壁面的熱流密度為181.5 W，熱量損失占36.3%.沿著圓管反應器軸的熱流密度分布曲線，如圖3所示.沿著軸向先增高后減少，并且在x=0.12 m達到最大值6 900 W·m-2.

圖2 反應器熱流密度分布云圖圖3 沿X軸向熱流密度

2.2 還原階段分析

(1)反應器溫度分析

t=1 s和100 s時反應器溫度分布云圖，如圖4所示.從左側圖可看出反應器壁面溫度最高為1 880 K，t=100 s時超過1 860 K溫度的范圍明顯少于t=1 s時的，這是因為CeO2的還原反應是吸熱反應，隨著反應的經行，反應器溫度下降.

圖4 t=1 s和100 s時反應器溫度分布

反應器壁面軸向溫度在不同時刻的溫度分布和多孔介質的體積平均溫度隨時間的變化，如圖5、圖6所示.壁面溫度和體積平均溫度都隨著反應的進行，不斷減小，進一步驗證了CeO2的還原反應是吸熱反應.

圖5 反應器壁面軸向溫度在不同時刻的溫度分布圖6 多孔介質體積平均溫度隨時間的變化

(2)氧空位濃度δ與O2產率

t=1 s時和t=100 s多孔介質區域的δ分布云圖，如圖7所示.在t=100 s時δ最高到0.041，且最高位置在出現在剛進入多孔介質區域的地方.這是因為氣體剛進入多孔介質時，由于不含氧氣，故此處氧氣的分壓為零，遠低于平衡壓力，CeO2中的氧原子在高溫條件下能夠以氧氣的形式釋放出來，故此處的氧空位最高.從圖7中還可看出，靠近出口處多孔介質的δ最小，這是因為氧氣的不斷釋放，氧氣在下游的質量分數達到最高，即氧氣分壓也達到最高，故此處的CeO2中的氧原子受氧氣分壓的抑制，釋放速率較慢.

圖7 不同時刻δ在多孔介質區域的分布云圖

沿軸向δ是不斷減小的，且隨著時間增加，δ的增長變緩，多孔介質上游區域的δ增長率明顯高于下游，如圖8所示.由圖9可知，氧氣的質量分數沿軸向不斷變大，隨之時間增加，氧氣的質量分數不斷降低.這是因為當CeO2剛開始反應時，由于多孔介質區域各處的氧氣分壓均為零，O2的產率很高，隨著反應的進行，氧氣分壓升高，特別是多孔介質下游區域，其氧氣分壓明顯高于上游，更接近平衡壓力，當氧氣分壓越接近平衡壓力，δ增長越慢，氧氣的產率越小，相反地上游區域氧氣分壓始終較低，δ增長始終較快，所以x=0 m處的δ值明顯高于x=0.374 m處的.

圖8 軸向氧空位濃度δ分布(不同時刻)圖9 軸向氧氣質量分數分布(不同時刻)

出口氧氣的摩爾流率，如圖10所示.在剛通入N2的時候，出口氧氣的摩爾流率可達到9.42×10-5mol·s-1，且隨著反應時間的增加，氧氣流率逐漸下降，這是因為總體上氧氣的分壓越來越接近平衡壓力(大約為240 Pa)，δ增長減緩，從CeO2中釋放出的氧氣變少.

圖10 出口氧氣的摩爾流率

2.3 氧化階段分析

(1)反應器內溫度分布

t=101 s和200 s時反應器溫度分布云圖，如圖11所示.在整個時間段內，反應器的最高溫度為1 770 K，低于還原階段的1 880 K，而從101 s到200 s這段時間，反應器內高溫區域變大.這是因為通入CO2后，CO2質量流量比N2的質量流量大，且CO2氣體的比熱容大于N2的比熱容，故在一定的熱流密度下，氧化階段的流場溫度低于還原階段時的溫度.此外，在整個氧化階段，反應是放熱的，所以隨著時間的增加，反應器整體升溫.

圖11 不同時刻壁面和計算域內溫度分布云圖

反應器壁面軸向溫度在不同時刻的溫度分布，如圖12所示.從圖12可知，隨著軸向距離的增加，壁面溫度先增高后下降.因為氣體被壁面不斷加熱，隨著流動距離溫度變高，在臨近出口時，由于與入口低溫氣體(1 625 K)換熱，所以溫度出現了下降，由于反應放熱，所以隨著時間增加，壁面溫度升高.從圖13可知，隨著時間的增加，多孔介質體積平均溫度不斷增加，進一步說明了氧化階段發生的反應是放熱的.

圖12 反應器壁面軸向溫度分布圖13 多孔介質體積平均溫度隨時間的變化

(2)氧空位濃度δ與CO產率

從圖14可看出，隨著軸向距離增加，氧空位濃度δ逐漸降低，且隨著時間的增加，δ也在減少，特別是x=0 m處，δ下降最快，在t=200 s時多孔介質區域δ的分布云圖，如圖15所示.可見在x=0 m處，δ約為0.03，體積平均δ約為0.020，與還原階段的初始值較接近，這是因為氧化階段，CO2高溫下發生熱分解，釋放出大量O2，使得此時的氧氣分壓大于平衡壓力，CeO2將O2吸收來填補自身的氧空位，所以δ出現下降趨勢，且δ越大，氧氣平衡壓力就越小，CeO2對O2的吸收速率更高，因此在x=0 m處δ下降最快.

圖14 軸向氧空位濃度δ分布情況(不同時刻)圖15 t=200s時δ分布云圖

從圖16可看出，隨著x的增加，CO質量分數先增加后減少，且隨著時間的增加，CO質量分數不斷下降，這是因為當通入CO2時，CO2熱分解釋放的O2被二氧化鈰CeO2吸收用以填補氧空位，進而促進CO2的熱分解，所以從x=0 m到x=0.15 m處，CO質量分數不斷升高；而CO質量分數在x=0.15 m之后出現下降趨勢，是因為這段區域的氧空位濃度δ較小，對O2的吸收能力有限，導致CO2熱分解產生的O2沒有被完全吸收，與CO再結合，導致CO質量分數下降.

CO摩爾流率隨時間的變化曲線，如圖17所示.隨時間增加，CO摩爾流率逐漸下降，這是因為隨著反應的進行，氧空位濃度δ逐漸減少，對CO2熱分解反應的促進表少，導致CO產率逐漸減少.

圖16 軸向CO質量分數分布 圖17 出口CO的摩爾流率

3 結 論

本文以Bader設計的反應器為基礎，建立了用于二氧化鈰小溫差下氧化還原循環分解CO2的太陽能高溫反應器的三維瞬態模型，考慮了流動過程和非化學計量反應的相互作用及其對溫度、壓力、組分濃度和燃料產率的影響.研究發現，反應器熱流密度集中分布在靠近光圈的一側，熱量損失取決于吸收腔壁面溫度和材料保溫性能，應盡量減少輻射能量直接作用于吸收腔壁面；反應器溫度與反應階段密切相關，還原階段溫度上升，氧化階段溫度下降；還原階段，二氧化鈰釋放氧氣，氧空位濃度δ受氧氣的抑制作用，在氣體入口處最大，出口處最小；氧化階段，二氧化碳在高溫下熱分解，二氧化鈰吸收分解后的氧原子，促進了CO2的熱分解作用；氧氣釋放量和二氧化碳熱分解作用與氧空位濃度、氣體流率、反應時間和溫度密切相關，在太陽能反應器設計階段需綜合考慮上述因素.