加氫空冷多相流熱耦合數值分析

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

0 引言

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

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

1 工況及設備

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

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

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

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

2 數理模型及邊界條件

2.1 數理模型

采用FLUENT系統展開數值模擬分析。

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

體積守恒：

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

Vm——各相介質體積，m3。

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

各相相率：

相間黏度不一，不同類型相態之間存在界面并發生相對滑移，則有：

這表示相間的相對速度為該相介質單獨通過流場區域的流速差值，m，l表示不同相態的流體。

同時，介質還需滿足流動控制方程［10-11］：

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

x——幾何坐標，m；

i，j——方向；

u——流速，m/s；

p——壓力，Pa；

τ——剪切應力，Pa；

Fi——體積力，N。

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

其中湍流動能k方程：

式中 k——湍流動能，m2/s2；

Gk——由層流速度梯度而產生的湍流動能，m2/s2；

Gb——由浮力產生的湍流動能，m2/s2；

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

YM——在可壓縮湍流中，擴散產生的波動，m。

湍流動能k可反映湍流脈動的劇烈程度。

湍流耗散率ε方程：

式中 C1ε，C2ε，C3ε——常數，C1ε=1.44，C2ε=1.92，C3ε=0.99。

ε可反映湍流的尺度，但依賴于經驗常數。因此，k-ε為半定量半經驗公式。

空冷器受風扇降溫使介質冷卻為等壓放熱過程［11］，單位體積介質的內能：

式（8）中第二項為零，則內能只包括焓值和動能。

流體能量方程為：

式（9）左邊表示流體內能、動能及勢能等對時間和空間尺度的偏導數，右邊則包括熱傳導、流體輸運和黏性阻力產生的熱量。

2.2 幾何建模及邊界條件

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

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

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

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

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

表1 介質物性參數及進口邊界條件Tab.1 Medium physical parameters and inlet boundary conditions

各管束末端設置壓力出口，湍流強度為1%，進出口水力直徑分別按150 mm及19 mm；壁面流動條件為固定型，無滑移。針對多相流和湍流分別按式（1）～（9）求解。同時，因為要模擬空冷器與外界空氣間的對流換熱，又采用能量方程（8）（9）按壁面對流耦合換熱分析。管束壁面設置對流耦合條件，傳熱系數按 haircoolerwall≈ 20.5 W/（m2·K）［12］。

3 分析結果

傳熱過程多由不穩定趨向穩定，因此，在數值模擬時采用非穩態一階隱式時間推進，流動控制方程采用Coupled（耦合）算法，對流項等則采用了一階迎風格式。設置殘差為10-5，當計算結果低于這一水平時達到收斂，此時迭代時間達到6 s。

3.1 流動特征

3.1.1 時程變化

多相流自進口流入管箱并在管箱中發生作用的機制屬于射流卷吸［13］。以下首先按1 s時刻的全場流線和管箱內的湍流動能分別展開討論，如圖5所示。

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

瞬時流線如圖5（a）所示，當時間為1 s時流速水平最高可達29 m/s，進口1，2射流核心區線型，流速低于5.8 m/s且保持較低水平，而位于兩側底部的管束，尤其進口2一側管束因為更靠近中間位置流場受到壁面擠壓，因此數值較高。而管箱中介質的湍動能則數值較低，僅在P1排部分管束位置存在湍動，介質在管箱內大范圍尚未形成明顯的介質脈動。進一步考察1 s時刻管箱內部各相介質分布，如圖6所示。

如圖6所示，1 s時管箱中的介質多為氣體，氣體介質占據了管箱中大部分空間，且在來流促進下不斷向兩側以擴張，氣相介質層間阻力不大，只有當介質由管箱進入管束時會明顯受阻，因此，1 s時湍流動能高值區域多集中在管束附近。

圖7～9示出了2，4，6 s時刻的瞬時流線及湍動能。

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

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

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

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

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

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

據能量方程，湍流動能的變化率實際上與黏度v成反比，因此，當液態水介質含量較高的區域不斷擴展時多相流湍動能的高值區位置也會不斷擴大。另外，發生過腐蝕泄漏的P2排管，氣相組分含量較高，在壓力不變的前提下，腐蝕性氣體NH3、HCl等含量更高易反應生成銨鹽，參照API 932B中關于銨鹽腐蝕的評價細則，這一區域確有較高的銨鹽流動腐蝕風險。

3.1.2 空間分布

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

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

P1-2排管束的流速在不同管束處存在明顯波動，其中P1排管束的算術平均值為3.11 m/s，P2排為2.93 m/s。API 932B認為，當設備中介質流速低于3 m/s時銨鹽結晶并腐蝕管束可能性較大，因此，流速水平也反映了有關管束區域較高的腐蝕風險［14-15］。

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

綜上，因不同管束位置而發生波動的流速是進口1、2射流在管箱中形成卷吸協同作用機制導致的結果，也是湍流隨機性及耗散性的表現，隨著數值模擬結果達到收斂，各管束中介質流速水平就已經達到“穩定”的波動狀態，根據API 932B避免管束發生銨鹽流動腐蝕需令來流速度更高，并保持均衡［16］。

3.2 溫度場結果

3.2.1 時程變化

6 s內多相流的溫度場變化如圖12所示。

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

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

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

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

圖13 5 s時刻全場管束的溫度與流速統計Fig.13 Temperature and velocity results of all tube bundles at the moment of 5 s

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

式中 ti——每個管束的溫度，℃；

ui——每個管束的流速，m/s；

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

4 結論

（1）數值模擬時考慮非穩態過程中不同時間流態的變化，隨著時間推移多相流介質在射流卷吸機制的作用下，脈動趨于緩和，流速、湍動能值都將逐步降低。

（2）位于下方的P2排管束存在較高的氣相組分比例和較多的液態水，同時，流速平均值也低于API 932B 規定的下限，因此存在較大的銨鹽流動風險。

（3）隨著多相流介質湍動趨緩，能量方程中的黏性耗散項也會減小，因此，各管流速水平同溫度具有較高的相關性。工程實際中可參照不同位置管束紅外測溫結果，估算流速，科學預測并有效防范空冷管束銨鹽流動腐蝕。