劉 延，劉曉晶，何 輝

(上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240)

反應(yīng)堆堆芯長(zhǎng)期處于高溫、高壓、高放射性的環(huán)境下，燃料性能是評(píng)價(jià)堆芯安全性和經(jīng)濟(jì)性的關(guān)鍵性指標(biāo)。反應(yīng)堆運(yùn)行時(shí)，蒸汽發(fā)生器等持續(xù)受到冷卻劑的沖刷腐蝕，形成大量以鐵、鎳離子及其氧化物為主的氧化腐蝕產(chǎn)物。這些腐蝕產(chǎn)物在過冷沸騰(SNB)的驅(qū)動(dòng)下，沉積在堆芯上部的燃料包殼表面，形成一層薄的疏松多孔的結(jié)垢層，被稱為氧化腐蝕產(chǎn)物沉積層(chalk river unidentified deposit, CRUD)[1]。一方面，由于腐蝕產(chǎn)物的沉積，燃料棒徑向熱阻增大，傳熱被惡化。另一方面，沉積層內(nèi)呈現(xiàn)疏松多孔的結(jié)構(gòu)，SNB被強(qiáng)化，使得硼元素不斷濃縮直至析出并被沉積層吸附，導(dǎo)致硼元素在燃料棒軸向分布不均勻。由于10B具有顯著的吸中子能力，因此沉積層內(nèi)的硼吸附會(huì)改變堆芯軸向的中子分布，使得堆芯功率分布向底部畸變，引發(fā)堆芯功率漂移現(xiàn)象(CRUD-induced power shift, CIPS)[2]。

為了評(píng)估CRUD內(nèi)的傳熱傳質(zhì)現(xiàn)象，許多學(xué)者開發(fā)了一系列軟件，對(duì)在壓水堆運(yùn)行條件下CRUD內(nèi)的傳熱傳質(zhì)進(jìn)行了模擬研究[3-11]。例如，EPRI的BOA(boron-induced offset anomaly)[12]和CASL的MAMBA-3D[13]。但關(guān)于CRUD內(nèi)煙囪沸騰和形貌結(jié)構(gòu)影響的研究仍相對(duì)較少。為此本文建立一套二維的傳熱傳質(zhì)預(yù)測(cè)方法，實(shí)現(xiàn)CRUD內(nèi)的傳熱、壓降、流速、溶質(zhì)擴(kuò)散和化學(xué)反應(yīng)等多物理現(xiàn)象的耦合，并獲得冷卻劑在煙囪表面蒸發(fā)對(duì)傳熱和硼吸附的影響，以及傳熱傳質(zhì)隨沉積層孔隙率、厚度、煙囪半徑、煙囪直徑等形貌參數(shù)的變化規(guī)律。

1 建模方法

圖1 CRUD內(nèi)傳熱、流動(dòng)、傳質(zhì)和化學(xué)反應(yīng)等多物理現(xiàn)象(a)和計(jì)算區(qū)域及邊界(b)示意圖Fig.1 Schematics of multi-physical phenomena within CRUD depositions (a) and simulation domain and boundary (b)

各種溶質(zhì)之間保持化學(xué)平衡。為區(qū)分CRUD內(nèi)蒸汽和水的不同位置，將煙囪區(qū)域稱為干區(qū)，而環(huán)繞它周圍的圓柱殼體稱為濕區(qū)[8,15]。煙囪表面的沸騰增強(qiáng)了溶質(zhì)的濃縮，不斷濃縮的溶質(zhì)反過來影響著溶質(zhì)間的化學(xué)平衡和飽和溫度。選擇1個(gè)濕區(qū)作為計(jì)算區(qū)域，以沉積層厚度方向z(燃料棒徑向方向)和垂直于z的r方向(燃料棒軸向方向)建立二維圓柱坐標(biāo)，如圖1b所示。

1.1 傳熱模塊

CRUD內(nèi)部的傳熱可分為多孔介質(zhì)的導(dǎo)熱、燃料包殼的傳熱和煙囪邊界沸騰帶走的熱量3個(gè)部分，由能量守恒得到的控制方程及邊界條件如式(1)～(5)所示。

(1)

(2)

(3)

(4)

(5)

式中：kCRUD為CRUD濕區(qū)內(nèi)固體和液體加權(quán)后的總導(dǎo)熱系數(shù)，受孔隙率和溫度等參數(shù)的影響；qclad為包殼與CRUD交界面的熱流密度；d為CRUD厚度；hc為冷卻劑與CRUD表面的對(duì)流換熱系數(shù)，在計(jì)算中取定值12 000 W/(m·K)[10]；Tf為冷卻劑主流區(qū)溫度；rc和Rc分別為煙囪半徑和煙囪中心到濕區(qū)對(duì)稱邊界的距離；he為煙囪表面水的蒸發(fā)換熱系數(shù)，由Pan[9]提出，用于計(jì)算煙囪表面發(fā)生SNB帶走的能量，如式(6)所示；Tsat為飽和溫度，受溶質(zhì)濃度等參數(shù)的影響。

(6)

1.2 毛細(xì)流動(dòng)模塊

CRUD毛細(xì)熱管內(nèi)的流速較小，因此內(nèi)部的流動(dòng)符合達(dá)西定律。假設(shè)CRUD內(nèi)冷卻劑是不可壓縮流體，控制方程及邊界條件如式(7)～(11)所示。

(7)

(8)

p|z=d=pf

(9)

(10)

(11)

式中：ε為CRUD孔隙率；ρw為水的密度；κ為CRUD滲透率；μw為冷卻劑動(dòng)力黏度；pf為冷卻劑主流區(qū)壓力。

根據(jù)達(dá)西定律，水在多孔介質(zhì)中的速度與壓力梯度呈正比，如式(12)所示。

(12)

1.3 溶質(zhì)擴(kuò)散與化學(xué)反應(yīng)模塊

冷卻劑在毛細(xì)作用下以u(píng)l的速度進(jìn)入CRUD，同時(shí)將H3BO3和Li+等溶質(zhì)帶入CRUD中。在冷卻劑流動(dòng)帶來的對(duì)流、化學(xué)勢(shì)驅(qū)動(dòng)下的擴(kuò)散和電場(chǎng)對(duì)帶電粒子作用下的遷移三者的共同作用下，溶質(zhì)不斷被濃縮，同時(shí)維持著溶質(zhì)之間的化學(xué)平衡[11,16]。各溶質(zhì)的濃度控制方程和邊界條件如式(13)～(17)所示。

(13)

(14)

C|z=d=Cf

(15)

(16)

(17)

式中：控制方程右側(cè)第1項(xiàng)代表化學(xué)反應(yīng)引起的濃度變化，第2項(xiàng)代表輸運(yùn)引起的濃度變化；對(duì)于穩(wěn)態(tài)情況下，式(13)左右側(cè)均等于0；J為溶質(zhì)摩爾流量；Cf為冷卻劑主流區(qū)溶質(zhì)濃度；nR為徑向方向單位向量。

輸運(yùn)過程由能斯特-普朗克方程控制，如式(18)所示。

(18)

式中：D為擴(kuò)散系數(shù)，需要基于分形理論對(duì)純水中的擴(kuò)散系數(shù)進(jìn)行修正[7]，如式(19)所示，其中τ為曲折度；z為溶質(zhì)電荷數(shù)；F為法拉第常數(shù)；Φ為電勢(shì)。

(19)

(20)

式中：kf和kr分別為正、逆反應(yīng)速率常數(shù)；對(duì)于穩(wěn)態(tài)，式(20)左右側(cè)等于0，此時(shí)反應(yīng)物和生成物滿足化學(xué)平衡，如式(21)所示。

(21)

式中：m為質(zhì)量摩爾濃度；γ為活度系數(shù)，受到溫度和濃度的影響；ζ為溶液體積修正項(xiàng)，是濃度的函數(shù)。因此，化學(xué)平衡狀態(tài)下的正逆反應(yīng)速率常數(shù)K可以表示為式(21)。

本文所考慮的化學(xué)反應(yīng)主要包括水的電離、硼酸的形式轉(zhuǎn)化和沉淀的析出和溶解。

水的電離：

(22)

硼酸的形式轉(zhuǎn)化：

(23)

(24)

(25)

沉淀的析出和溶解：

4H3BO3+2Li+

(26)

1.4 數(shù)值計(jì)算方法

采用有限體積法進(jìn)行柱坐標(biāo)下的離散，使用溫度、壓力和各溶質(zhì)濃度中的最大相對(duì)誤差作為收斂條件，允許的最大相對(duì)誤差小于10-7。上述三模塊之間互相影響，在耦合計(jì)算過程中互相反饋直至最終計(jì)算結(jié)果達(dá)到收斂。CRUD傳熱傳質(zhì)計(jì)算流程圖如圖2所示。傳熱模塊獲得溫度分布，并為毛細(xì)流動(dòng)模塊和溶質(zhì)擴(kuò)散與化學(xué)反應(yīng)模塊更新物性。毛細(xì)流動(dòng)模塊得到的壓力和流速分別用于更新蒸汽物性和確定溶質(zhì)傳輸?shù)牧魉?。溶質(zhì)擴(kuò)散與化學(xué)反應(yīng)模塊得到的濃度更新飽和溫度進(jìn)而影響溫度和壓力分布。

圖2 CRUD傳熱傳質(zhì)計(jì)算流程圖Fig.2 Flow chart of simulation process of CRUD heat and mass transfer

2 對(duì)比驗(yàn)證

Haq等[10]在Cohen多孔介質(zhì)模型和Pan模型的基礎(chǔ)上，實(shí)現(xiàn)了CRUD內(nèi)的溫度場(chǎng)、流動(dòng)場(chǎng)和濃度場(chǎng)的耦合計(jì)算。使用Haq等相同的數(shù)據(jù)(表1)，獲得的結(jié)果與Haq等的結(jié)果進(jìn)行比較，如圖3、4所示。對(duì)比顯示，本文獲得的結(jié)果在趨勢(shì)上與Haq等的計(jì)算結(jié)果一致，但在數(shù)值上略有差別，最大相對(duì)誤差為0.83%。這主要是由于導(dǎo)熱系數(shù)上的差別，算例基于分形理論確定了CRUD的導(dǎo)熱系數(shù)，而Haq等在算例中的CRUD導(dǎo)熱系數(shù)采用了定值0.506 W/(m·K)。

表1 對(duì)比驗(yàn)證計(jì)算參數(shù)Table 1 Value of parameter used in validation analysis

圖3 多孔殼體中心軸向方向溫度變化Fig.3 Axial temperature in shell center

圖4 CRUD與包殼表面徑向方向溫度變化Fig.4 Radial temperature at CRUD and cladding interface

計(jì)算結(jié)果顯示，CRUD內(nèi)溫度最大位置出現(xiàn)在多孔沉積層中心靠近包殼的區(qū)域，略高于煙囪表面區(qū)域的溫度，這主要是由于毛細(xì)流在煙囪表面發(fā)生沸騰帶走了部分熱量。整個(gè)CRUD區(qū)域內(nèi)，煙囪表面靠近包殼區(qū)域的壓力最小，但硼酸和鋰濃度都在此區(qū)域內(nèi)達(dá)到最大值。

3 結(jié)果討論

為評(píng)估CRUD形貌參數(shù)對(duì)CRUD的傳熱傳質(zhì)的影響，使用多組CRUD形貌參數(shù)和典型壓水堆熱工參數(shù)進(jìn)行了計(jì)算，如表2所列。

表2 輸入?yún)?shù)Table 2 Value of input parameter

采用表中標(biāo)準(zhǔn)輸入?yún)?shù)得到的溫度、壓力、硼酸濃度分布如圖5所示。CRUD的溫度在與包殼交界面處達(dá)到最大值(634.60 K)，與相同r位置CRUD與冷卻劑交界面的溫差達(dá)到18.74 K。r方向上溫差均較小，最大值僅為1.20 K。最大SNB熱流密度發(fā)生在煙囪表面靠近包殼區(qū)域內(nèi)，達(dá)到0.31 MW/m2。同時(shí)這一區(qū)域也是最大壓降的位置，達(dá)到1.8 kPa，使得r方向流速達(dá)到0.52 cm/s。將總能量分為通過煙囪表面SNB帶走的能量和通過CRUD傳導(dǎo)至CRUD與冷卻劑表面通過冷卻劑對(duì)流換熱帶走的能量?jī)刹糠?，將兩者占總能量的比例分別稱為SNB占比和熱傳導(dǎo)占比，在分析單元內(nèi)，SNB占比達(dá)到75.07%。

圖5 溫度(a)、壓力(b)和硼酸濃度(c)分布Fig.5 Distributions of temperature (a), pressure (b), and boric acid concentration (c)

3.1 厚度的影響

圖6 硼質(zhì)量和SNB占比隨厚度的變化Fig.6 Boron mass and SNB ratio vs. CRUD thickness

壓水堆包殼表面發(fā)現(xiàn)的CRUD厚度變化范圍為20～90 μm，與壓水堆所處的熱工和水化學(xué)條件、運(yùn)行時(shí)間和冷卻劑雜質(zhì)去除情況有關(guān)。隨CRUD厚度變化的單位面積內(nèi)硼質(zhì)量和SNB占比如圖6所示。隨著CRUD厚度的增加，硼也不斷得到積累，質(zhì)量不斷增加，在厚度達(dá)到40 μm時(shí)開始析出Li2B4O7，這與EPRI[2]報(bào)道的硼析出大致發(fā)生在35～42 μm相符合。在發(fā)生硼析出后，可溶性硼的增加趨勢(shì)變緩，更多硼以沉淀的形式析出。同時(shí)，厚度達(dá)到40 μm后SNB占比有微弱增大。

3.2 孔隙率的影響

針對(duì)壓水堆內(nèi)CRUD的表征分析表明CRUD的孔隙率約為0.4～0.8，并且會(huì)隨著硼和金屬離子的沉淀反應(yīng)、金屬氧化物和氫氣的氧化還原反應(yīng)等變化，但在現(xiàn)有模擬中一般不考慮這些反應(yīng)對(duì)孔隙率的影響。隨孔隙率變化的最大溫差、SNB占比和最大硼酸濃度如圖7所示。由于水的導(dǎo)熱系數(shù)小于構(gòu)成CRUD多孔骨架的金屬氧化物，因此孔隙率越高，導(dǎo)熱系數(shù)越低，溫差越大。SNB占總能量的比例與孔隙率呈近似線型比例關(guān)系，這可能是由于孔隙率越大，毛細(xì)熱管內(nèi)的冷卻劑越多，相應(yīng)地通過SNB帶走的熱量也更多。采用分形理論修正后的擴(kuò)散系數(shù)隨孔隙率變大而變大，在冷卻劑毛細(xì)流動(dòng)流速較小且變化不大的情況下，更大的擴(kuò)散系數(shù)僅需較小的擴(kuò)散梯度即可與硼酸的對(duì)流建立平衡，因此孔隙率越高硼酸濃度越低。

3.3 煙囪內(nèi)徑的影響

沉積層內(nèi)的煙囪對(duì)于CRUD傳熱具有極其重要的意義，因?yàn)榻^大部分的蒸發(fā)相變都發(fā)生在煙囪表面。研究表明，CRUD煙囪的直徑范圍一般為2～5 μm。隨煙囪內(nèi)徑變化的最大溫差、SNB占比、最大壓降和最大硼酸濃度如圖8所示。對(duì)于固定的煙囪密度，煙囪內(nèi)徑的增大意味著“濕區(qū)”所占面積份額的減小和煙囪表面相變面積的增大，因此相應(yīng)的SNB占比也相應(yīng)增大，熱傳導(dǎo)占比減小，最大溫差也減小。SNB的增加，也引起了最大壓降的增大，毛細(xì)流動(dòng)流速增大，這最終引起了最大硼酸濃度的增大。

圖7 最大溫差(a)、SNB占比(b)和最大硼酸濃度(c)隨孔隙率的變化Fig.7 Maximum temperature difference (a), SNB ratio (b), and boric acid concentration maximum (c) vs. porosity

圖8 最大溫差(a)、SNB占比(b)、最大壓降(c)和最大硼酸濃度(d)隨煙囪內(nèi)徑的變化Fig.8 Maximum temperature difference (a), SNB ratio (b), maximum pressure drop (c), and boric acid concentration maximum (d) vs. chimney diameter

圖9 最大溫差(a)、SNB占比(b)、最大壓降(c)和最大硼酸濃度(d)隨煙囪密度的變化Fig.9 Maximum temperature difference (a), SNB ratio (b), pressure drop maximum (c), and boric acid concentration maximum (d) vs. chimney density

3.4 煙囪密度的影響

CRUD內(nèi)煙囪密度的變化范圍較大，在不同反應(yīng)堆、堆芯不同位置甚至同一燃料棒不同方位上都有所不同。一般來說，煙囪的密度范圍為330～4 500 mm-2。隨煙囪密度變化的最大溫差、SNB占比、最大壓降和最大硼酸濃度如圖9所示。隨著單位面積內(nèi)煙囪數(shù)量的增大，CRUD的最大溫差減小，SNB占比、最大壓降、最大硼酸濃度都相應(yīng)增大，這一趨勢(shì)也與煙囪內(nèi)徑增大所引起的趨勢(shì)相同，均與單位面積內(nèi)SNB增強(qiáng)有關(guān)。

4 結(jié)論

本研究建立了一種用于CRUD內(nèi)傳熱傳質(zhì)現(xiàn)象的預(yù)測(cè)方法，實(shí)現(xiàn)了CRUD內(nèi)傳熱、流動(dòng)、溶質(zhì)輸運(yùn)和化學(xué)反應(yīng)等多物理現(xiàn)象的耦合計(jì)算。此方法合理地預(yù)測(cè)了CRUD內(nèi)的溫度、壓力、流速和濃度分布。

在所分析的多物理現(xiàn)象中，在CRUD獨(dú)特結(jié)構(gòu)煙囪表面發(fā)生的沸騰相變對(duì)傳熱傳質(zhì)均具有重要影響，SNB越強(qiáng)，CRUD的最大溫差越小，最大壓降越大，硼酸濃度也越大。通過參數(shù)敏感性分析，明晰了隨厚度的增大，硼的存在形式從可溶性硼向Li2B4O7轉(zhuǎn)化的趨勢(shì)；解明了影響CRUD內(nèi)傳熱傳質(zhì)的孔隙率、煙囪內(nèi)徑和煙囪密度對(duì)SNB均具有正相關(guān)性，而對(duì)硼濃度分別具有負(fù)、正、正相關(guān)性。該研究有利于進(jìn)一步理解CRUD內(nèi)硼吸附現(xiàn)象，可為堆芯功率漂移現(xiàn)象CIPS的預(yù)測(cè)提供技術(shù)支持。