基于CFD模擬的甲烷裂解太陽能管式反應器結構優化

肖凡，賈勝坤，羅祎青，袁希鋼，2

（1天津大學化工學院，天津 300354；2化學工程聯合國家重點實驗室（天津大學），天津 300354）

引 言

氫能作為21世紀最具有發展潛力的清潔能源之一受到普遍關注[1]。目前制氫方法主要有電解水法[2-4]、煤氣化法[5-6]、甲醇裂解法[7]等。電解水制氫能量消耗大；煤氣化法制氫會產生大量二氧化碳，加重溫室效應，不利于“節能減排”戰略的實施；甲醇裂解制氫，產物不純，并伴有大量的碳氧化物，后續分離工藝煩瑣，能耗大。而化石燃料（低碳烷烴）經高溫直接裂解成氫氣和炭黑不產生二氧化碳，近年來受到研究者的廣泛關注。太陽輻射清潔易得，輻射集中后產生的高溫可以作為裂解制氫的熱能來源。因此，太陽能高溫裂解甲烷制氫研究具有重要的意義和發展前景[8]。

太陽能高溫裂解甲烷過程，主要是利用太陽光聚焦產生的高溫熱能加熱反應器，促使在氬氣保護下的甲烷直接發生高溫裂解反應，生成氫氣和炭黑[9]。太陽能的利用減少了化石能源的消耗，副反應較少，并且生成的炭黑是重要的化工原材料，因此優勢顯著。目前，高溫裂解甲烷反應器主要有顆粒流反應器[10]、渦流反應器[11]和管式反應器[12]等形式。其中，管式反應器主要包括單管式[13]和多管式反應器[14]，因其結構簡單而受到關注。在管式反應器的研究中，發現影響甲烷轉化率和氫氣收率的不僅包括進料氣體流速及甲烷濃度、反應溫度、催化劑，還包括反應器管徑、入口結構[15]以及太陽能加熱功率[16]等因素。與實驗研究對比，計算流體力學（CFD）模擬研究具有成本低和容易實施等優點，越來越多地被用于探究各種因素對管式反應器甲烷轉化率以及氫氣產率的影響[17]，且CFD模擬可對實際實驗進行指導[18]。

在已有研究中，為保證反應器內較低的流速和長停留時間，反應器內的流體大多為層流狀態[15-19]，以實現較高甚至百分之百的轉化率。但層流反應器難以滿足實際工業的高通量需求，因此需要提高反應物的流速來提高反應通量。較大的流速會導致反應物停留時間更短，而受到太陽能反應器透鏡尺寸等因素的限制，很難通過增加反應管的長度來提高停留時間。因此，為了保證甲烷轉化率，在給定反應器管長內，通過改變流場實現反應過程的強化具有重要意義。Cao等[20]在湍流狀態下，以反應過程的熵產最小為目標函數，通過虛擬外力對太陽能甲烷裂解反應器內的流場進行了優化，提高了湍流狀態下的甲烷轉化率，但如何將虛擬外力實物化尚未解決。

本文針對太陽能甲烷裂解反應器，通過CFD模擬的方法，開展通過主動和被動方式實現流場的優化研究，提出反應器具體建構方案，并針對實際反應器的優化設計給出定量的參考。本文首先建立了甲烷裂解反應器的湍流CFD模型，在湍流反應擴散模型的基礎上，引入了碳顆粒的生成和聚集模型，實現對碳顆粒輻射吸收系數的更準確計算，利用離散坐標（DO）模型[19]對太陽能輻射場進行求解，并通過已有的實驗結果對模擬結果進行了驗證。然后，通過引入翅片擋板[21-22]以及射流[23]措施對流場進行優化調節，對射流流速及角度、翅片擋板高度進行優化，實現反應器內流場的優化，從而強化反應過程[24]，提高甲烷轉化率。最后，以增加的反應轉化率和代表強化成本的黏性耗散為指標，篩選出離散的Pareto最優解，并利用機器學習中的支持向量機回歸（SVR）算法[25]對離散的Pareto最優解進行插值，得到連續的最優操作曲線以及與之相對應的射流角度和流速等操作參數。使用SVR算法進行插值相對于傳統插值方法具有預測誤差小、對原始數據分布要求低、可以在數據量較小時取得較好的擬合效果等優點。

1 考慮炭顆粒及其輻射的CFD模型

1.1 流體流動控制方程

假設顆粒流速和氣體主體流速保持一致，且對氣體的物理性質沒有影響，故此反應可視為均相反應。流動模型采用基于Navier-Stokes方程的均相模型，湍流模型采用標準k-ε模型進行封閉。

流體流動控制方程主要分為連續性方程、動量傳遞方程以及能量傳遞方程。二維均相流體流動模型：

式中，μe為有效黏度，為氣體混合物的黏度與湍流黏度之和，混合物黏度采用質量平均組分黏度，湍流黏度采用“標準k-ε”模型進行計算。

1.2 能量方程

甲烷熱裂解反應在高溫環境下具有強吸熱效應，其標準狀態下的反應焓變為75 kJ·mol-1[26-27]。因此，在CFD模擬中必須考慮能量方程，即：

式中，hi是組分i的焓變；Schem是反應熱；Srad是輻射熱。

沿路程S輻射傳遞方程為：

沿路程S影響輻射強度[28]大小的量包含四部分：吸收引起的衰減、發射引起的增強、散射引起的減弱以及射入散射引起的增強。式（5）采用DO模型進行求解。

1.3 組分輸運方程

太陽能高溫熱裂解甲烷反應過程，涉及到甲烷、氫氣、氬氣、炭顆粒等幾種組分。為了更好地模擬炭顆粒粒徑對輻射系數的影響，將炭顆粒粒徑劃分為不同的尺寸。反應組分（甲烷、氫氣、炭顆粒）的傳遞通過組分輸運方程進行描述：

式中，Ri為反應源項；下角標i=1~5表示5種不同尺寸的炭顆粒組分，i=6和7則代表CH4和H2;Di,e為有效擴散系數：

式中，Sct為湍流Schmidt數，取為0.7，即假設湍流擴散系數與湍流黏度系數成正比。

1.4 反應源項

反應源項包含甲烷、氫氣以及炭顆粒等幾種組分。甲烷熱裂解反應過程中，伴隨著氫氣和亞微米級炭黑顆粒生成，這些炭黑顆粒在氣體混合物中懸浮，對流體流動、傳熱以及反應過程都會產生影響。實際反應中，會有少量乙烷、乙烯、乙炔等副產物生成[22,24]，但含量很小，在本文的模擬過程中予以忽略。

甲烷太陽能熱裂解反應方程式可表示為：

無催化劑條件下，反應動力學方程[18,26-27]為：

式中，α=0.84，Ea=259 kJ·mol-1，K0=7.68×108s-1。

甲烷直接裂解生成氫氣和炭黑顆粒，組分輸運方程中炭顆粒組分的源項[18]Ri為：

炭顆粒分為生成與聚集兩部分；組分運輸方程中甲烷和氫氣組分的源項分別為：

1.5 炭顆粒生成與聚集模型

在實際的甲烷裂解反應中，炭顆粒的粒徑主要分布在1~9μm[18-19]。因此，考慮計算的精度與負荷的平衡，將炭顆粒組分分為5種不同的尺寸等級，對應的尺寸分別為1、3、5、7、9μm。第i等級炭顆粒的質量分數為：

式中，mi是第i階單個炭顆粒的質量；Ni是單位體積內i階炭顆粒的個數。

炭顆粒會通過反應不斷生成，且顆粒間也會發生聚集合并。不同尺寸的炭顆粒來源包括均相反應[N1]hom、非均相反應[Ni]het[29-30]以及顆粒的聚集[Ni][31-33]coag三部分。均相反應（成核反應）速率采用Arrhenius反應速率表達形式：

非均相反應主要發生在炭顆粒的表面（炭顆粒成長反應）：

式中，Hhet,i為甲烷擴散到級別i炭顆粒上的傳質系數，具體計算公式為：

均相反應和非均相反應會形成粒徑較小的初級粒子，這些初級粒子形成后會經過布朗碰撞發生聚集生成新的粒子。布朗碰撞采用菲克聚集系數模型[18-19]。其中，初級聚集單體的來源為：

多級聚集單體的來源為：

湍流聚集系數采用菲克聚集系數：

炭顆粒的平均輻射吸收系數：

式中，Qabs指的是炭顆粒的平均光譜吸收率。本節計算得到的輻射系數將用在式（5）中，利用DO模型來求解輻射傳熱。

2 反應邊界條件與模型驗證

2.1 幾何模型和邊界條件

太陽能反應器幾何模型如圖1所示。反應管內徑為10 mm，長度為100 mm。上下壁面溫度為2100 K，壁面輻射發射率設為0.9，與高溫下的石墨發射率保持一致。甲烷和氬氣以摩爾比1∶1的比例混合后均勻進入反應器，氣體流速為2 m·s-1，入口氣體溫度1800 K。采用充分發展的出口邊界條件。

圖1 太陽反應器單元模型Fig.1 Schematic diagram for a cell of solar reactor

2.2 數值模擬求解

選用ANSYS FLUENT 16.0TM商業軟件求解器對控制方程進行求解。其中，壓力速度耦合采用SIMPLE Scheme進行求解，梯度選擇Least Squares Cell Based進行求解，壓力選擇Standard進行求解。連續性方程（1）、動量方程（2）、能量傳遞方程（3）、組分輸運方程（6）選擇QUICK Scheme進行求解。輻射傳遞方程（5）選擇二階迎風格式（second order upwind scheme）進行求解。CH4、H2以及炭顆粒組分選擇一階迎風格式（first order upwind scheme）求解。DO模型求解時每個八分圓中的天頂角和方位角采用4×4進行離散化求解。產生10套不同密度的計算網格，在相同的工況條件下進行計算，以反應管出口溫度變化率不超過1%為標準進行網格獨立性驗證，考慮計算負荷與計算精度，選擇10萬個結構化網格數進行模擬。

模擬過程選擇穩態模擬，甲烷轉化率計算公式如下：

式中，min為進口甲烷質量；mout為出口甲烷質量。

2.3 模型驗證

為驗證模型的可靠性，本文以Abanades等[14]的實驗為基礎，以該實驗太陽能反應器單元管為原型進行CFD模擬。參照文獻[14]實驗條件進行邊界條件設置，控制方程采用式（1）~式（25）。模擬結果如表1所示，實驗值與模擬值吻合很好，證明了本文模型準確。

表1 模型驗證數據Table 1 Model validation data

3 太陽能管式反應器甲烷熱裂解流場優化

3.1 反應管內結構優化

Cao等[20]以熵產極值為目標，獲得了湍流條件下太陽能高溫甲烷裂解制氫管式反應器中最優流速場，并證明該流速場可作為該反應過程優化的熱力學極限。該最優流速場主要特征是，在靠近管的中間位置有一個很接近管徑尺度的非軸對稱渦，這種流型實現了在給定黏性耗散（即流體阻力）條件下系統的熵產速率最小，即熱力學效率最高[20]。本文根據文獻[20]給出的最優流場結構特征，通過在管壁增加射流以及翅片擋板的方式重現該最優流場結構。本文采取的射流和翅片擋板等措施如圖2所示。上壁面距離進口33.5~35.5 mm位置處設置射流入口，上壁面52.5 mm和下壁面60 mm位置處設置翅片擋板。上述結構尺寸通過多次試算得到，其目的是獲得盡量靠近文獻[20]得到的流場。

圖2 太陽能反應器單元結構優化Fig.2 Optimization schematic diagram for a cell of solar reactor

3.2 優化結果與討論

本研究對射流流速V、擋板高度H以及射流角度θ三個參數進行優化。

3.2.1 射流流速優化 在得到基本優化結構后，本文固定擋板高度為1.5 mm，射流角度為0°，射流氣體組分為氬氣，溫度為2100 K，僅對射流速度進行優化。射流速度變化范圍是1~15 m·s-1。以射流流速9 m·s-1為例進行說明，優化前、后的流線、溫度、甲烷質量分數對比分別如圖3~圖5所示。

比較圖3(a)、(b)可以看出，射流的引入導致流場發生擾動，流線圖與文獻[20]中給出的最優流場接近。比較圖4(a)、(b)可以看出，在反應管中段，經過射流擾動后的溫度明顯高于優化前。甲烷裂解為吸熱反應，溫度升高有利于甲烷裂解反應的進行，當射流流速為9 m·s-1時，甲烷轉化率可以提升約2%左右。比較圖5(a)、(b)可以看出，由于前半段甲烷裂解反應更為充分，優化后的甲烷濃度在反應管后半段更低。

圖3 無優化(a)和優化后(b)太陽能反應器內的流線圖Fig.3 Streamlines in the solar reactor for the base case(a)and optimization case(b)

圖4 無優化(a)和優化后(b)太陽能反應器內的溫度分布Fig.4 Temperature profile in the solar reactor for the base case(a)and optimization case(b)

圖5 無優化(a)和優化后(b)太陽能反應器內的甲烷質量分數分布Fig.5 CH4 mass fraction profile in the solar reactor for the base case(a)and optimization case(b)

圖6(a)~(d)分別為反應器中x=34.5、52.5、60和90 mm位置處，不同射流速度下的溫度分布對比。圖7(a)~(d)是相同位置處的甲烷摩爾分數分布。由圖6能夠明顯觀測到，在y=3~6 mm范圍內優化后的溫度有明顯的增加，且隨著流速增大，溫度也隨之升高。優化后反應器的中后部溫度升高，使得化學反應推動力增大，引起甲烷轉化率的增加。

對比圖6和圖7可以發現，甲烷的摩爾分數分布與對應位置的溫度分布大致呈負相關趨勢。隨著射流速度的增加，在y=2~8 mm間，甲烷摩爾分數隨之降低，甲烷熱裂解反應充分進行。

圖6 無優化（直管）和優化后反應器內的溫度分布對比Fig.6 Comparison of temperature distribution for the base case(straight tube)and the optimization case in the reactor

圖7 無優化（直管）和優化后反應器內的甲烷摩爾分數分布對比Fig.7 Comparison of CH4 mole fraction distribution for the base case(straight tube)and the optimization case in the reactor

3.2.2 擋板高度優化 本文保持管內反應氣體組分比例與質量通量不變，在不同射流流速下，對擋板高度進行優化。射流速度變化范圍是1~9 m·s-1，擋板高度分別取1.0、1.5和2.0 mm，結果如表2所示。模擬發現，在確定的氣體流速下，擋板高度達到1.5 mm后，增加擋板高度對甲烷轉化率影響不明顯。

表2 不同擋板高度在不同射流流速下甲烷增加的轉化率Table 2 The increase of methane conversion rate for different baffle heights at different jet flow rates

3.2.3 射流角度優化 射流對反應器流場分布有顯著影響，而擋板高度的影響甚微，因此本節對射流角度進行優化，同時擋板高度不變，設定為1.5 mm。y方向分量射流速度為1~9 m·s-1，考慮到工業氣體流速極限，實際氣體流速不宜超過40 m·s-1。設θ為射流方向與豎直方向的夾角，順時針旋轉為正，射流角度變化范圍取-30°~85°，以5°為間隔變換角度考察計算結果。

圖8給出不同射流速度和角度下甲烷轉化率增加量的計算結果。可以看出，正角度的射流會提高甲烷轉化率，而負角度射流會降低甲烷轉化率。其原因應為正角度射流進入反應管內產生逆向擾動，流體湍動程度更為劇烈，反應更加充分，甲烷轉化率會有所提高；在負角度下，射流僅會加速反應氣體流出，反應停留時間減少，導致反應不充分。圖8表明，在相同角度下，隨著射流流速的增加，甲烷轉化率呈現先上升后趨于平緩且下降趨勢。這表明確定角度下射流流速并不是越大越好，存在最佳的流速，且這一最佳射流流速隨射流角度的增加而增加。圖8還表明，射流流速與角度對甲烷轉化率增加量的影響呈現較為復雜的曲線交叉現象，在交叉區域，射流角度增加的同時射流流速需要同時增加才能提高甲烷轉化率。例如，當射流速度為10 m·s-1時，射流角度80°給出的甲烷轉化率增加量反而低于65°、70°和75°對應的甲烷轉化率增加量。

圖8 不同角度及不同射流流速下甲烷增加的轉化率（相同流速下，以1.5 mm擋板、0°夾角結構下甲烷轉化率為基準進行比較）Fig.8 Increase of methane conversion rate in different angles and different jet flow rates

圖9(a)~(c)分別為反應器中x=52.5、60和90 mm位置處射流角度優化前后溫度分布對比（射流流速5 m·s-1）。圖9表明，射流角度的變化會影響反應管內的溫度分布。在y=2~8 mm的管內中心區域，溫度隨著射流角度的增加而逐漸升高。射流角度的提升會加大流場的擾動，導致高溫氣體在管內混合較為充分，管內溫度中半段會有所提升，甲烷裂解為吸熱反應，溫度升高會加快反應進行，相同時間下，反應會進行得更加徹底，甲烷轉化率會有所提高。

圖9 流速5 m·s-1下無優化（0°）和優化后（30°,45°,60°）反應器內的溫度分布對比Fig.9 Comparison of temperature distribution in the reactor without optimization(0°)and after optimization(30°,45°,60°)at a flow rate of 5 m·s-1

3.3 最優操作參數的確定

3.2節結果表明射流流速及角度的改變、擋板的引入均會強化反應過程。但是這些措施均導致流體的湍流強度的增加，因而都是以增加黏性耗散[34]，即流體流動的阻力為代價的。圖10為不同射流角度和射流流速下的黏性耗散云圖。圖10表明，相同射流角度下，反應器內黏性耗散會隨噴射流速變大而增加；相同射流流速下，反應器內的黏性耗散會隨射流角度的增大而增加。黏性耗散增加意味著流體阻力以及相應的流體輸送能量的增加。因此，反應器設計需要權衡反應轉化率和黏性耗散，來確定最優的射流角度和流速參數。圖11為甲烷增加的轉化率與黏性耗散關系，可以看出，固定黏性耗散的情況下，隨著角度的增大，甲烷的轉化率大體呈現增加的趨勢。但與圖8相類似，不同射流角度的曲線會存在交叉，這意味著，相同黏性耗散下并不總是角度越大效果越好。此外，當射流角度為75°、80°和85°時，當黏性耗散增加到一定程度后，甲烷增加的轉化率反而開始下降，這說明存在給定射流角度下，射流流速存在一個最優值，該最優值可使甲烷轉化率的增加值達到最大，而繼續增加射流流速則黏性耗散會繼續增加，甲烷轉化率的增加值反而下降。

圖10 不同射流角度和射流流速下黏性耗散云圖Fig.10 Contour of viscous dissipation under different injection angles and injection flow rates

圖11 甲烷增加的轉化率與黏性耗散關系Fig.11 Relationship between increased methane conversion rate and viscous dissipation

以提高轉化率和降低黏性耗散為目標函數，對所有模擬結果進行離散Pareto最優解的篩選，結果如圖12(a)所示。因模擬工作得到的數據點有限，相應的Pareto最優解并不連續，因此，本文采取插值的方法將離散解連續化。在現有插值方法中，SVR算法插值具有預測誤差小、對原始數據分布要求低、可以在數據量較小時取得較好的擬合效果等優點。本文采用SVR算法對數據點進行插值，利用SVR算法公式對黏性耗散和增加的轉化率進行插值擬合，并以此得到連續的Pareto最優解曲線，如圖12(b)所示。圖12(b)表明，當黏性耗散低于0.7 W時，甲烷增加的轉化率對黏性耗散的增加較為敏感；當黏性耗散高于0.7 W時，這一敏感度顯著降低，即繼續增加轉化率則黏性耗散會急劇增加。因此可以選取該點附近某點作為最優設計點。射流角度為85°，垂直射流速度為3.50 m·s-1時，甲烷增加的轉化率8.8%為最大值。考慮到甲烷提高的轉化率為主要目標函數，取最大甲烷增加轉化率的90%（增加的轉化率為7.92%），作為建議操作點M（0.68 W,7.92%）。該操作點對應的垂直射流流速為2.54 m·s-1，射流角度為85°。

圖12 甲烷增加的轉化率與黏性耗散Pareto最優解Fig.12 Pareto optimal solution of increased methane conversion rate and viscous dissipation

4 結 論

本文引入炭顆粒的生成和聚集模型，考慮反應器內的輻射傳熱，給出了太陽能甲烷高溫裂解過程更加嚴格的模型，并以此模型為基礎對甲烷高溫裂解管式反應器流場結構進行模擬。引入射流與翅片擋板對反應器管內優化流場重現來提高甲烷轉化率。通過模擬，本文對射流角度、射流流速以及翅片擋板高度對甲烷轉化率增加量的影響開展了研究。結果表明，翅片擋板對反應過程強化影響不明顯，而射流角度和流速則對甲烷轉化率有顯著影響。在考察射流角度和流速對甲烷轉化率影響的同時，考察了射流角度和流速對相應的黏性耗散的影響，并以甲烷轉化率增加量和黏性耗散為目標，給出了兩目標Pareto前沿，并據此給出了確定最優射流角度和射流流速的分析方法。在最優操作曲線上，以甲烷增加的轉化率為主要目標函數，同時兼顧較小的黏性耗散，給出了建議操作點，通過本文方法設計的反應器與傳統直管反應器相比甲烷轉化率可提高達7.92%。

應該指出，本文研究中采用的最優流場是基于文獻[20]的最優流場經重現而得，雖然引入炭顆粒產生的動力學以及輻射傳熱方程，但仍采用了較為簡單的單相流模型。因而今后采用更加嚴格的多相流模型使模擬更加接近實際反應過程具有重要研究價值。

符號說明

aλ——吸收系數，m-1

CCH4——甲烷濃度，mol·L-1

DCH4,g——甲烷的擴散系數，m2·s-1

Di——組分i的擴散系數，m2·s-1

Dt——湍流擴散系數，m2·s-1

di——粒子直徑，m

Ea——活化能，kJ·kmol-1

hi——組分焓變，kJ·mol-1

iλ——輻射強度，W·m-2

K——湍流聚集系數，m3·s-1

K0——反應速率常數，s-1

m——單個炭顆粒的質量，kg

N——粒子濃度，m-3

NA——阿伏伽德羅常數，6.02×1023mol?1

[N]——粒子源項，m-3·s-1

p——壓力，Pa

R——熱力學常數，8.314 J·(mol·K)-1

Schem——反應熱，J·(m3·s)-1

Srad——輻射傳熱，J·(m3·s)-1

Sct——湍流Schmidt數

Shg-p——Sherwood數

T——溫度，K

u——速度，m·s-1

Y——質量分數

βCH4,i——達涅克修正系數

εk——能量耗散率

λe——有效熱導率，W·(m·K)-1

μe——有效黏度，kg·(m·s)-1

ν——動力黏度，kg·(m·s)-1

ρ——密度，kg·m-3

σsλ——散射系數，m-1