基于無網格方法的四葉草型催化劑體相溫度分布數值模擬

王 闊，趙 波，秦 波，杜艷澤

(1.中國石化撫順石油化工研究院，遼寧 撫順 113001；2.遼寧石油化工大學理學院)

基于無網格方法的四葉草型催化劑體相溫度分布數值模擬

王 闊1，趙 波2，秦 波1，杜艷澤1

(1.中國石化撫順石油化工研究院，遼寧 撫順 113001；2.遼寧石油化工大學理學院)

以四葉草型及圓柱型加氫裂化催化劑三維體相環境為計算實體，以模擬工業運行溫度的函數作為邊界條件，采用無網格數值方法求解傅里葉傳熱方程，并根據計算結果分析加氫裂化反應過程中工況條件及催化劑尺寸對于催化劑內部溫度分布的影響。結果表明：在加氫裂化反應中，圓柱型及四葉草型催化劑顆粒內部的最高溫度及平均溫度均隨反應放熱量、空速、原料密度、催化劑半徑及催化劑長度的增加而增大，其中反應放熱量和催化劑半徑對催化劑內部溫度的影響更顯著；即使在理想的宏觀等溫反應條件下，兩種催化劑內部仍然是非等溫區域；在工況條件及催化劑半徑和長度均相同時，四葉草型催化劑內部的最高溫度、平均溫度均低于圓柱型催化劑，具有比圓柱型催化劑更均勻的溫度分布。

加氫裂化催化劑 傳熱 無網格法 數值模擬 溫度分布

油品的加氫裂化反應整體上為強放熱反應，溫度和熱量是該類反應最重要的影響因素。有關加氫裂化反應體系的能量恒算以及節能研究的文獻眾多[1-2]，但一般僅限于較大空間尺度或單純反應熱層面[1-5]。對于小型實驗裝置，整個反應系統基本可以近似視為等溫反應系統；而工業裝置的反應系統則可以視為一個近似的絕熱系統。在反應過程中，反應器內部的溫度分布會影響催化劑床層不同位置的反應速率，進而影響該位置的熱量釋放，熱量的傳遞又會引起體系溫度分布的重新調整。這一溫度的正反饋過程在加氫裂化工業生產中十分重要。

加氫裂化反應過程的熱量集中于催化劑顆粒，每個催化劑顆粒都相當于反應體系的一個“微熱源”。這些微熱源的溫度分布情況直接或間接影響反應體系的質量傳遞、熱量傳遞、反應進程以及產品分布。同時，催化劑體相內部的溫度分布對催化劑壽命也有極大的影響。因此，對催化劑體相內部的溫度分布進行研究十分重要。傳統的煉化工業一般采用熱電偶測定溫度，但只能得到反應器內部宏觀點的溫度近似值，無法對反應過程的實時微觀熱效應進行分析。近年來，雖然 “非接觸式”的紅外測溫技術得到廣泛發展[6]，但由于催化劑內部微小而復雜的形態結構，目前沒有任何溫度測量設備能夠對其內部溫度進行測量；同時，由于催化劑外觀幾何形態以及反應條件的復雜性，導致常規導熱模型體系的經驗公式難于適用。本課題基于傅里葉傳熱定律的數值模擬方法來解決這一類問題，以四葉草型和圓柱型加氫裂化催化劑三維體相環境為計算實體，以模擬工業運行溫度的函數作為邊界條件，采用無網格數值方法求解傅里葉傳熱方程，并根據計算結果分析加氫裂化反應過程中工況條件及催化劑尺寸對催化劑內部溫度分布的影響。

1 模型建立與求解過程

1.1 傅里葉傳熱方程

對一般傳熱過程，當存在穩定熱源時，體系溫度對于空間的分布形態可由一個穩態傅里葉傳熱方程確立。對于一個三維體系，一般的數學形式如式(1)所示。

(1)

式中：F(x，y，z)表示體系內溫度對于空間的函數，當空間點處于體相內部(v)時，整個體系由一個偏微分方程描述，其中qv(x，y，z)是體系的熱源函數，該函數表示當體系存在熱源或熱匯時該熱源或熱匯對空間的分布形式；λ表示計算體系的導熱系數，由該體系組成材料的物化性質決定。當溫度點處于體相表面(σ)時，溫度可以用第一型迪里克萊邊界條件描述，可將表面溫度分布表達為函數G(x，y，z)。對極少數具有比較簡單的熱源或熱匯函數形式及比較規則的幾何外觀形態和邊界條件的導熱體系而言，函數F(x，y，z)存在解析形式。由于拉普拉斯型方程的初值問題的解是不穩定的[7]，對于絕大多數復雜實際工程的傳熱體系，對模型方程體系的適定性問題討論極為復雜，式(1)只能由數值方法求解。

1.2 加氫反應過程熱量衡算及相關方程參數的確立

整個加氫裂化過程的精確反應熱計算比較困難，根據已有經驗及文獻，加氫裂化反應的放熱量一般為200～500 kJkg，原料密度一般為0.85～1.05 kgL，反應過程的體積空速一般為0.9～1.5 h-1。通過以上數據可對反應體系進行熱量衡算及方程參數估計。

取極端原料性質及操作條件進行方程參數估計。假定加氫裂化反應放熱量Qm=500 kJkg，油品密度ρ=1.05 kg/L，則每加工1 L油品的放熱量q=Qm×ρ=500 kJ/kg×1.05 kg/L=525 kJ/L。假定體積空速LHSV=1.5 h-1，則每升催化劑每小時的放熱量Qv=q×LHSV=525×1.5=787.5 kJ/(L·h)，也即單位體積熱功率P=787.5 kJ/(L·h)=218.75 W/L=2.187 5×10-4W/mm3。根據文獻[8]，一般催化劑的導熱系數為0.25～0.32 W/(m·K)，因此計算中導熱系數λ取2.58×10-4W/(mm·K)，式(1)中相關熱功率與導熱系數的比值qv(x，y，z)/λ=(2.187 5×10-4W/mm3)/[2.58×10-4W/(mm·K)]=0.85 K/mm2。

1.3 無網格方法及計算過程

目前主要采用數值方法對較為復雜的工程偏微分方程進行求解，常規的數值方法包括有限差分方法(FDM)、有限元方法(FEM)。前者對于復雜邊界條件的適應性很差，在現代工程計算應用中受到了一定的限制；后者自20世紀50年代開始興起，堪稱數值方法領域的一個里程碑，然而它也存在缺點和局限性，例如該方法需要在計算前對計算體系進行相關的網格劃分，而形成網格的計算成本較高，特別是對于復雜的三維體系，生成復雜的高質量網格相當困難[9]。產生這類問題的根源在于有限元計算過程中使用了生成網格單元的相關連接信息。因此，近年來無網格類計算方法應運而生。無網格方法的特點以及與其它計算方法結合使用的案例參見文獻[10-13]。

本課題使用無網格計算方法對三維傅里葉傳熱方程進行數值求解，并對具有不同幾何形態的圓柱型及四葉草型催化劑的體相溫度空間分布形式進行描述和分析，主要計算過程為：對指定采樣點確定相關的三維長方體計算支持域，在支持域中選擇具有局域性質的局域函數，該局域函數以多項式函數為基本函數；首先計算支持域內的計算節點局域矩陣，并選擇三維長方體計算支持域中的三次樣條權重函數；然后在支持域內構造權重矩陣，依據已定義的局域矩陣以及權重矩陣產生對于無網格計算采樣點的形函數矩陣；將采樣點形函數以及形函數的導數依次代入式(1)，構造出以計算點溫度數值為未知向量的線性方程組，該線性方程組的系數矩陣及常數項矢量分別稱為力矩陣和載荷矩陣，在載荷矩陣中可根據具體控制方程及邊界條件加載相應的數據，根據計算需要，力矩陣的行數與列數(即采樣點個數與計算點個數)可相同也可不同；最后用線性最小二乘法求解力矩陣對于載荷矩陣的未知向量，即得到計算節點溫度函數的近似數值，再用計算節點溫度函數的數值與生成的形函數矩陣重構整個采樣節點溫度函數的數值，便可獲得原偏微分方程體系的近似數值解。

1.4 計算參數及邊界條件的設定

相比于傳統圓柱型催化劑，四葉草型催化劑具有較大的宏觀外表面積以及在反應過程中具有較小的系統壓降，因此得到了廣泛應用。四葉草型催化劑的截面外觀如圖1所示。由圖1可知，四葉草型催化劑包含兩類半徑，即4個相切的小圓柱體的半徑r、與4個小圓柱體相切的外接大圓的半徑R。相比于傳統的半徑為R的圓柱型催化劑，四葉草型催化劑具有更小的體積和更大的外表面積。本課題對具有不同大圓半徑R及長度的圓柱型催化劑和四葉草型催化劑內部溫度場進行計算和分析。

圖1 四葉草型催化劑的截面外觀

圓柱型催化劑采樣點的軸向及徑向分布如圖2所示。在整個計算過程中定義的權函數支持域尺寸為2。對于半徑為2 mm、長度為10 mm的圓柱型催化劑的計算，分別在圓柱體的徑向取6個采樣點、平面角方向取20個采樣點、軸向取41個采樣點。

圖2 圓柱型催化劑的采樣點及計算點分布●—采樣點； ●—計算點。 圖3同

四葉草型催化劑的截面外觀相對復雜，對稱性遠較圓柱型催化劑差，在徑向平面上具有更為豐富的信息。四葉草型催化劑采樣點軸向及徑向分布如圖3所示。在整個計算過程中定義的權函數支持域尺寸為2。對于大圓半徑為2 mm、長度為10 mm的四葉草型催化劑的計算，分別在徑向取6個采樣點、徑向平面角方向取40個采樣點、軸向取41個采樣點。

圖3 四葉草型催化劑的采樣點及計算點分布

整個計算過程以外界溫度370 ℃作為邊界條件，通過計算獲得熱點均勻分布[14]的不同形狀催化劑的溫度空間分布情況。

1.5 計算體系

采用分布式計算系統對模型進行計算，整套設備包含常規處理器計算機5臺，共160個核心。計算體系由24口交換機完成計算數據的傳輸和交互[15-16]。整個無網格計算過程在本質上具有較好的并發性，采樣點對應的形函數及其偏導數的計算以及力矩陣及載荷矩陣的計算組裝過程均采用多機并行方式實現，由體系的160個核心并行參與。

2 結果與討論

采用上述模型對不同反應工況及催化劑尺寸條件下熱點分布均勻的圓柱型及四葉草型催化劑在外表面溫度為370 ℃時的內部溫度分布進行計算，所得催化劑體系內部的最高溫度、平均溫度以及溫度分布的標準差見表1。根據圓柱型催化劑在工況7下的三維溫度場分布繪制二維軸向及徑向截面圖，如圖4所示；根據四葉草型催化劑在工況15下的三維溫度場分布繪制二維軸向及徑向截面圖，如圖5所示。

2.1 圓柱型催化劑內部溫度場分析

由表1可知：裂化反應的放熱量、原料密度、反應空速及催化劑半徑和長度均影響圓柱型催化劑的內部溫度分布；催化劑顆粒內部最高溫度及平均溫度均隨反應放熱量、原料密度、反應空速以及催化劑半徑、催化劑軸向長度的增加而增大；反應放熱量及催化劑顆粒半徑對催化劑最高溫度及平均溫度的影響顯著大于催化劑長度、反應空速及原料密度。

對于熱點分布均勻的等溫反應條件，雖然不同計算工況對應的最高溫度及平均溫度數值有所不同，但溫度分布形式基本類似。由圖4可知，圓柱型催化劑內部溫度分布的等勢面在徑向上是嚴格的圓型分布，在軸向剖面上近似橢圓分布，催化劑外表面溫度為370 ℃，由外到內溫度逐漸增加，其徑向溫度梯度大于軸向溫度梯度，在圓柱型催化劑的幾何中心處溫度最高。因此，可以認為雖然在理論上圓柱型催化劑內部并非等溫環境，但其內部的反應環境與等溫環境類似。而當催化劑半徑及長度較大時，催化劑圓柱體中心的最高溫度與表面溫度相差較大，其反應環境相對于等溫反應偏離較大。因此，僅從熱效應上考慮，制備圓柱型催化劑時，其半徑及長度不宜過大。

表1 圓柱型及四葉草型催化劑的內部溫度分布

圖4 圓柱型催化劑在工況7下的內部溫度分布截面圖

圖5 四葉草型催化劑在工況15下的內部溫度分布截面圖

2.2 四葉草型催化劑內部溫度場分析

由表1可知：與圓柱型催化劑相似，裂化反應的放熱量、原料密度、反應空速以及催化劑半徑均影響四葉草型催化劑內部的溫度分布；催化劑顆粒內部的最高溫度及平均溫度均隨反應放熱量、原料密度、反應空速、催化劑外接大圓半徑以及催化劑長度的增加而增大；反應放熱量及催化劑顆粒外接大圓半徑對催化劑最高溫度及平均溫度的影響顯著大于催化劑長度、反應空速及原料密度。

通過對比圓柱型和四葉草型催化劑的內部溫度分布數據可知，在所有工況條件及催化劑半徑和長度均相同時，四葉草型催化劑內部的最高溫度及平均溫度均低于圓柱型催化劑，兩者溫度一般相差0.1～0.9 ℃，且四葉草型催化劑內部溫度的波動范圍小于相應的圓柱型催化劑。

不同工況下四葉草型催化劑內部的溫度分布比較相似。在熱點分布均勻的等溫反應條件下，四葉草型催化劑的截面較圓柱型催化劑復雜，其內部的溫度分布形式也與圓柱型催化劑相差較大，尤其體現出了較為復雜的截面溫度分布信息。由圖5可知，四葉草型催化劑溫度分布的等勢面在截面每個葉上的徑向方向近似呈貝殼型分布，沿催化劑對應的大圓圓心與小圓圓心連線方向的溫度梯度變化最小，沿大圓圓心與小圓公切點方向的溫度梯度變化最大，由于催化劑體系有對稱的特征，最高溫度區域與圓柱型催化劑相似，仍在催化劑的中心軸線的中點位置。在軸向剖面上，其溫度等勢面呈近似圓角矩形分布，催化劑的外表面溫度為370 ℃，由外到內溫度逐漸增加，其徑向的溫度梯度大于軸向的溫度梯度。因此，同樣認為雖然四葉草型催化劑的內部并非等溫環境，但其內部的反應環境與等溫環境類似。與圓柱型催化劑相似，當四葉草型催化劑外接大圓半徑及長度較大時，催化劑中心的最高溫度與表面溫度相差較大，其反應環境相對于等溫反應偏離較大。另外，由于計算機硬件系統內存的限制，沿著大圓圓心與小圓公切點方向及其鄰域的計算點與采樣點數目有限，導致對該方向及其鄰域內溫度變化的描述具有一定的誤差。

3 結 論

在加氫裂化反應中，圓柱型及四葉草型催化劑顆粒內部的最高溫度及平均溫度均隨反應放熱量、空速、原料密度、催化劑半徑及催化劑長度的增加而增大，其中反應放熱量和催化劑半徑對催化劑內部溫度的影響更顯著。即使在理想的宏觀等溫反應條件下，兩種催化劑內部仍然是非等溫區域。在工況條件及催化劑半徑和長度均相同時，四葉草型催化劑內部的最高溫度、平均溫度均低于圓柱型催化劑，具有比圓柱型催化劑更均勻的溫度分布。

[1] 尹兆林.煉油企業全廠用能分析及節能優化[J].石油煉制與化工，2012，43(10)：86-91

[2] 董兆海，袁永新，王明傳.加氫裂化裝置能耗及節能分析[J].齊魯石油化工，2011，39(2)：87-91

[3] 劉小波，毛羽，王娟，等.基于多孔介質加氫裂化反應器多相流數值模擬[J].石油學報(石油加工)，2012，28(2)：260-267

[4] 龍軍，邵潛，賀振富，等.規整結構催化劑及反應器研究進展[J].化工進展，2004，23(9)：925-931

[5] 梁文杰，闕國和.石油化學[M].2版.東營：中國石油大學出版社，2011：356

[6] 鄭忠，何臘梅.紅外測溫技術及在鋼鐵生產中的應用[J].工業加熱，2005，34(3)：25-29

[7] 廖玉麟.數學物理方程[M].上海：華東理工大學出版社，1995：12-13

[8] 吳建民，張海濤，應衛勇，等.鈷基催化劑固定床有效導熱系數[J].過程工程學報，2010，10(1)：29-34

[9] 劉桂榮，顧元通.無網格法理論及程序設計[M].濟南：山東大學出版社，2007：24

[10]劉欣.無網格方法[M].北京：科學出版社，2011：9

[11]貝新源，岳宗五.三維SPH程序及其在斜高速碰撞問題中的應用[J].計算物理，1997，14(2)：155-166

[12]張雄，劉巖.無網格方法[M].北京，清華大學出版社，2001：42

[13]鐘萬勰.彈性力學求解新體系[M].大連：大連理工大學出版社，1995：35

[14]莫爾比代利M，加夫里迪斯A，瓦爾馬A.催化劑設計[M].王安杰，李翔，趙蓓，等，譯.北京：化學工業出版社，2004：128

[15]王曉丹，吳崇明.基于MATLAB的系統分析與設計[M].西安：西安電子科技大學出版社，2000：125

[16]劉維.實戰Matlab之并行程序設計[M].北京：北京航空航天大學出版社，2012：135-137

NUMERICAL SIMULATION OF TEMPERATURE DISTRIBUTION IN FOUR-LEAF SHAPE CATALYST BY MESHLESS METHOD

Wang Kuo1， Zhao Bo2， Qin Bo1， Du Yanze1

(1.SINOPECFushunResearchInstituteofpetroleumandPetrochemical，Fushun，Liaoning113001； 2.LiaoningShihuaUniversity，CollegeofScience)

Based on three-dimensional bulk of real four-leaf type hydrocracking catalyst， the meshfree calculation to solve the fourier partial differential equation for heat transmission was used to simulate the influence of external temperature fluctuation on the internal temperature distribution in the catalyst using industrial operating temperature and catalysts size as the boundary conditions. The analysis results show that the reactions in the catalyst do not occur under real isothermal surroundings. The maximum and average temperature increases in the cylindrical and four-leaf type catalyst as the increase of reaction heat release， space velocity， feed density， catalyst radius and length of the catalyst. Among them， the heat release and catalyst radius have much more influence. Even under ideal apparent isothermal reaction conditions， the non-isothermal area still exists in the two catalysts. When the operation conditions and catalyst size are the same， the four-leaf type catalyst has the highest internal temperature， while its average temperature and distribution are lower and better than cylindrical one.

hydrocracking catalyst； heat transmission； meshless method； numerical simulation； temperature distribution

2016-02-29； 修改稿收到日期： 2016-05-18。

王闊，碩士，工程師，研究方向為工業催化及分子模擬。

王闊，E-mail：wangkuo.fshy@sinopec.com。

中國石油化工股份有限公司項目(JQ-011308)。