加氫空冷多相流熱耦合數(shù)值分析

（合肥通用機(jī)械研究院有限公司 國家壓力容器與管道安全工程技術(shù)研究中心，合肥 230031）

0 引言

加氫空冷器在石化生產(chǎn)中的地位舉足輕重，且易發(fā)銨鹽流動(dòng)腐蝕，管束一旦出現(xiàn)這種損傷后果不堪設(shè)想。美國石油學(xué)會(huì)的推薦標(biāo)準(zhǔn)API 932B［1］認(rèn)為空冷介質(zhì)中含有的 NH3、HCl、H2S 等腐蝕性氣體在相應(yīng)壓力、溫度及流速等條件下將形成銨鹽并堵塞管束造成設(shè)備腐蝕泄漏。因此，準(zhǔn)確分析并預(yù)測介質(zhì)在空冷管束中的流速、溫度和分壓等特征參量對確保裝長周期安全運(yùn)行具有重要意義［2-6］。

針對加氫空冷器的銨鹽流動(dòng)腐蝕問題，國內(nèi)外研究機(jī)構(gòu)展開了積極探索：其中，美國多家企業(yè)與研究機(jī)構(gòu)于2000年就銨鹽腐蝕制定了預(yù)測指南，具體成果主要體現(xiàn)在NH4HS對碳鋼等6種材料的腐蝕影響；Sun等［7］則以物理化學(xué)方法，建立了NH4Cl離子平衡反應(yīng)模型，并測試了不同溫度、濕度下碳鋼、不銹鋼等材料的腐蝕特性；偶國富等也通過工藝仿真與流體計(jì)算建立了高壓空冷多相流沖蝕預(yù)測模型［8］。然而綜合分析卻發(fā)現(xiàn)以上成果基本上都將空冷介質(zhì)的流動(dòng)特性與溫度作為獨(dú)立特征參量考察，對流場-溫度場間的關(guān)聯(lián)性缺乏系統(tǒng)闡釋。因此，針對在役工況加氫空冷多相流介質(zhì)的流動(dòng)特征和溫度場展開熱流耦合條件數(shù)值分析，以期探明流場-溫度場間關(guān)聯(lián)及其作用機(jī)制，并為準(zhǔn)確分析及預(yù)測空冷設(shè)備銨鹽腐蝕，預(yù)防設(shè)備失效提供有力支撐。

1 工況及設(shè)備

分析對象為某石化中壓加氫裝置的空冷設(shè)備，設(shè)備所在單元的工藝流程如圖1所示。

圖1 加氫冷卻系統(tǒng)工藝流程Fig.1 Process flow chart of hydrogenation plant cooling system

裝置產(chǎn)能為3.5×103kg/a，其中熱高分空冷器共有16臺(tái)，型號(hào)為A101A-P。每臺(tái)空冷器共2個(gè)進(jìn)口、2個(gè)出口，規(guī)格DN150，換熱管束共6排225根，進(jìn)出口接管和管束的排布如圖2（a）所示?？绽浣橘|(zhì)熱高分氣包含氣相的H2，C1～C4等，少量重組分油C5+和液態(tài)水。空冷器運(yùn)行壓力11.5 MPa、進(jìn)口溫度125 ℃、出口溫度55 ℃，其換熱原理如圖2（b）所示，其下方存在高功率風(fēng)扇進(jìn)行散熱，溫度較高的油氣介質(zhì)在圓柱形管束內(nèi)流動(dòng)，而溫度較低的空氣則自下向上地對其進(jìn)行冷卻。根據(jù)傳熱學(xué)理論，風(fēng)扇以對流方式帶走熱高分氣熱量，而介質(zhì)則通過管束向外界傳熱，形成對流［9］。

圖2 空冷器A101A-P管束分布Fig.2 Distribution of tube bundles of air cooler A101A-P

2 數(shù)理模型及邊界條件

2.1 數(shù)理模型

采用FLUENT系統(tǒng)展開數(shù)值模擬分析。

針對介質(zhì)中存在的氣、油和液態(tài)水多相流按mixture模型分析，有如下控制方程。

體積守恒：

式中 Vn——多相流總體積，m3；

Vm——各相介質(zhì)體積，m3。

n取3，m=1，2，3時(shí)分別表示氣、油和水三相，下同。

各相相率：

相間黏度不一，不同類型相態(tài)之間存在界面并發(fā)生相對滑移，則有：

這表示相間的相對速度為該相介質(zhì)單獨(dú)通過流場區(qū)域的流速差值，m，l表示不同相態(tài)的流體。

同時(shí)，介質(zhì)還需滿足流動(dòng)控制方程［10-11］：

式中 ρm——密度，kg/m3；

x——幾何坐標(biāo)，m；

i，j——方向；

u——流速，m/s；

p——壓力，Pa；

τ——剪切應(yīng)力，Pa；

Fi——體積力，N。

針對過程工業(yè)領(lǐng)域的湍流問題，RANS（雷諾平均Navier-Stokes）模型即可以獲得足夠準(zhǔn)確結(jié)果并節(jié)約大量計(jì)算資源，因此，即采用RANS模型中的RNG k-ε方程模擬湍流：

其中湍流動(dòng)能k方程：

式中 k——湍流動(dòng)能，m2/s2；

Gk——由層流速度梯度而產(chǎn)生的湍流動(dòng)能，m2/s2；

Gb——由浮力產(chǎn)生的湍流動(dòng)能，m2/s2；

ε——湍流耗散率，m2/s3；

YM——在可壓縮湍流中，擴(kuò)散產(chǎn)生的波動(dòng)，m。

湍流動(dòng)能k可反映湍流脈動(dòng)的劇烈程度。

湍流耗散率ε方程：

式中 C1ε，C2ε，C3ε——常數(shù)，C1ε=1.44，C2ε=1.92，C3ε=0.99。

ε可反映湍流的尺度，但依賴于經(jīng)驗(yàn)常數(shù)。因此，k-ε為半定量半經(jīng)驗(yàn)公式。

空冷器受風(fēng)扇降溫使介質(zhì)冷卻為等壓放熱過程［11］，單位體積介質(zhì)的內(nèi)能：

式（8）中第二項(xiàng)為零，則內(nèi)能只包括焓值和動(dòng)能。

流體能量方程為：

式（9）左邊表示流體內(nèi)能、動(dòng)能及勢能等對時(shí)間和空間尺度的偏導(dǎo)數(shù)，右邊則包括熱傳導(dǎo)、流體輸運(yùn)和黏性阻力產(chǎn)生的熱量。

2.2 幾何建模及邊界條件

16臺(tái)空冷器中A101M曾發(fā)生過管束腐蝕泄漏，因此針對其展開研究，該設(shè)備腐蝕泄漏時(shí)的狀況如圖3所示。

圖3 A101M的腐蝕概況Fig.3 Corrosion status of air cooler A101M

泄漏管束位于入口管箱處因此僅針對此范圍展開討論，空冷管束長度約4 200 mm，接近管束總長的1/2。采用六面體網(wǎng)格進(jìn)行網(wǎng)格劃分，并對管束壁面和管束-管箱連接部位進(jìn)行加密。無關(guān)性驗(yàn)證時(shí)分別采用 2.5×106、4×106及 6.2×106的網(wǎng)格數(shù)量，網(wǎng)格數(shù)為4×106時(shí)相對誤差近1.5%，因此視作達(dá)到了網(wǎng)格無關(guān)性要求。流場幾何模型與網(wǎng)格如圖4所示。

圖4 A101M流場幾何模型及網(wǎng)格Fig.4 Geometric model and meshes of flow field of air cooler A101M

根據(jù)企業(yè)DCS數(shù)據(jù)展開工藝模擬得出介質(zhì)物性參數(shù)及進(jìn)口邊界條件。其中進(jìn)口設(shè)置為速度進(jìn)口，進(jìn)口1，2處流速存在差距，具體參數(shù)見表1。

表1 介質(zhì)物性參數(shù)及進(jìn)口邊界條件Tab.1 Medium physical parameters and inlet boundary conditions

各管束末端設(shè)置壓力出口，湍流強(qiáng)度為1%，進(jìn)出口水力直徑分別按150 mm及19 mm；壁面流動(dòng)條件為固定型，無滑移。針對多相流和湍流分別按式（1）～（9）求解。同時(shí)，因?yàn)橐M空冷器與外界空氣間的對流換熱，又采用能量方程（8）（9）按壁面對流耦合換熱分析。管束壁面設(shè)置對流耦合條件，傳熱系數(shù)按 haircoolerwall≈ 20.5 W/（m2·K）［12］。

3 分析結(jié)果

傳熱過程多由不穩(wěn)定趨向穩(wěn)定，因此，在數(shù)值模擬時(shí)采用非穩(wěn)態(tài)一階隱式時(shí)間推進(jìn)，流動(dòng)控制方程采用Coupled（耦合）算法，對流項(xiàng)等則采用了一階迎風(fēng)格式。設(shè)置殘差為10-5，當(dāng)計(jì)算結(jié)果低于這一水平時(shí)達(dá)到收斂，此時(shí)迭代時(shí)間達(dá)到6 s。

3.1 流動(dòng)特征

3.1.1 時(shí)程變化

多相流自進(jìn)口流入管箱并在管箱中發(fā)生作用的機(jī)制屬于射流卷吸［13］。以下首先按1 s時(shí)刻的全場流線和管箱內(nèi)的湍流動(dòng)能分別展開討論，如圖5所示。

圖5 1 s時(shí)刻的流場特征Fig.5 Characteristic of flow field at the moment of 1 s

瞬時(shí)流線如圖5（a）所示，當(dāng)時(shí)間為1 s時(shí)流速水平最高可達(dá)29 m/s，進(jìn)口1，2射流核心區(qū)線型，流速低于5.8 m/s且保持較低水平，而位于兩側(cè)底部的管束，尤其進(jìn)口2一側(cè)管束因?yàn)楦拷虚g位置流場受到壁面擠壓，因此數(shù)值較高。而管箱中介質(zhì)的湍動(dòng)能則數(shù)值較低，僅在P1排部分管束位置存在湍動(dòng)，介質(zhì)在管箱內(nèi)大范圍尚未形成明顯的介質(zhì)脈動(dòng)。進(jìn)一步考察1 s時(shí)刻管箱內(nèi)部各相介質(zhì)分布，如圖6所示。

如圖6所示，1 s時(shí)管箱中的介質(zhì)多為氣體，氣體介質(zhì)占據(jù)了管箱中大部分空間，且在來流促進(jìn)下不斷向兩側(cè)以擴(kuò)張，氣相介質(zhì)層間阻力不大，只有當(dāng)介質(zhì)由管箱進(jìn)入管束時(shí)會(huì)明顯受阻，因此，1 s時(shí)湍流動(dòng)能高值區(qū)域多集中在管束附近。

圖7～9示出了2，4，6 s時(shí)刻的瞬時(shí)流線及湍動(dòng)能。

圖6 1 s時(shí)刻熱高分氣三相介質(zhì)分布Fig.6 Distribution of hot high pressure separation gas at the moment of 1 s

圖7 2 s時(shí)刻的流場特征Fig.7 Characteristic of flow field at the moment of 2 s

圖8 4 s時(shí)刻的流場特征Fig.8 Characteristic of flow field at the moment of 4 s

圖9 6 s時(shí)刻的流場特征Fig.9 Characteristic of flow field at the moment of 6 s

圖 7（a）～9（a）分析了第 2，4，6 s時(shí)刻的瞬時(shí)流線。2～4 s瞬時(shí)流線所得出的最高流速保持下降的趨勢，4～6 s又有所升高，大量流線也在2～4 s時(shí)不斷發(fā)生糾纏并出現(xiàn)交叉，同時(shí)，漩渦在4 s時(shí)刻開始形成，隨后擴(kuò)大，6 s則逐漸縮小，之前發(fā)生交織的流線也漸漸稀疏。相應(yīng)地，湍流動(dòng)能值也發(fā)生著明顯的變化，如圖7（b）～9（b）。其中，2 s時(shí)湍流動(dòng)能分布情況與1 s時(shí)接近，4，6 s時(shí)各管束尤其P2排管束區(qū)域的湍動(dòng)能高值區(qū)逐漸相連，但最高值會(huì)不斷降低。湍動(dòng)能逐漸下降標(biāo)志著流場脈動(dòng)程度逐漸降低，但湍動(dòng)能高值區(qū)的擴(kuò)展則意味著越來越多的介質(zhì)在發(fā)生湍動(dòng)。之所以存在這一現(xiàn)象，是因?yàn)閬砹魅肷浜笫紫葘⒉粩鄾_刷管箱、管束壁面隨后再同設(shè)備中已經(jīng)存在的介質(zhì)發(fā)生黏阻作用，隨后先前進(jìn)入管箱的介質(zhì)又被來流卷進(jìn)入射區(qū)域，從而形成了入射-卷吸機(jī)制耗散了介質(zhì)動(dòng)能，而又帶動(dòng)了越來越大范圍介質(zhì)的脈動(dòng)。自2～6 s時(shí)刻，多相流中氣、油兩相分布變化不大，而水相分布變化明顯，如圖10所示。

圖10 不同時(shí)刻管箱內(nèi)的氣、水組分等值線Fig.10 Contour of gaseous phase and liquid-water in the channel box at different moments

據(jù)能量方程，湍流動(dòng)能的變化率實(shí)際上與黏度v成反比，因此，當(dāng)液態(tài)水介質(zhì)含量較高的區(qū)域不斷擴(kuò)展時(shí)多相流湍動(dòng)能的高值區(qū)位置也會(huì)不斷擴(kuò)大。另外，發(fā)生過腐蝕泄漏的P2排管，氣相組分含量較高，在壓力不變的前提下，腐蝕性氣體NH3、HCl等含量更高易反應(yīng)生成銨鹽，參照API 932B中關(guān)于銨鹽腐蝕的評(píng)價(jià)細(xì)則，這一區(qū)域確有較高的銨鹽流動(dòng)腐蝕風(fēng)險(xiǎn)。

3.1.2 空間分布

6 s時(shí)刻全場管束的平均流速分布如圖11所示。

圖11 6 s時(shí)全場管束的平均流速Fig.11 Average velocities in all tube bundles at the moment of 6 s

P1-2排管束的流速在不同管束處存在明顯波動(dòng)，其中P1排管束的算術(shù)平均值為3.11 m/s，P2排為2.93 m/s。API 932B認(rèn)為，當(dāng)設(shè)備中介質(zhì)流速低于3 m/s時(shí)銨鹽結(jié)晶并腐蝕管束可能性較大，因此，流速水平也反映了有關(guān)管束區(qū)域較高的腐蝕風(fēng)險(xiǎn)［14-15］。

若按進(jìn)口1，2接管中軸線為參考位置，軸線正下方屬于射流核心區(qū)，這部分區(qū)域的湍流動(dòng)能耗散與生成可達(dá)到平衡，介質(zhì)將直接沖擊底板令動(dòng)能大幅下降，因此P1-10、P2-9等管束流速較低；接管兩翼的區(qū)域流場則會(huì)因?yàn)槿肷?卷吸機(jī)制不斷擴(kuò)大湍動(dòng)范圍，從而令各管流速不斷波動(dòng)，并在介質(zhì)沖刷左右側(cè)管箱壁面后急遽下降。

綜上，因不同管束位置而發(fā)生波動(dòng)的流速是進(jìn)口1、2射流在管箱中形成卷吸協(xié)同作用機(jī)制導(dǎo)致的結(jié)果，也是湍流隨機(jī)性及耗散性的表現(xiàn)，隨著數(shù)值模擬結(jié)果達(dá)到收斂，各管束中介質(zhì)流速水平就已經(jīng)達(dá)到“穩(wěn)定”的波動(dòng)狀態(tài)，根據(jù)API 932B避免管束發(fā)生銨鹽流動(dòng)腐蝕需令來流速度更高，并保持均衡［16］。

3.2 溫度場結(jié)果

3.2.1 時(shí)程變化

6 s內(nèi)多相流的溫度場變化如圖12所示。

圖12 1～6 s時(shí)刻的溫度場Fig.12 Temperature field results at the moments of 1 s to 6 s

溫度場最高值一直在125 ℃左右，但最低溫度卻隨時(shí)間不斷升高。另一方面，70 ℃以上區(qū)域也在逐漸擴(kuò)展。各管束溫度在介質(zhì)流出方向呈現(xiàn)波浪式分布。存在這一現(xiàn)象，與介質(zhì)流速不一有關(guān)，介質(zhì)流出較快則換熱時(shí)間較短，溫度較高，反之溫度較低。

3.2.2 空間分布及熱-流耦合分析

在6 s時(shí)刻，按不同管束位置考察平均溫度及流速如圖13所示。

圖13 5 s時(shí)刻全場管束的溫度與流速統(tǒng)計(jì)Fig.13 Temperature and velocity results of all tube bundles at the moment of 5 s

由圖13可見，P1，P2排管束溫度高低起伏同流速較為一致。考察2個(gè)參量相關(guān)系數(shù)：

式中 ti——每個(gè)管束的溫度，℃；

ui——每個(gè)管束的流速，m/s；

P1排管相關(guān)系數(shù)rtu約為0.75，P2排則可達(dá)0.85。根據(jù)式（9），左邊流場能量的空間變化率等于右邊的溫度二階偏導(dǎo)和其他能量值，脈動(dòng)程度較低時(shí)，能量的耗散項(xiàng)和輸運(yùn)項(xiàng)等基本可以忽略，右邊近似等于介質(zhì)傳熱量偏導(dǎo)，而左邊的流體能量包括速度項(xiàng)，因此，速度空間偏導(dǎo)同傳熱量偏導(dǎo)存在較高相關(guān)性，多相流介質(zhì)速度同溫度高度相關(guān)。

4 結(jié)論

（1）數(shù)值模擬時(shí)考慮非穩(wěn)態(tài)過程中不同時(shí)間流態(tài)的變化，隨著時(shí)間推移多相流介質(zhì)在射流卷吸機(jī)制的作用下，脈動(dòng)趨于緩和，流速、湍動(dòng)能值都將逐步降低。

（2）位于下方的P2排管束存在較高的氣相組分比例和較多的液態(tài)水，同時(shí)，流速平均值也低于API 932B 規(guī)定的下限，因此存在較大的銨鹽流動(dòng)風(fēng)險(xiǎn)。

（3）隨著多相流介質(zhì)湍動(dòng)趨緩，能量方程中的黏性耗散項(xiàng)也會(huì)減小，因此，各管流速水平同溫度具有較高的相關(guān)性。工程實(shí)際中可參照不同位置管束紅外測溫結(jié)果，估算流速，科學(xué)預(yù)測并有效防范空冷管束銨鹽流動(dòng)腐蝕。