核電機組二回路碳鋼給水管道節(jié)流孔板下游流動加速腐蝕數(shù)值模擬研究

宗凌風(fēng)， 周克毅， 黃軍林， 司曉東， 徐青藍

(1. 東南大學(xué) 能源熱轉(zhuǎn)換及其過程測控教育部重點實驗室， 南京 210096；2. 中國核動力研究設(shè)計院 核反應(yīng)堆系統(tǒng)設(shè)計技術(shù)重點實驗室， 成都 610213)

核電機組二回路的汽水輸運管道大量使用碳鋼，流動加速腐蝕(FAC)是影響管道安全性的重要因素之一[1]。FAC過程中，管道內(nèi)壁保護性氧化膜被流動工質(zhì)加速分解和破壞，無法有效抑制管壁基體的持續(xù)氧化[2]。管道持續(xù)減薄，并可能在機組運行過程中突然破裂，嚴重影響機組安全經(jīng)濟運行[3]。

FAC問題出現(xiàn)的早期，學(xué)者將其歸因于單純的機械沖刷。采取增加管徑和減少彎角等措施后，F(xiàn)AC問題得到了緩解，之后一段時間內(nèi)FAC問題未受到關(guān)注。1986年美國Surry核電站發(fā)生的管道破裂事故是FAC研究的轉(zhuǎn)折點。事故中，給水泵下游彎管在機組運行時突然破裂，造成4死4傷的悲劇[4]。FAC致使局部管壁嚴重減薄是造成該起事故的主因。對FAC的研究自此受到高度重視，美國電力研究院(RPRI)、法國電力公司(EDF)、德國西門子公司(Siemens/KWU)等隨即開展相關(guān)研究。然而，2004年和2007年，日本Mihama核電站和美國Iatan核電站又接連發(fā)生了嚴重的管道破裂事故，再次造成人員傷亡的悲劇[5]。分析顯示，節(jié)流孔板及閥門下游的FAC是這些事故的主因。因此，迫切需要深入探究節(jié)流孔板下游的FAC機理，并提出有效抑制和預(yù)防節(jié)流孔板下游FAC的方法[6]。

在彎管、節(jié)流孔板和T形三通管等流場劇烈變化的奇點下游，最易發(fā)生FAC事故[7]。工質(zhì)流經(jīng)這些管道部件時，流向的劇烈變化及二次流不穩(wěn)定性發(fā)展等因素，致使FAC速率顯著增加[8]。CRAWFORD N M等[9]發(fā)現(xiàn)孔板和閥門的結(jié)構(gòu)特性促使下游湍流靠近壁面，提高了壁面的傳質(zhì)速率。

開展FAC的試驗研究需要搭建高溫高壓循環(huán)回路，并且對循環(huán)工質(zhì)的pH、溫度和溶氧[10]等有嚴格要求。隨著計算流體動力學(xué)(CFD)的不斷發(fā)展，其計算能力可以滿足傳質(zhì)和腐蝕研究的要求，而精確的近壁網(wǎng)格和合理的近壁湍流模型可保證壁面?zhèn)髻|(zhì)數(shù)據(jù)的計算準確性。學(xué)者采用CFD方法模擬了節(jié)流孔板下游流場并分析了孔板形狀對下游FAC的影響。YIN G等[11]結(jié)合k-ω剪應(yīng)力輸運湍流模型，模擬了孔板厚度和孔徑比對孔板下游壓力、流量和渦量分布的影響。AHMED W H等[12]采用模擬和試驗驗證，研究孔徑比為0.25、0.50和0.75的孔板下游FAC，并用氫氧化鈣制成試驗段驗證了數(shù)值結(jié)果。

國內(nèi)外學(xué)者對直管、彎管等的FAC問題已開展了較為深入的研究，對管內(nèi)流場影響FAC的機理有了一定的認識。相比之下，雖然對流場更為復(fù)雜的節(jié)流孔板下游FAC問題也有了一定的研究，但還未能形成完整的體系。部分學(xué)者針對孔板孔徑比對FAC的影響展開了研究，但對節(jié)流孔板的倒角角度、厚度等幾何結(jié)構(gòu)因素影響近壁流場進而影響FAC的具體機制仍不明確，亦未在此基礎(chǔ)上提出有效的緩解FAC的具體措施。筆者基于Sanchez-Caldera模型，采用數(shù)值模擬的方法，研究節(jié)流孔板下游FAC速率和位置的分布。在此基礎(chǔ)上，重點分析了不同入口流速和倒角角度對節(jié)流孔板下游FAC分布的影響，以期為設(shè)計人員進行核電機組二回路給水管道節(jié)流孔板的設(shè)計選型提供參考。

1 FAC機理和模型

1.1 FAC機理

圖1是FAC過程，即氧化膜形成和溶解的機理，其中δ為氧化膜的厚度，δy為傳質(zhì)邊界層的厚度，δd為飽和層的厚度。FAC主要分為3個過程[13]：(1)金屬基體表面失去電子生成Fe2+；(2)一部分Fe2+與水發(fā)生水合反應(yīng)生成Fe3O4，另一部分Fe2+通過氧化膜的空隙傳遞到主體溶液中；(3)氧化膜進一步溶于水并擴散到主體溶液中。

圖1 氧化膜形成和溶解機理

第1個過程發(fā)生在基體/氧化膜界面處。在陽極，F(xiàn)e失去電子生成Fe2+，發(fā)生氧化反應(yīng)；在還原性溶液中，陰極處H2O得到電子，生成OH-和H2，發(fā)生還原反應(yīng)。結(jié)合陽極和陰極的反應(yīng)，第1個過程的電化學(xué)反應(yīng)為：

(1)

此時，F(xiàn)e2+和Fe(OH)2的濃度可以根據(jù)反應(yīng)平衡常數(shù)獲得。

第2個過程的反應(yīng)是Fe(OH)2通過水合反應(yīng)生成Fe3O4，即

(2)

該過程生成了底層致密、外層疏松的多孔氧化膜層[14]，剩余的可溶性含F(xiàn)e組分通過氧化膜/主體溶液界面擴散到主體溶液中，并被不斷流動的主體溶液帶走，使得可溶性含F(xiàn)e組分的溶解過程持續(xù)進行。

第3個過程，氧化膜在氧化膜/水界面處溶解，生成可溶性含F(xiàn)e組分[15]。該反應(yīng)過程由多個過程組合而成，可總結(jié)為以下化學(xué)方程式：

(3)

式中：b=0、1、2和3。

FAC過程由電化學(xué)腐蝕過程和傳質(zhì)過程共同控制。當(dāng)主體溶液流速增加時，傳質(zhì)過程速度加快，體系內(nèi)濃度差增大，基體Fe和氧化膜的活化和溶解反應(yīng)速率加快，進而加快了FAC速率。

1.2 FAC預(yù)測模型

基于FAC機理，結(jié)合Sanchez-Caldera[16]模型，可獲得FAC速率計算式為：

(4)

式中：KFAC為FAC速率，mm/a；Ceq為可溶性含F(xiàn)e組分的溶解度，mol/L；C∞為主體溶液中可溶性含F(xiàn)e組分的濃度，mol/L；K*為反應(yīng)速率常數(shù)，m/s；k為傳質(zhì)系數(shù)，m/s；f為基體氧化膜界面Fe2+轉(zhuǎn)化為氧化膜的比例；θ為氧化膜的孔隙率；DFe為可溶性含F(xiàn)e組分在水中的擴散系數(shù)，m2/s。

BERGE P H等[17]認為，有一半的可溶性含F(xiàn)e組分轉(zhuǎn)化為氧化膜，因此f=0.5；主體溶液中，可溶性含F(xiàn)e組分的濃度很低，因此C∞可忽略；活化反應(yīng)速率常數(shù)遠大于氧化層向主體溶液的傳質(zhì)系數(shù)，因此1/K*亦可忽略。式(4)可簡化為：

(5)

當(dāng)流體溫度T低于100 ℃時，δ約為5 nm；當(dāng)T高于120 ℃時，δ約為1 μm。當(dāng)T≤422 K時，θ=0.05%；當(dāng)422 K

1.3 溶解度模型

高溫低氧條件下，基體和溶液界面生成的Fe3O4氧化膜具有微溶性和滲透性，在主體溶液流過時會不斷溶解。Fe在溶液中的溶解度可表示為Fe的可溶性組分的溶解度之和，即

Ceq=[Fe2+]+[Fe(OH)+]+[Fe(OH)2]+

(6)

式(3)可以分解為以下氧化物在氧化膜/主體溶液界面的溶解反應(yīng)：

(7)

上述方程的反應(yīng)平衡常數(shù)記為Kx(x=0、1、2、3)，則

(8)

因此，式(6)可簡化為：

Ceq=(K010-2m+K110-m+K2+

K310m)[H2]1/3

(9)

式中：m為含F(xiàn)e溶液的pH；[H2]為7.80×10-4～3.25×10-3mol/L[18]。

式(7)的化學(xué)平衡常數(shù)可表示為：

(10)

式中：R為理想氣體常數(shù)，J/(mol·K)；A、B、C均為常數(shù)，其取值見表1[16]。

表1 參數(shù)A、B和C的取值

1.4 傳質(zhì)系數(shù)模型

由FAC機理可知，傳質(zhì)系數(shù)是表征FAC速率的重要參數(shù)之一。試驗研究中，常用3個無因次準則數(shù)——雷諾數(shù)(Re)、舍伍德數(shù)(Sh)和施密特數(shù)(Sc)的組合形式來表示傳質(zhì)系數(shù)[19-20]：

(11)

式中：kFe為Fe的傳質(zhì)系數(shù)，m/s；d為特征長度，在管道內(nèi)部，d的取值是管道的內(nèi)徑，m。

通過ANSYS軟件求解連續(xù)性方程、動量方程和能量方程，計算孔板下游流速分布，基于流速和Re的關(guān)系式求出Re的分布，并代入式(11)即可得到孔板下游傳質(zhì)系數(shù)的分布。同時，計算采用如下假設(shè)：流動為定常不可壓縮流動；忽略重力的影響；不考慮管子的散熱。

2 計算結(jié)果與分析

2.1 溶解度

水化學(xué)條件選取核電機組二回路管道運行的溫度壓力范圍，結(jié)合FAC高發(fā)的溫度范圍。具體取值如下：溫度為150 ℃，壓力為0.95 MPa，常溫條件(25 ℃)下流體的pH為9.4，150 ℃條件下流體的pH為6.86，H2濃度為1.5×10-3mol/L。將上述取值代入式(9)，可以得到可溶性含F(xiàn)e組分的飽和溶解度為6.95×10-8mol/L。

2.2 流場

孔板的幾何模型見圖2，管道內(nèi)徑d為50 mm，孔板直徑d0為25 mm，孔板厚度為3 mm，孔板上游長度為10d，孔板下游長度為20d，即上游長度為500 mm，下游長度為1 000 mm，流體到達孔板時已經(jīng)充分發(fā)展。將孔板出口平面作為坐標原點，則此處z/d=0(z為孔板下游到孔板的距離)，孔徑比d0/d=0.5，孔板出水側(cè)倒角角度β分別設(shè)置為15°、30°、45°、60°和90°。當(dāng)?shù)菇墙嵌葹?0°時，孔板出水側(cè)不存在倒角，即倒角角度為0°。

圖2 孔板幾何模型

在計算前處理軟件ICEM中對幾何模型進行網(wǎng)格劃分，采取O形剖分的方法劃分結(jié)構(gòu)性網(wǎng)格，對孔板的倒角進行二次剖分，并對孔板倒角和孔板壁面進行局部網(wǎng)格加密處理。為保證計算精度，第一層網(wǎng)格的高度為0.3 mm，網(wǎng)格總數(shù)量為360萬，網(wǎng)格質(zhì)量分布在0.7～1.0，符合FLUENT軟件求解器的計算要求。網(wǎng)格劃分和孔板網(wǎng)格加密情況見圖3。

圖3 網(wǎng)格劃分

在FLUENT軟件求解器中，采用基于壓力的求解方式，選用Realizablek-ε湍流模型，標準壁面函數(shù)，邊界條件采用速度入口、壓力出口。入口流速分別為0.8 m/s、2.0 m/s和5.0 m/s，管內(nèi)壓力為0.95 MPa。求解算法采用速度壓力耦合的SIMPLE算法；動量方程、湍動能、湍動能耗散和能量方程采用二階迎風(fēng)離散格式，松弛因子采用默認值，收斂標準為1×10-4，對該模型進行數(shù)值計算。

將倒角角度為45°、入口流速為0.8 m/s的數(shù)值模擬結(jié)果和實驗結(jié)果進行對比，結(jié)果見圖4。由圖4(a)可知：數(shù)值模擬結(jié)果較好地符合了壁面壓力經(jīng)過孔板時陡降的情況，與試驗數(shù)據(jù)[12]符合較好。由圖4(b)可知：數(shù)值模擬結(jié)果顯示了流體在經(jīng)過孔板時軸線速度陡增的現(xiàn)象，與試驗結(jié)果一致。因此，可認為模擬的管內(nèi)流場是準確的。

圖4 孔板下游壁面壓力和軸線速度分布

圖5為倒角角度為45°時，入口流速為2.0 m/s時孔板下游橫截面和管道中心截面的流場分布圖。當(dāng)流體流經(jīng)孔板時，由于流道瞬間縮小，根據(jù)質(zhì)量守恒定律，流體的流速增大；在通過倒角后，流道恢復(fù)到初始的橫截面積。在孔板后流體發(fā)生分離的現(xiàn)象，在兩側(cè)孔板后形成對稱的漩渦。由于兩側(cè)漩渦的存在，管道中心位置高速流體流通區(qū)域縮小，流出倒角后的流體中心速度繼續(xù)增大；分離流體后和壁面重新附著，軸線速度也在達到最大值以后逐漸減小，最后趨于穩(wěn)定。隨著與中心處高速流體發(fā)生動量交換，管道外圍流體的流速也逐漸增大，最后趨于穩(wěn)定。

圖5 入口流速為2.0 m/s時孔板下游橫截面及管道中心截面流場分布圖

由圖5可知：管道下游流場呈中心對稱分布，在z/d=0時管道外圍流場分布不均勻，隨后到z/d=0.8時趨于穩(wěn)定。管內(nèi)漩渦從中心位置向外圍逐漸擴散，在z/d=4.0左右時，漩渦消失，流體重新附著。管道軸線處的高速流通區(qū)域面積先增大后減小，隨著流體的流動，管道中心流體與外圍流體發(fā)生動量交換，流速先增大后減小。

圖6是入口流速為0.8 m/s時不同倒角角度孔板下游的流場分布圖。

圖6 不同倒角孔板下游流場分布

圖6中截取了孔板下游-0.5d～5.5d流場變化比較劇烈的局部放大圖。由圖6可知：當(dāng)?shù)菇墙嵌葹?5°和90°時，孔板下游的流場分布較為相似；而倒角角度為45°時，孔板下游流場中軸線位置的高速流體區(qū)域較為狹長，外圍兩側(cè)的漩渦的中心位置基本不變，而末端分離流體的附著點位置向下游移動。當(dāng)?shù)菇墙嵌葹?5°和90°時，孔板下游的流場變化較大，使得中心流體的流速迅速增大并達到峰值，并與周圍流體發(fā)生激烈的動量交換，隨后迅速減小后趨于穩(wěn)定。而倒角角度為45°時，下游的流場變化相對穩(wěn)定，由于與周圍流體的動量交換較少，中心流體的流速在增大到峰值后緩慢降低，流場穩(wěn)定速度較慢。因此，倒角角度會對孔板下游流場產(chǎn)生一定的影響。

圖7是入口流速為0.8 m/s時不同倒角角度孔板下游的速度矢量圖。由圖7可見：孔板外圍的漩渦呈軸對稱旋轉(zhuǎn)分布，且在該漩渦的上游、孔板的正后方還存在一個尺度較小的漩渦。隨著倒角角度的增大，小漩渦逐漸增大。當(dāng)?shù)菇墙嵌葹?5°時，孔板下游的流體流動方向與管道軸線呈現(xiàn)出一定的角度，這減緩了孔板的存在對流場的影響，因此下游流體平均流速較小，流場變化較為緩慢。而當(dāng)?shù)菇墙嵌葹?5°和90°時，孔板下游的流體流動方向與管道軸線夾角較小，表明此時倒角角度對減緩孔板影響的作用不明顯。

圖7 不同倒角角度下孔板下游速度矢量分布圖

2.3 FAC速率

取入口流速為0.8 m/s、2.0 m/s和5.0 m/s，分別計算出倒角角度為15°、30°、45°、60°和90° 5種孔板下游的流場分布，并計算出Sh。根據(jù)Sh與傳質(zhì)系數(shù)的關(guān)系，計算出傳質(zhì)系數(shù)。

倒角角度為45°時不同入口流速孔板下游的傳質(zhì)系數(shù)分布和入口流速為0.8 m/s時不同倒角角度孔板下游的傳質(zhì)系數(shù)分布見圖8。

圖8 孔板下游傳質(zhì)系數(shù)分布

由圖8可知：當(dāng)?shù)菇墙嵌纫欢〞r，孔板下游傳質(zhì)系數(shù)隨入口流速的增大而增大；當(dāng)入口流速一定時，倒角角度從15°變化到90°時，傳質(zhì)系數(shù)整體先減小后增大。倒角角度從30°變化到45°再到60°時的變化率遠小于從15°變化到30°以及從60°變化到90°時的變化率。由于流體入口流速的增大，軸線位置的流體在相同時間內(nèi)流到了孔板下游更遠的位置，同時加強自身與周圍流體的動量交換，傳質(zhì)系數(shù)增大，并且峰值向遠離孔板的位置偏移。當(dāng)流體入口流速一定時，倒角角度為45°能減緩孔板對流場帶來的變化，下游流場的變化較為平穩(wěn)，孔口的流體能更好地和周圍流體發(fā)生動量交換，因此下游的動量交換減小，傳質(zhì)系數(shù)減小。當(dāng)?shù)菇墙嵌刃∮诨蛘叽笥?5°時，孔板下游流場變化增大，孔口的動量交換減小，下游的動量交換加劇，因此下游傳質(zhì)系數(shù)增加。倒角角度的改變對傳質(zhì)系數(shù)峰值位置的影響不大。隨著倒角角度的增大，傳質(zhì)系數(shù)峰值位置遠離孔板。

圖9為入口流速為0.8 m/s、倒角角度為45°時孔板下游FAC速率分布模擬結(jié)果和相同條件下與文獻[12]試驗結(jié)果的對比圖。由圖9可知，F(xiàn)AC速率的數(shù)值模擬結(jié)果和試驗結(jié)果符合較好。與圖8所示的傳質(zhì)系數(shù)分布相似，F(xiàn)AC速率在孔板下游急劇增加并且在z/d=1.5～2.5處達到峰值；隨著下游流動的發(fā)展，F(xiàn)AC速率緩慢降低，直至趨于穩(wěn)定。

圖9 FAC速率的數(shù)值模擬結(jié)果與試驗結(jié)果對比

圖10顯示了倒角角度為15°、45°和90°時，孔板下游在3種入口流速下的FAC速率。

圖10 不同流速下孔板下游FAC速率模擬結(jié)果

由圖10可知，孔板下游FAC速率隨著入口流速的增大而顯著增加。此外，當(dāng)入口流度從0.8 m/s變化到2.0 m/s時，F(xiàn)AC速率峰值位置無明顯變化；而當(dāng)入口流速度從2.0 m/s變化到5.0 m/s時，F(xiàn)AC速率峰值位置明顯向下游移動。這是因為隨著流量增大，在流體慣性作用下，動量交換最劇烈的部位向遠離孔板的位置移動，所以FAC速率峰值位置向下游移動。

圖11為入口流速為0.8 m/s時不同倒角角度孔板下游的FAC速率。

圖11 入口流速為0.8 m/s時孔板下游FAC速率模擬結(jié)果

由圖11可知，隨著倒角角度的增大，F(xiàn)AC速率先減小后增大，倒角角度為45°時FAC速率峰值最小。當(dāng)?shù)菇墙嵌葹?5°時，孔板下游流場變化較為緩慢，孔口中心流體與周圍流體得到充分的動量交換，下游流場變化相對穩(wěn)定。此時，倒角減緩了孔板對流場的影響，因此FAC速率最小。隨著倒角角度的增加，F(xiàn)AC速率峰值位置緩慢向下游移動。對比圖10和圖11，入口流速對FAC速率的影響更顯著。

倒角角度和入口流速均影響了節(jié)流孔板下游FAC速率峰值的大小和位置。表2為FAC速率峰值位置隨入口流速和倒角角度的變化，表3為FAC速率峰值隨入口流速和倒角角度的變化。

表2 FAC速率峰值位置隨入口流速和倒角角度的變化

表3 FAC速率峰值隨入口流速和倒角角度的變化

結(jié)合表2、表3可以看出：當(dāng)?shù)菇墙嵌葹?5°時，F(xiàn)AC的高發(fā)區(qū)位于1.55d～1.70d；當(dāng)?shù)菇墙嵌葹?5°時，F(xiàn)AC的高發(fā)區(qū)位于1.70d～1.85d；當(dāng)?shù)菇墙嵌葹?0°時，F(xiàn)AC的高發(fā)區(qū)位于1.80d～2.00d。當(dāng)?shù)菇墙嵌葟?5°變化到45°時，F(xiàn)AC速率峰值逐漸減小，當(dāng)?shù)菇墙嵌葟?5°變化到90°時，F(xiàn)AC速率峰值逐漸增大。FAC速率峰值位置隨倒角角度的增加向遠離孔板的位置移動。

當(dāng)給定的入口流速一定時，以0.8 m/s為例，當(dāng)?shù)菇墙嵌葟?5°變化到60°時，F(xiàn)AC速率峰值的增量較小，相對變化量為1.18%；當(dāng)?shù)菇墙嵌葟?0°變化到90°時，相對變化量為5.81%。因此，當(dāng)?shù)菇墙嵌仍?5°～60°變化時，F(xiàn)AC速率峰值的變化量小于倒角角度在60°～90°的FAC速率峰值變化量。此外，倒角角度在30°～45°變化時，F(xiàn)AC速率峰值的變化量小于倒角角度在15°～30°的FAC速率峰值變化量。

當(dāng)孔板出水側(cè)存在一定角度的倒角時，隨著倒角角度的改變，流體的流向會發(fā)生一定的改變，導(dǎo)致下游FAC速率峰值位置改變；當(dāng)?shù)菇墙嵌葹?0°時，其下游FAC速率明顯大于倒角角度為30°、45°和60°的孔板；倒角角度為30°、45°和60°的孔板下游的FAC速率比較接近，當(dāng)?shù)菇墙嵌葹?5°時，F(xiàn)AC速率最小。

3 結(jié)語

(1) 基于Sanchez-Caldera模型，借助ANSYS軟件研究了節(jié)流孔板下游FAC速率的大小和位置分布，數(shù)值模擬結(jié)果和文獻試驗結(jié)果數(shù)據(jù)吻合較好。

(2) 對于給定的倒角角度，入口流速越大時，孔板下游的FAC速率越大。當(dāng)入口流速從0.8 m/s變化到2.0 m/s時，F(xiàn)AC速率峰值位置變化不明顯；當(dāng)入口流速從2.0 m/s變化到5.0 m/s時，F(xiàn)AC速率峰值位置明顯向下游移動。

(3) 對于給定的入口流速，倒角角度越小，F(xiàn)AC速率峰值位置越靠近孔板，并且隨倒角角度的增加，F(xiàn)AC速率峰值先減小后增大。倒角角度在30°～45°及45°～60°變化時，下游FAC速率峰值的變化率小于倒角角度在15°～30°及60°～90°的變化率。

(4) 倒角角度為90°時，下游FAC速率最大；當(dāng)?shù)菇墙嵌葹?5°時，F(xiàn)AC速率最小。因此，當(dāng)?shù)菇墙嵌葹?5°時，能夠有效減小下游的FAC速率。