王 闊

(中國石化 撫順石油化工研究院， 遼寧 撫順 113001)

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

王 闊

(中國石化 撫順石油化工研究院， 遼寧 撫順 113001)

以真實(shí)圓柱型加氫裂化催化劑三維體相環(huán)境為計(jì)算實(shí)體，以模擬工業(yè)運(yùn)行溫度的定值函數(shù)作為邊界條件，采用無網(wǎng)格數(shù)值方法求解傅里葉傳熱方程，并使用計(jì)算結(jié)果分析加氫裂化催化反應(yīng)過程中外界溫度對于催化劑內(nèi)部溫度分布的影響。結(jié)果表明，實(shí)際加氫裂化反應(yīng)過程中，催化劑顆粒內(nèi)部并非等溫，外界的反應(yīng)溫度以及催化劑顆粒內(nèi)部的熱點(diǎn)分布對于催化劑團(tuán)簇內(nèi)部的溫度場分布有一定的影響。催化劑體相內(nèi)部的平均溫度也隨著反應(yīng)體系放熱情況、催化劑粒徑、原料油密度、反應(yīng)空速以及催化劑內(nèi)部熱點(diǎn)分布情況的不同而有所變化。

加氫裂化； 傳熱； 無網(wǎng)格； 數(shù)值模擬； 分布計(jì)算

一般來講，加氫裂化反應(yīng)是一種整體具有強(qiáng)烈放熱效應(yīng)的反應(yīng)，溫度和熱量是該類反應(yīng)最重要的影響因素。對于加氫裂化反應(yīng)體系的能量恒算以及節(jié)能研究眾多[1-2]，但一般僅限于較大空間尺度或單純反應(yīng)熱層面[1-3]。在科學(xué)實(shí)驗(yàn)和工業(yè)生產(chǎn)實(shí)踐中，對于小型實(shí)驗(yàn)裝置而言，整個(gè)反應(yīng)系統(tǒng)基本可以近似視為等溫反應(yīng)系統(tǒng)；對工業(yè)裝置而言，反應(yīng)系統(tǒng)可以被視為一個(gè)比較嚴(yán)格的絕熱系統(tǒng)。在反應(yīng)過程中，反應(yīng)器內(nèi)部的溫度分布必然會影響到催化劑床層不同位置的反應(yīng)速率，進(jìn)而影響到相關(guān)位置的熱量釋放。對于反應(yīng)系統(tǒng)尤其是與工業(yè)裝置緊密相關(guān)的絕熱系統(tǒng)而言，反應(yīng)產(chǎn)生熱量的傳遞必然引起體系溫度分布的重新調(diào)整。這一溫度的正反饋過程在加氫裂化工業(yè)生產(chǎn)中十分重要。

審視整個(gè)加氫裂化溫度正反饋過程不難發(fā)現(xiàn)，在絕熱反應(yīng)體系中，整個(gè)過程的熱量正是來自于催化劑顆粒自身，每個(gè)催化劑顆粒都相當(dāng)于反應(yīng)體系的1個(gè)“微熱源”。這些微熱源的溫度分布情況直接或間接的影響到反應(yīng)體系的質(zhì)量傳遞、熱量傳遞和反應(yīng)進(jìn)程以及產(chǎn)品分布，催化劑體相內(nèi)部的溫度分布對于催化劑的壽命也有著極大的影響。因此，對于催化劑內(nèi)部體相溫度分布的描述意義十分重大。

對于催化劑內(nèi)部體相溫度分布的研究十分困難。傳統(tǒng)的煉化工業(yè)一般采用熱電偶測定溫度，但只能測定得到反應(yīng)器內(nèi)部的宏觀點(diǎn)的溫度近似值，無法分析反應(yīng)過程的實(shí)時(shí)微觀熱效應(yīng)。雖然近年來“非接觸式”的紅外測溫技術(shù)得到廣泛發(fā)展[4]，但由于催化劑內(nèi)部微小而復(fù)雜的形態(tài)結(jié)構(gòu)，還沒有任何溫度測量設(shè)備能夠在如此微觀尺度進(jìn)行溫度測量；同時(shí)，由于催化劑外觀幾何形態(tài)以及反應(yīng)溫度條件的復(fù)雜性，導(dǎo)致常規(guī)導(dǎo)熱模型體系的經(jīng)驗(yàn)公式難于適用。因此，基于傅里葉傳熱定律的數(shù)值模擬方法成為解決這一類問題的有效途徑。

1 圓柱型催化劑體相溫度分布數(shù)值模擬理論及計(jì)算部分

1.1 傅里葉傳熱方程簡介

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

(1)

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

1.2 加氫反應(yīng)過程熱量衡算及相關(guān)方程參數(shù)確立

在一般的加氫裂化反應(yīng)過程中，催化劑體系內(nèi)部的溫度分布可由三維的傅里葉傳熱方程描述。但實(shí)際工業(yè)反應(yīng)條件下的熱源函數(shù)以及催化劑表面溫度分布十分復(fù)雜，因此，對三維的傅里葉傳熱方程采用無網(wǎng)格數(shù)值方法進(jìn)行求解。

一般而言，整個(gè)加氫裂化過程的精確反應(yīng)熱計(jì)算比較困難，根據(jù)已有經(jīng)驗(yàn)及相關(guān)文獻(xiàn)，加氫裂化反應(yīng)放熱范圍一般為200～500 J/g之間，而原料油密度一般在0.85～1.05 kg/L之間，反應(yīng)過程的體積空速一般在0.9～1.5 h-1之間。通過這些數(shù)據(jù)，可對反應(yīng)體系進(jìn)行熱量衡算及相關(guān)方程參數(shù)估計(jì)。

例如，假定反應(yīng)放熱Qm=500 J/g，油品密度ρ=1.05 kg/L，每加工1 L油品放熱量q=Qm·ρ=500 J/g×1.05 kg/L=525 kJ/L。假定相關(guān)反應(yīng)空速LHSV=1.5 h-1，則每升催化劑每小時(shí)放熱Qv=q·LHSV=525 kJ/L×1.5 h-1=787.5 kJ/(L·h)，也即其單位體積熱功W=787.5 kJ/(L·h)=218.75 W/L=2.1875×10-4W/mm3。而一般催化劑的導(dǎo)熱系數(shù)λ在0.25～0.32 W/(m·K)范圍[10]。在此取導(dǎo)熱系數(shù)為2.58×10-4W/(mm·K)，因此，式(1)中相關(guān)熱功率與導(dǎo)熱系數(shù)的比值qv(x,y,z)/λ=(2.1875×10-4W/mm3)/(2.58×10-4W/(mm·K))= 0.85 K/mm2。

1.3 無網(wǎng)格方法簡介及計(jì)算過程介紹

較為復(fù)雜的工程偏微分方程目前主要采用數(shù)值方法進(jìn)行求解，常規(guī)的數(shù)值方法包括有限差分方法(FDM)和有限元方法(FEM)。

FDM應(yīng)用最早，但其對于復(fù)雜邊界條件的適應(yīng)性很差。因此，在現(xiàn)代工程計(jì)算應(yīng)用中受到了一定的限制。而FEM自20世紀(jì)50年代開始興起，堪稱數(shù)值方法領(lǐng)域的一個(gè)里程碑。該方法的主導(dǎo)思想是將1個(gè)復(fù)雜形態(tài)的連續(xù)體劃分成有限個(gè)單元，各個(gè)單元通過一定的規(guī)則連接起來。正是基于這一指導(dǎo)思想，使得FEM在處理復(fù)雜幾何形體的各種線性和非線性問題中表現(xiàn)出很強(qiáng)的靈活性和魯棒性，在工程領(lǐng)域得到廣泛應(yīng)用，同時(shí)出現(xiàn)了大量的商業(yè)軟件包。

然而FEM也有其固有的缺點(diǎn)和局限性，例如有限元計(jì)算需要在計(jì)算前對計(jì)算體系進(jìn)行相關(guān)的網(wǎng)格劃分。而形成相關(guān)計(jì)算網(wǎng)格計(jì)算成本比較高，特別是對于復(fù)雜的三維體系而言，生成復(fù)雜的高質(zhì)量三維網(wǎng)格相當(dāng)困難[11]。產(chǎn)生這類問題的根源在于在有限元計(jì)算過程中使用了生成網(wǎng)格單元的相關(guān)連接信息。形成一種不使用單元與網(wǎng)格劃分的方法十分重要，因此，近年來無網(wǎng)格類計(jì)算方法應(yīng)運(yùn)而生[14-16]。

在本研究中，嘗試使用無網(wǎng)格計(jì)算方法對三維傅里葉傳熱方程進(jìn)行數(shù)值求解，進(jìn)而實(shí)現(xiàn)對于具有不同熱點(diǎn)分布形態(tài)的圓柱型催化劑體相溫度空間分布形式的描述和分析。以下介紹整個(gè)計(jì)算的大略求解過程。

首先，對于指定采樣點(diǎn)確定相關(guān)的三維長方體計(jì)算支持域。在支持域中選擇具有局域性質(zhì)的局域函數(shù)，該局域函數(shù)用多項(xiàng)式函數(shù)為基，相關(guān)基底函數(shù)L可以表示為三維二次基底函數(shù)集合如式(2)所示。其基底函數(shù)對應(yīng)的矩陣規(guī)模為1×n。在計(jì)算過程中，首先計(jì)算在支持域內(nèi)的計(jì)算點(diǎn)局域矩陣PT，如式(3)所示。

L={1,x,xy,xz,y,yz,z,x2,y2,z2}

(2)

(3)

式(3)中，PT的規(guī)模為m×n。m代表相關(guān)采樣點(diǎn)計(jì)算支持域內(nèi)所包含的計(jì)算點(diǎn)數(shù)目，n則代表選定的三維二次基底函數(shù)的項(xiàng)數(shù)，在本研究中為10。

其次，選擇在相關(guān)的三維長方體計(jì)算支持域中的3次樣條權(quán)重函數(shù)ωx(Rx)，如式(4)所示。式(4)中的Rx由式(5)計(jì)算。

(4)

(5)

在式(5)中，xi為在支持域中的每個(gè)計(jì)算點(diǎn)的x軸坐標(biāo)，而x為采樣點(diǎn)的x軸坐標(biāo)；Rω，x為定義的權(quán)函數(shù)支持域在x軸方向的尺寸。在支持域外部的計(jì)算點(diǎn)權(quán)重值定義為0。

按此規(guī)則定義函數(shù)ωy(Ry)以及ωy(Ry)，最終定義權(quán)函數(shù)如式(6)所示。

ω(x-xi,y-yi,z-zi)=ωx(Rx)ωy(Ry)ωz(Rz)

(6)

按照式(6)在支持域內(nèi)構(gòu)造權(quán)重矩陣Ω， 如式(7)所示。

(7)

式(7)中，Ω的規(guī)模為m×m。m代表相關(guān)采樣點(diǎn)計(jì)算支持域所包含的計(jì)算點(diǎn)數(shù)目。

依據(jù)已定義的局域矩陣PT以及權(quán)函數(shù)矩陣Ω分別形成兩類矩陣A及B，如式(8)、(9)所示。

A=PΩPT

(8)

B=PΩ

(9)

式(8)、(9)中，A的規(guī)模為n×n；B的規(guī)模為n×m。

接下來，定義后續(xù)無網(wǎng)格計(jì)算使用的任意采樣點(diǎn)的形函數(shù)矩陣如式(10)所示。

Φ=LA-1B

(10)

依據(jù)相關(guān)定義，形函數(shù)矩陣規(guī)模為1×m，形函數(shù)本身是采樣點(diǎn)位置坐標(biāo)的函數(shù)。

最后確立式(1)中函數(shù)F(x,y,z)近似表達(dá)式的具體形式為式(11)。

(11)

由于函數(shù)F(x,y,z)近似表達(dá)式中形函數(shù)為包含位置采樣點(diǎn)變量x，y以及z的函數(shù)，可將生成的形函數(shù)分別對x，y以及z進(jìn)行偏導(dǎo)求取，其結(jié)果如式(12)～(14)所示。

(12)

(13)

(14)

將原函數(shù)式(11)及求導(dǎo)結(jié)果式(12)～(14)依次以相關(guān)矩陣的形式代入式(1)中的控制方程及邊界條件方程，經(jīng)過整理獲得以計(jì)算點(diǎn)溫度數(shù)值為未知向量的線性方程組。對應(yīng)線性方程組的系數(shù)矩陣及常數(shù)項(xiàng)矢量分別稱為 “力矩陣”和“載荷矩陣”。在“載荷矩陣”中可以根據(jù)具體的控制方程及邊界條件加載相關(guān)數(shù)據(jù)；根據(jù)計(jì)算需要，相關(guān)的“力矩陣”行數(shù)與列數(shù)，即采樣點(diǎn)個(gè)數(shù)與計(jì)算點(diǎn)個(gè)數(shù)可以相同也可以不同，采樣點(diǎn)位置與計(jì)算點(diǎn)位置可以重合也可以不重合，這也是無網(wǎng)格方法的靈活性之一。

最后，用線性最小二乘法求解力矩陣對于載荷矩陣的未知向量，即得到計(jì)算點(diǎn)溫度函數(shù)的近似數(shù)值。再用計(jì)算點(diǎn)溫度函數(shù)的數(shù)值與生成的形函數(shù)矩陣重構(gòu)整個(gè)采樣點(diǎn)溫度函數(shù)數(shù)值，即可獲得原偏微分方程體系的近似數(shù)值解。

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

由于圓柱型催化劑幾何實(shí)體形狀的對稱性比較好，在計(jì)算過程中，將直角坐標(biāo)系映射為柱坐標(biāo)系。分別在柱體的半徑方向取6個(gè)徑向采樣點(diǎn)、徑向截面角方向取20個(gè)采樣點(diǎn)、在軸向方向取41個(gè)采樣點(diǎn)，同時(shí)生成5個(gè)徑向計(jì)算點(diǎn)、由于圓柱型催化劑幾何實(shí)體形狀的對稱性比較好，在計(jì)算過程中，將直角坐標(biāo)系映射為柱坐標(biāo)系。分別在柱體的半徑方向取6個(gè)徑向采樣點(diǎn)、徑向截面角方向取20個(gè)采樣點(diǎn)、在軸向方向取41個(gè)采樣點(diǎn)，同時(shí)生成5個(gè)徑向計(jì)算點(diǎn)、20個(gè)徑向截面角方向計(jì)算點(diǎn)和40個(gè)軸向計(jì)算點(diǎn)，其空間分布如圖1所示。在整個(gè)計(jì)算過程中，定義的權(quán)函數(shù)支持域尺寸為2。

圖1 圓柱體采樣點(diǎn)及計(jì)算點(diǎn)分布Fig.1 Distribution of sampling and calculation points about cylinder

在催化劑工作過程中，催化劑內(nèi)部溫度較其鄰域更高的區(qū)域稱為“熱點(diǎn)”。實(shí)際催化劑內(nèi)部熱點(diǎn)的成因及分布十分復(fù)雜！由于實(shí)際催化劑制備過程及催化反應(yīng)的極端復(fù)雜性，一般不能保證催化劑完全浸漬均勻；有時(shí)因特定的反應(yīng)，甚至需要專門制備活性組分為蛋黃型或蛋殼型分布的非均勻催化劑[12]。也就是說，催化劑圓柱體內(nèi)部很可能出現(xiàn)活性組分分布不均勻的現(xiàn)象，該現(xiàn)象導(dǎo)致了催化劑內(nèi)部熱點(diǎn)分布不均勻，或出現(xiàn)突出熱點(diǎn)的情況。因此，相關(guān)的數(shù)值模擬應(yīng)盡量考慮到實(shí)際工業(yè)催化劑的對應(yīng)情況。為簡化問題，在整個(gè)計(jì)算過程中模擬了熱點(diǎn)均勻分布、熱點(diǎn)蛋黃型分布以及熱點(diǎn)蛋殼型分布3類不同類型的催化劑的溫度分布情況。

3類模擬計(jì)算均假定催化劑外部溫度為恒溫370℃。同時(shí)確保整個(gè)催化劑圓柱型顆粒的放熱功率與混合均勻情況下的放熱功率一致，但不同熱點(diǎn)分布形式對應(yīng)的催化劑熱點(diǎn)的功率強(qiáng)度和位置則發(fā)生變化。也即調(diào)整相關(guān)參數(shù)使得式(15)保持不變。

(15)

在整個(gè)計(jì)算過程中，假定圓柱型催化劑長度為20 mm，半徑為2 mm，功率值W為0.1202 W。

第1類情況假定，催化劑內(nèi)部熱點(diǎn)完全均勻，催化劑熱功率空間分布形式對于每種反應(yīng)工況為定值，但隨具體反應(yīng)工況的不同而有所變化。其相關(guān)催化劑體系溫度分布計(jì)算序號為1-16。

第2類情況假定，催化劑內(nèi)只有1個(gè)蛋黃型分布的熱點(diǎn)區(qū)域。且熱點(diǎn)位置不發(fā)生變化，催化劑熱功率空間分布形式由式(16)描述。在計(jì)算過程中，調(diào)整相關(guān)參數(shù)C、D使得式(15)成立。所選3種條件對應(yīng)計(jì)算序號及相關(guān)參數(shù)列于表1，相關(guān)熱功率函數(shù)圖像如圖2所示。

qv(x,y,z)=Ce(-D((x-x1)2+(y-y1)2))

(16)

第3類情況假定，催化劑內(nèi)部熱點(diǎn)分布呈蛋殼型，且最高熱點(diǎn)位置發(fā)生變化，催化劑熱功率空間分布形式由式(17)描述。

qv(x,y,z)=A1e(-B1((x-x1)2+(y-y1)2))-
A1e(-B2((x-x1)2+(y-y1)2))

(17)

圖2 熱點(diǎn)蛋黃型分布的催化劑的熱功率函數(shù)Fig.2 Power function of the catalyst with hotspots eggyolk distribution

同樣在計(jì)算過程中調(diào)整相關(guān)參數(shù)A1、B1及B2使得式(15)成立，所選3種條件對應(yīng)計(jì)算序號及相關(guān)參數(shù)列于表2，相關(guān)熱功率函數(shù)圖像如圖3所示。

表1 熱點(diǎn)蛋黃型分布的催化劑相關(guān)功率函數(shù)的計(jì)算參數(shù)Table 1 Calculation parameters of power function about the catalyst with hot spots eggyolk distribution

表2 熱點(diǎn)蛋殼型分布的催化劑相關(guān)功率函數(shù)的計(jì)算參數(shù)Table 2 Calculation parameters of power function about the catalyst with hot spots eggshell distribution

圖3 熱點(diǎn)蛋殼型分布的催化劑的熱功率函數(shù)Fig.3 Power function of the catalyst with hotspots egg shell distribution

1.5 模型所涉及的計(jì)算體系

整個(gè)計(jì)算過程使用的是分布式的高密度計(jì)算系統(tǒng)。其整套設(shè)備包含常規(guī)處理器計(jì)算機(jī)5臺，共計(jì)16整個(gè)計(jì)算過程使用的是分布式的高密度計(jì)算系統(tǒng)。其整套設(shè)備包含常規(guī)處理器計(jì)算機(jī)5臺，共計(jì)160核心。整個(gè)計(jì)算集群包含顯示體系和計(jì)算體系。顯示體系通過KVM切換器有遠(yuǎn)程主機(jī)控制顯示。計(jì)算體系則有24口交換機(jī)完成計(jì)算數(shù)據(jù)的傳輸和交互。計(jì)算體系的每個(gè)計(jì)算點(diǎn)可以單獨(dú)并行使用，也可以整體分布使用。該系統(tǒng)十分適合高密度計(jì)算占優(yōu)的系統(tǒng)體系。整個(gè)計(jì)算過程采用并行計(jì)算方式進(jìn)行，形函數(shù)及其偏導(dǎo)數(shù)計(jì)算以及力矩陣及載荷矩陣的組裝過程采用多機(jī)并行計(jì)算，共有體系的160核心并行參與。

2 結(jié)果與討論

2.1 相關(guān)計(jì)算結(jié)果

對總計(jì)22種不同的工況條件進(jìn)行計(jì)算，并對其中7種情況的三維溫度分布分別作二維切面圖和三維等勢圖。

2.2 等溫反應(yīng)條件下熱點(diǎn)均勻型催化劑內(nèi)部溫度場分析

計(jì)算了不同反應(yīng)工況及催化劑尺度條件下熱點(diǎn)分布均勻的催化劑在外表面370℃進(jìn)行反應(yīng)時(shí)的催化劑內(nèi)部溫度分布情況，將計(jì)算所得催化劑體系內(nèi)部的最高溫度、平均溫度以及溫度分布的標(biāo)準(zhǔn)差列于表3。將7號計(jì)算對應(yīng)的三維溫度場分布作二維軸向及徑向切面圖，結(jié)果示于圖4。

表3 熱點(diǎn)分布均勻的圓柱型催化劑表面370℃反應(yīng)條件下的體相溫度性質(zhì)Table 3 Bulk phase temperature properties of the cylindrical catalyst with hot spots homogeneousdistribution and uniform temperature of 370℃ on surface

圖4 序號7條件催化劑內(nèi)部溫度分布切片圖Fig.4 The internal temperature distribution of the catalyst in the serial number 7

對表3的分析發(fā)現(xiàn)，相關(guān)裂化反應(yīng)的反應(yīng)熱、原料油密度、反應(yīng)空速以及催化劑半徑均直接或間接地影響催化劑內(nèi)部的溫度分布，且所對應(yīng)的相關(guān)不同工藝條件下催化劑內(nèi)部的溫度分布比較相似。僅就單因素分析來看，催化劑顆粒內(nèi)部最高溫度及平均溫度均隨反應(yīng)熱、原料油密度、反應(yīng)空速以及催化劑半徑和催化劑的軸向長度的增加而增大；就其影響效果而言，反應(yīng)熱及催化劑顆粒半徑對于催化劑最高溫度及平均溫度的影響要顯著的大于催化劑顆粒長度、反應(yīng)空速及原料油密度的影響。

對于熱點(diǎn)分布均勻的等溫反應(yīng)條件，雖然不同計(jì)算工況對應(yīng)的最高溫度及平均溫度數(shù)值有所不同，但其溫度分布形式基本類似?？梢哉J(rèn)為，溫度分布的等勢面在徑向方向是嚴(yán)格的圓型分布，而在軸向剖面上比較近似橢圓分布；在催化劑的外表面溫度為370℃時(shí)，催化劑內(nèi)部溫度沿徑向和軸向逐漸增加，且徑向方向溫度梯度大于軸向方向的，在圓柱型催化劑中心的溫度最高?？梢哉J(rèn)為，雖然在理論上催化劑內(nèi)部并非等溫環(huán)境，但其內(nèi)部的反應(yīng)環(huán)境與等溫環(huán)境比較類似。但當(dāng)催化劑顆粒半徑及長度較大時(shí)，催化劑圓柱體中心最高溫度與表面溫度則相差較大，在該催化劑內(nèi)部的反應(yīng)相對于等溫反應(yīng)偏離較大。因此，僅就熱效應(yīng)而言，制備圓柱型催化劑時(shí)，其顆粒半徑及長度一般不宜過大。

2.3 等溫條件下不同熱點(diǎn)分布催化劑內(nèi)部溫度場分析

模擬計(jì)算獲得的熱點(diǎn)蛋黃型分布催化劑的溫度場信息及分布圖分別列于表4及圖5。

表4 熱點(diǎn)蛋黃型分布催化劑的體相系溫度信息Table 4 Bulk phase temperature properties of the cylindrical catalyst with hot spots eggyolk distribution

The uniform temperature of 370℃ on surface

圖5 熱點(diǎn)蛋黃分布催化劑對應(yīng)的體相溫度軸向及徑向分布圖Fig.5 The axial and radial temperature distribution of the catalyst with hot spots eggyolk distributionSerical number: (a) 17; (b) 18; (c) 19

從表4及圖5可以發(fā)現(xiàn)，第17～19號計(jì)算條件所對應(yīng)蛋黃型催化劑熱點(diǎn)分布由集中逐步過渡到分散，所對應(yīng)的熱功率空間分布形式逐漸呈均勻分布，第19號計(jì)算條件所對應(yīng)熱功率已接近完全均勻分布；催化劑體系對應(yīng)的最高溫度及平均溫度都逐漸降低，最終接近熱點(diǎn)均勻的催化劑體系；催化劑內(nèi)部的軸向及徑向溫度分布的梯度逐漸降低，表征催化劑內(nèi)部溫度波動(dòng)的標(biāo)準(zhǔn)差逐漸減小，催化劑內(nèi)部更接近等溫環(huán)境。

將模擬計(jì)算獲得的熱點(diǎn)蛋殼型分布的催化劑的溫度場信息及分布圖分別列于表5及圖6。

從表5及圖6可以發(fā)現(xiàn)，第20～22號計(jì)算條件所對應(yīng)蛋殼型催化劑熱點(diǎn)分布由催化劑徑向中心逐步過渡到徑向邊緣，同時(shí)相關(guān)的熱功率函數(shù)峰值也有所降低；隨之催化劑體系對應(yīng)的最高溫度及平均溫度都逐漸降低，最終甚至低于熱點(diǎn)均勻的催化劑體系；同時(shí)，催化劑內(nèi)部的軸向及徑向溫度分布的梯度逐漸降低，催化劑內(nèi)部接近等溫的區(qū)域逐漸增加，表征催化劑內(nèi)部溫度波動(dòng)的標(biāo)準(zhǔn)差逐漸減小，催化劑內(nèi)部較蛋黃分布催化劑更接近等溫環(huán)境。

表5 熱點(diǎn)蛋殼型分布催化劑的體相溫度信息Table 5 Bulk phase temperature properties of the cylindrical catalyst with hot spots eggshell distribution

Uniform temperature of 370℃ on surface

圖6 熱點(diǎn)蛋殼分布催化劑對應(yīng)的溫度軸向及徑向分布圖Fig.6 The axial and radial temperature distributions of the catalyst with hot spot eggshell distributionSerical number: (a) 20; (b) 21; (c) 22

3 結(jié) 論

(1)即使在理想的宏觀等溫反應(yīng)條件下，圓柱型催化劑顆粒內(nèi)部仍然是具有非恒定溫度場分布的非等溫反應(yīng)區(qū)域。

(2)在常規(guī)加氫裂化反應(yīng)中，較高的反應(yīng)熱、較高的處理空速、較大的原料油密度以及較大的催化劑顆粒半徑和軸向長度均導(dǎo)致較高的催化劑內(nèi)部溫度。其中，反應(yīng)熱和催化劑顆粒半徑對于催化劑內(nèi)部溫度的影響更顯著。

(3)熱點(diǎn)蛋殼型分布的圓柱型催化劑內(nèi)部溫度梯度及溫度波動(dòng)相對最低，熱點(diǎn)均勻分布的圓柱型型催化劑次之，熱點(diǎn)蛋黃型分布的圓柱型催化劑的熱點(diǎn)更為集中。熱點(diǎn)分布的過度集中更容易導(dǎo)致催化劑的飛溫現(xiàn)象。

符號說明：

A——無網(wǎng)格計(jì)算使用的中間生成矩陣；

A1——蛋殼型催化劑功率計(jì)算的調(diào)制參數(shù)；

B——無網(wǎng)格計(jì)算使用的中間生成矩陣；

B1——蛋殼型催化劑功率計(jì)算的調(diào)制參數(shù)；

B2——蛋殼型催化劑功率計(jì)算的調(diào)制參數(shù)；

C——蛋黃型催化劑功率計(jì)算的調(diào)制參數(shù)；

D——蛋黃型催化劑功率計(jì)算的調(diào)制參數(shù)；

F(x,y,z)——體系內(nèi)溫度對于空間的函數(shù)；

G(x，y，z)——體系的表面溫度函數(shù)；

L——無網(wǎng)格計(jì)算使用的三維二次基底函數(shù)集合；

LHSV ——加氫裂化反應(yīng)空速， h-1；

PT——無網(wǎng)格計(jì)算使用的計(jì)算點(diǎn)局域矩陣；

Qm——加氫裂化反應(yīng)熱， J/g；

Qv——單位體積催化劑單位時(shí)間放熱，kJ/(L· h)；

qv(x,y,z)——體系的熱源函數(shù)；

x，y，z——體系的空間坐標(biāo)；

W——催化劑單位體積熱功率，W/mm3；

λ——催化劑顆粒導(dǎo)熱系數(shù)，W/(m·K)；

ρ——加氫裂化反應(yīng)原料油密度， kg/L；

Φ(x,y,z)——無網(wǎng)格生成的形函數(shù)矩陣；

ωx(Rx)——無網(wǎng)格計(jì)算使用的x方向的3次樣條權(quán)重函數(shù)；

ωy(Ry)——無網(wǎng)格計(jì)算使用的y方向的3次樣條權(quán)重函數(shù)；

ωz(Rz)——無網(wǎng)格計(jì)算使用的z方向的3次樣條權(quán)重函數(shù)；

ω(x-xi,y-yi,z-zi)——無網(wǎng)格計(jì)算使用的3維權(quán)重樣條函數(shù)；

Ω(x,y,z)——無網(wǎng)格計(jì)算使用的三維權(quán)重矩陣。

[1] 尹兆林. 煉油企業(yè)全廠用能分析及節(jié)能優(yōu)化[J].石油煉制與化工， 2012， 43(10): 86-91. (YIN Zhaolin. Energy consumption analysis and optimization of a refining enterprise[J].Petroleum Processing and Petrochemicals, 2012， 43(10): 86-91.)

[2] 董兆海, 袁永新, 王明傳. 加氫裂化裝置能耗及節(jié)能分析[J].齊魯石油化工, 2011, 39(2): 87-91. (DONG Zhaohai， YUAN Yongxin， WANG Mingchuan. Analysis on energy consumption and saving of hydrocracking unit[J].Qilu Petrochemical Technology, 2011, 39(2): 87-91.)

[3] 劉小波， 毛羽， 王娟， 等. 基于多孔介質(zhì)加氫裂化反應(yīng)器多相流數(shù)值模擬[J].石油學(xué)報(bào)(石油加工)， 2012， 20(2)： 206-267. (LIU Xiaobo， MAO Yu， WANG Juan， et al. Computational fluid dynamics of multiphase flow in hydrocracking reactor based on porous media[J].Acta Petrolei Sinica (Petroleum Processing Section), 2012， 20(2)： 206-267.)

[4] 鄭忠， 何臘梅. 紅外測溫技術(shù)及在鋼鐵生產(chǎn)中的應(yīng)用[J].工業(yè)加熱， 2005, 34(3)： 25-29. (ZHENG Zhong， HE Lamei. Infrared temperature measurement technology and its application to steel making process[J].Industrial Heating, 2005, 34(3)： 25-29.)

[5] 王曉丹， 吳崇明. 基于MATLAB的系統(tǒng)分析與設(shè)計(jì)-圖像處理[M].西安： 西安電子科技大學(xué)出版社， 2000.

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

[7] 劉維. 實(shí)戰(zhàn)Matlab之并行程序設(shè)計(jì)[M].北京： 北京航空航天大學(xué)出版社， 2012： 135-137.

[8] 龍軍, 邵潛, 賀振富, 等.規(guī)整結(jié)構(gòu)催化劑及反應(yīng)器研究進(jìn)展[J].化工進(jìn)展, 2004, 23(9): 925-931. (LONG Jun, SHAO Qian, HE Zhenfu, et al. Monolithic catalysts and monolithi rectors[J].Chemical Industry and Engineering Progess, 2004, 23(9): 925-931.)

[9] 梁文杰， 闕國和. 石油化學(xué)[M].東營： 中國石油大學(xué)出版社， 2011： 356.

[10] 吳建民， 張海濤， 應(yīng)衛(wèi)勇， 等. 鈷基催化劑固定床有效導(dǎo)熱系數(shù)[J].過程工程學(xué)報(bào)， 2010, 10(1)： 29-34. (WU Jianmin, ZHANG Haitao, YING Weiyong, et al. Effective thermal conductivity of fixed packing bed of cobalt-based catalyst[J].The Chinese Journal of Process Engineering， 2010, 10(1)： 29-34.)

[11] LIU G R， GU Y T. 無網(wǎng)格法理論及程序設(shè)計(jì)[M].濟(jì)南： 山東大學(xué)出版社， 2007.

[12] 莫爾比代利M， 加夫里迪斯A， 瓦爾馬A. 催化劑設(shè)計(jì)[M].北京： 化學(xué)工業(yè)出版社， 2004.

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

[14] 貝新源， 岳宗五. 三維SPH程序及其在斜高速碰撞問題中的應(yīng)用[J].計(jì)算物理， 1997， 14(2)： 155-166. (BEI Xinyuan, YUE Zongwu. A study on 3 dimension SPH[J]. Chinese Journal of Computation Physics, 1997， 14(2)： 155-166.)

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

[16] 鐘萬勰. 彈性力學(xué)求解新體系[M].大連： 大連理工大學(xué)出版社， 1995.

The Numerical Simulation About Temperature Distribution ofCylinder Catalyst by Meshless Method

WANG Kuo

(FushunResearchInstituteofPetroleumandPetrochemical,SINOPEC,Fushun113001,China)

Based on three-dimensional environment of real cylinder hydrocracking catalyst， the meshless calculation to solve Fourier partial differential equation was provided with industrial operating temperature as the boundary condition. The influence of external temperature on the internal temperature distribution of the catalyst was analyzed with the calculated results. The analysis results showed that during the actual reactions there was not isothermal inside catalyst particle. At the same time, the temperature distribution in bulk catalyst was restricted by the reaction temperature of hydrocracking process outside and internal hotspot distribution of the catalyst particle. The average temperature of bulk catalyst phase was influenced by reaction enthalpy, catalyst particle size, material oil density, reaction space velocity and catalyst inner hotspot distribution.

hydrocracking; thermal conductivity; meshless; numerical simulation; distributed computing

2016-04-20

中國石油化工股份有限公司總和課題(JQ-011308)資助 作者簡介： 王闊，男，工程師，碩士，研究方向?yàn)楣I(yè)催化與分子模擬；E-mail:wangkuo.fshy@sinopec.com

1001-8719(2017)02-0310-10

TQ 116.2

A

10.3969/j.issn.1001-8719.2017.02.016