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

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

基于無網(wǎng)格方法的四葉草型催化劑體相溫度分布數(shù)值模擬

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

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

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

加氫裂化催化劑 傳熱 無網(wǎng)格法 數(shù)值模擬 溫度分布

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

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

1 模型建立與求解過程

1.1 傅里葉傳熱方程

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

(1)

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

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

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

取極端原料性質(zhì)及操作條件進行方程參數(shù)估計。假定加氫裂化反應放熱量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。根據(jù)文獻[8]，一般催化劑的導熱系數(shù)為0.25～0.32 W/(m·K)，因此計算中導熱系數(shù)λ取2.58×10-4W/(mm·K)，式(1)中相關熱功率與導熱系數(shù)的比值qv(x，y，z)/λ=(2.187 5×10-4W/mm3)/[2.58×10-4W/(mm·K)]=0.85 K/mm2。

1.3 無網(wǎng)格方法及計算過程

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

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

1.4 計算參數(shù)及邊界條件的設定

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

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

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

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

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

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

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

1.5 計算體系

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

2 結果與討論

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

2.1 圓柱型催化劑內(nèi)部溫度場分析

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

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

表1 圓柱型及四葉草型催化劑的內(nèi)部溫度分布

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

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

2.2 四葉草型催化劑內(nèi)部溫度場分析

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

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

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

3 結 論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[16]劉維.實戰(zhàn)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，F(xiàn)ushun，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。

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

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

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