基于灰色關聯(lián)度的棒束通道內液態(tài)鉛鉍合金中顆粒物運動沉積研究

周 濤 唐劍宇 蔣 屹

1（東南大學 能源與環(huán)境學院核科學與技術系 南京 210096）

2（核熱工安全與標準化研究團隊 南京 210096）

3（大型發(fā)電裝備安全運行與智能測控國家工程研究中心 南京 210096）

4（華北電力大學 核科學與工程學院 北京 102206）

鉛鉍快堆作為國際第四代核能系統(tǒng)的優(yōu)選技術路線之一，具有中子經濟性好、嬗變率高、功率密度高等優(yōu)點；同時，作為冷卻劑的鉛鉍共晶具有傳熱性能強、抗輻照能力好的特點，這使得鉛鉍堆的系統(tǒng)組成和配套設施得以簡化，相比其他堆型更容易實現小型模塊化，成為構建小堆的良好選擇。雖然采用鉛鉍快堆構建小堆具有造價更低、更易建造、更安全靈活的優(yōu)點，但液態(tài)鉛鉍合金在其流動過程中會對其流經的冷卻劑通道產生腐蝕沖刷、腐蝕，以及管道流體加速腐蝕效應，產生的不溶性顆粒物將在冷卻劑通道、棒束通道內局部聚集，其沉積效應可能會導致堆芯過熱，嚴重時會導致堆芯熔化［1］。閆靜賢等［2］利用CFD方法從三維的對流-傳質擴散方程出發(fā)，研究了鉛鉍合金流體回路中，腐蝕所產生的顆粒物的沉積情況，并建立了腐蝕模型。王琛等［3］通過實驗和模擬相結合的方式，討論了液態(tài)鉛鉍流體在管內流動時，不同湍流普朗特數模型和不同湍流模型的選取對計算結果的影響，提出了不同普朗特數下推薦采用的物理模型。周濤等［4-5］對鉛鉍流體熱工力學及安全特性以及核電中顆粒物運動沉積方法進行了整理，并提出了未來進一步的研究方向。陳娟等［6］對液態(tài)鉛鉍合金中的金屬納米顆粒的熱泳現象進行研究，分析了不同種類的納米顆粒在液態(tài)鉛鉍中的熱泳速度。鄒文重等［7］對注入氣體后液態(tài)鉛鉍流動換熱特性進行了研究，得到了注入氣體后不同條件下液態(tài)鉛鉍的溫度與壓力分布。阮章順等［8］對鉛鉍快堆燃料包殼材料進行了研究，分析了鉛鉍合金的溫度與氧含量對包殼材料的腐蝕行為。劉鵬等［9］對鉛鉍合金系統(tǒng)的膨脹箱內的固態(tài)顆粒物運動進行研究，得出鉛鉍合金中氧化鋰顆粒物的運動規(guī)律。目前研究工作［10-13］主要集中于液態(tài)鉛鉍合金的流動換熱特性研究及液態(tài)鉛鉍合金與結構材料的相容性研究。而對于液態(tài)鉛鉍合金在小型快堆堆芯棒束通道內的顆粒物沉積運動研究不足。通過對不溶性粒狀腐蝕產物在堆芯燃料棒流域的沉積過程進行計算研究，得到小型鉛鉍快堆棒束通道內顆粒物沉積運動情況，以便了解其對小型鉛鉍快堆安全運行的影響，為探究鉛鉍反應堆安全設計提供參考。

1 研究對象

1.1 幾何模型

依據陳釗［14］的小型模塊化自然循環(huán)鉛鉍快堆SNCLFR-100的堆芯設計方案，其燃料組件通道分布如圖1所示。

圖1 燃料組件模型Fig.1 Fuel assembly model

從圖1可以看出，堆芯燃料組件內部的通道按相對位置和壁面條件分為三種類型：類三角形通道、類五邊形通道和類梯形通道。通道長度為400mm，相鄰棒心間距為7mm，燃料元件外徑為6mm。

1.2 網格劃分

使用ANSYS ICEM CFD 前處理軟件對通道模型進行網格劃分，計算區(qū)域為單純的流體域。為保證計算精度，采用六面體結構網格。為適應k-ε湍流模型，對邊界網格進行加密處理，保證y+值位于30～300，整體網格整體質量0.7 以上。通道網格情況如圖2所示。

圖2 內三角(a)、類五邊形(b)和類梯形(c)通道網格劃分Fig.2Mesh model of triangular channel (a), pentagonal channel (b) and trapezoidal channel (c)

1.3 參數設置

液相選擇鉛鉍合金流體，鉛鉍合金流體的物性參數選擇參考koji 等［15］、蘇子威等［16］的研究結果。在實際運行過程中，管道基體材料會氧化脫落［11］且在液態(tài)鉛鉍合金（Lead-Bismuth Eutectic，LBE）作為冷卻劑材料時常會采用LiH與O2反應來降低氧氣含量［9］。因此，會產生不同種類的顆粒物，綜合考慮后選擇Li2O、PbO、Fe3O4作為研究對象。由于LBE 中固態(tài)顆粒的溶解與沉積始終處于一種動態(tài)平衡過程中，整體顆粒物含量較少且處于稀相，因此，選擇進口顆粒物濃度為1%。相關參數設置如表1所示。

表1 參數設置Table 1 Parameter settings

在表1中，T指流體當前的溫度，K；鉛鉍流體的密度、定壓比熱容、熱導率、動力黏度均是與溫度相關的函數。流體湍流模型選取標準k-ε模型，近壁面處理采用標準壁面函數，求解方式采用SIMPLE 算法，其中壓力項采用標準差分方式，湍動能、動量和能量等均采用二階迎風差分格式，這種格式使得計算的速度較快，同時又具有較好的精度和收斂性；離散相與流體的耦合計算采用雙向耦合方式，以模擬顆粒物與流場之間的相互影響，步長采用默認步長。進口邊界條件為速度進口，并作為顆粒物入射面，出口邊界條件為壓力出口，并作為顆粒物出口面；燃料棒表面邊界條件為定熱流密度，組件壁面可認為與相鄰壁面溫度相同，不發(fā)生傳熱，設置為絕熱壁面，對于顆粒物，壁面為吸收面；同時編寫UDF 定義顆粒物于壁面的邊界條件，根據顆粒物的臨界速度判斷其是trap 還是reflect。通道中的對稱面設置為對稱邊界條件，顆粒物濃度采用體積占比來表示。

1.4 網格無關性檢驗

1） 流體網格無關性驗證

為了保證計算網格的獨立性以及時間步長的收斂性，根據§1.3的參數設置計算不同網格數量條件下三種管道的出口速度，得到的結果如圖3所示。

圖3 三種通道的出口速度Fig.3 Outlet velocity of different channels

由圖3可以得到，當網格數量大于450000時，內通道的流體出口速度幾乎保持不變，當網格數量大于500000時，角通道的流體出口速度幾乎保持不變，當網格數量大于550000時，邊通道的流體出口速度幾乎保持不變。

2） 顆粒物網格無關性驗證

根據§1.3 參數設置，選取兩種不同粒徑的顆粒物，計算在不同網格條件下管道出口濃度，結果如圖4所示。

圖4 不同粒徑三種通道的出口顆粒物濃度 (a) 內通道，(b) 角通道，(c) 邊通道Fig.4 Outlet particle concentration of different channels with different particle sizes(a) Triangular channel, (b) Pentagonal channel, (c) Trapezoidal channel

由圖4可以得到，對于內通道，網格數量大于450000時，不同粒徑下出口顆粒物濃度隨網格變化較小。對于角通道，網格數量大于450000時，不同粒徑下出口顆粒物濃度隨網格變化較小。對于邊通道，網格數量大于550000時不同粒徑下出口顆粒物濃度隨網格變化較小。綜合兩次驗證結果，考慮計算機運行效率，確定內通道模型網格節(jié)點數為495 662，邊通道模型網格節(jié)點數為587 894，角通道模型網格節(jié)點數為506 543。

2 計算模型

2.1 標準k-ε模型

標準k-ε模型屬于二方程模型，適用于完全湍流狀態(tài)且分子黏性可以忽略的流動模擬，是目前工程流場的計算中應用最多、適用范圍最廣的湍流模型。因此，模擬選用該模型進行計算。標準k-ε模型是以湍動能k方程為基礎，通過引入湍動能耗散率ε而建立的輸運方程。該模型中的ε定義為：

湍流黏度μt可表示成k與ε的函數：

標準k-ε模型的輸運方程為：

式中：μ為流體動力黏度，Pa·s；ρ為流體 密 度，kg·m-3；ui為時均速度，m·s-1；Gk為層流速度梯度產生的湍流動能，kg?m2·s-2；Gb為浮力產生的湍流動能，kg?m2·s-2；YM為可壓縮湍流中脈動擴張項，kg?m2·s-2；Cμ、C1ε、C2ε和C3ε為經驗常量；σk和σε分別為k方程和ε方程的湍流普朗特數；Sk和Sε為k方程和ε方程的用戶自定義源項，kg?m2·s-2；xi、xj、xk是方向坐標；x、y、z分別記作x1、x2、x3，其中i=1、2、3；j=1、2、3；k=1、2、3。

2.2 離散相計算模型

在離散相（Direct Part Marking，DPM）中，流場中顆粒的運動軌跡是通過對拉氏坐標系下的顆粒作用力微分方程進行積分得到的，顆粒的作用力平衡方程可表示為：

式中：FD(u-up)是單位顆粒質量所受到的流體曳力，其中：

式中：u為流體速度，m·s-1；up為顆粒速度，m·s-1；μ為流體動力黏度，Pa·s；ρ為流體密度，kg·m-3；ρp為顆粒密度，kg·m-3；dp為顆粒粒徑，m；CD為曳力系數；Rep為相對雷諾數；α為彎管角度；t為時間，s；gx為重力加速度，m·s-2；Fx為其他作用力，N·kg-1，包括浮力、壓力梯度力、附加質量力、沙夫曼力、熱泳力、湍流力等。

2.3 湍流普朗特數

因為液態(tài)鉛鉍的特殊性質，其湍流普朗特數較小，常規(guī)模型不再適合。目前，針對液態(tài)鉛鉍湍流普朗特數的計算模型有多種，考慮到定壁面熱流密度條件下，Cheng［17］提出的湍流普朗特數模型的數值模擬結果與實驗數據的吻合度更好。其模型表達為：

其中：

式中：Prt表示湍流普朗特數；Re表示雷諾數；Pr是分子普朗特數。

2.4 灰色關聯(lián)度

灰色關聯(lián)分析是將因素之間發(fā)展趨勢的相似或相異程度，作為衡量因素間關聯(lián)程度的一種方法，并通過建立灰色關聯(lián)分析模型，來衡量因素間的關聯(lián)程度大小。

X0為參考序列，表示系統(tǒng)中具有行為特征的數據序列；Xi為比較序列，表示影響系統(tǒng)行為的數據序列。為保證數據處理準確，兩種數據序列據需要進行無量綱化處理。這兩種數據［18］均為1-時距序列。比較序列在j時刻上的觀測數據為Xi（j），可表示為：

X0和Xi在區(qū)間［k-1，k］，（k=2，3，…，n）上的斜率所組成的序列號，可表示為：X0和Xi在各個對應時段的斜率比值所組成的序列表達式可寫為：

序列K0=(k01，k02，…，k0(n-1))的變異系數表達式可寫為：

其中：

序列K0i=(k01ki1，k02ki2，…，k0(n-1)ki(n-1))的廣義變異系數表達式可寫為：

其中：

綜合以上公式，可得到比較序列Xi、參考序列X0的灰色關聯(lián)度結果，可寫為：

計算序列之間的灰色關聯(lián)度，應將序列與序列之間的結果當作整體來看待，避免單獨的點與點之間的數據對比?？紤]到序列之間存在正負相關性，可通過加權平均取絕對值的方式抵消引起的正負行問題。并且最終得到的結果位于［-1，1］之間，根據曹明霞［19］所得模型，X0和Xi在每一段的斜率比值接近1，灰色關聯(lián)度就越大。

3 數值模擬結果及分析

3.1 各通道速度分布

選取流體進口速度2m·s-1，進口溫度533K，顆粒物種類為Fe3O4，加熱段熱流密度為80000 W·m-2的工況，各通道軸向速度分布如圖5所示，各通道充分湍流段速度分布云圖如圖6所示。

圖5 通道軸向速度分布Fig.5 Axial velocity distribution of channels

圖6 通道速度分布Fig.6Velocity distribution in turbulent section of channels

從圖5和圖6可以看出，由于流場受到堆芯流域結構的擠壓，通道中徑向速度的分布整體上沿著中心至壁面梯度呈下降趨勢。沿流體流動的方向，通道內流體中心速度逐步增大后趨于穩(wěn)定，邊界速度在進口處快速下降后上升，隨后趨于穩(wěn)定。其原因在于：進口處邊界處流體受壁面加熱以及自身黏滯力的影響，橫向流動較弱，而隨著軸向距離的增加，壁面黏滯力的影響逐步向流體內部傳遞，邊界層增厚轉變?yōu)橥牧鬟吔鐚?，流體的脈動強化了動量傳遞，使得速度的變化變得平緩。

3.2 通道湍動能分布

選取流體進口速度2m·s-1，進口溫度533K，顆粒物種類為Fe3O4，加熱段熱流密度為80000 W·m-2的工況，各通道軸向1mm、200mm、399mm處的鉛鉍合金湍動能分布如圖7所示。

圖7 內三角形通道(a)、類五邊形通道(b)、類梯形通道(c)湍動能分布Fig.7 Triangular channel (a), pentagonal channel (b), trapezoidal channel (c) turbulent kinetic energy distribution

從圖7可以看出，類三角形通道的湍動能呈中心對稱分布，谷值出現在流域和對稱面中部，類五邊形通道和類梯形通道流域中部的湍流速度漲落較小，從而出現湍動能谷值。從圖7中還可以看出，湍動能的峰值出現在燃料棒外表面近壁面區(qū)域，產生該現象的原因是湍流邊界層內的速度呈現的脈動特性使得近壁面流體的速度漲落較大，導致湍流強度增強。

3.3 顆粒物濃度分布

3.3.1 軸向顆粒物分布

選取流體進口速度2m·s-1，進口溫度533K，顆粒物種類為Fe3O4，加熱段熱流密度為80000 W·m-2的工況，探測點為距離加熱壁面0.1mm 的位置處，各通道內顆粒物軸向濃度分布如圖8所示。

圖8 軸向顆粒物濃度變化Fig.8Variation of particle concentration distribution in streamline direction

從圖8可以看出，三種通道內顆粒物的變化趨勢趨于一致，隨著軸向距離的加大，管道中顆粒物的濃度在入口處先迅速下降，并在管道中后部逐步趨于穩(wěn)定。表明通道內顆粒物的沉積主要發(fā)生在入口階段。原因在于湍流分為湍流核心區(qū)和湍流邊界層，湍流邊界層又可分為過渡區(qū)和黏性底層兩個流動區(qū)，顆粒物的沉積需要經過三個流動區(qū)。在湍流核心區(qū)，顆粒物容易向壁面運動，并進入到過渡區(qū)，而過渡區(qū)的顆粒物需要具有足夠慣性，才能穿過過渡區(qū)進入黏性底層并在該區(qū)域形成運動與沉積的動態(tài)平衡。此時顆粒物是否沉積與其受力大小相關。在入口處，近壁面流體速度梯度大，越靠近壁面流體的速度越低，顆粒物的流場跟隨性也隨之減弱，使得顆粒物在入口處大量沉積。而隨著流體湍流強度逐步增強，附壁層減薄，近壁面湍動能逐步增大使得顆粒物平均速度加大、不規(guī)則運動加劇，受湍流擴散作用克服慣性力而向著中心富集。

3.3.2 徑向顆粒物分布

選取流體進口速度2m·s-1，進口溫度533K，顆粒物種類為Fe3O4，加熱段熱流密度為80000 W·m-2的工況，探測點為通道軸向距離240mm 位置處，各通道內顆粒物徑向濃度分布如圖9所示。

圖9 徑向顆粒物濃度變化 (a) 類三角形通道，(b) 類五邊形通道，(c) 類梯形通道Fig.9Variation of particle concentration distribution in radius location(a) Triangular channel, (b) Pentagonal channel, (c) Trapezoidal channel

從圖9可以看出，在湍流充分發(fā)展段，湍動能的大小是影響顆粒物濃度分布的主要因素，湍動能的分布與顆粒物濃度分布成反比關系。顆粒物濃度的谷值出現在燃料棒的近壁面附近，這是由于燃料棒的近壁面附近流體的湍動能漲落較大，不利于顆粒物沉積。

3.3.3 壁面顆粒物分布

選取流體進口速度2m·s-1，進口溫度533K，加熱段熱流密度為80000 W·m-2，顆粒物種類為Fe3O4的工況，各通道壁面顆粒物濃度分布如圖10所示。

圖10 通道壁面顆粒物濃度分布 (a) 類三角形通道，(b) 類五邊形通道，(c) 類梯形通道Fig.10 Particle concentration distribution in channel wall(a) Triangular channel, (b) Pentagonal channel, (c) Trapezoidal channel

從圖10可以看出，隨著管道長度的增加，各通道的顆粒物沉積規(guī)律類似，表現為在進口段燃料棒表面發(fā)生大面積附著沉積現象，中段及后段以點狀沉積現象為主。在包殼表面點狀沉積區(qū)域會出現包殼點蝕現象，點蝕區(qū)域也會引起局部傳熱惡化，破壞包殼完整性；大面積附著沉積會改變包殼材料的熱導率和堆芯中子通量分布，引起堆芯軸向功率偏移。

3.3.4 顆粒物種類對沉積率的影響

選取流體進口速度2m·s-1，進口溫度533K，加熱段熱流密度為80000 W·m-2的工況，顆粒物種類選取PbO、Li2O、Fe3O4，粒徑2μm。通過模擬計算，不同通道內顆粒物沉積率如圖11所示。

圖11 顆粒物沉積率隨顆粒物種類的變化關系Fig.11 Relationship between particle deposition rate and particle type

從圖11可以看出，不同通道內沉積率變化趨勢相同，顆粒物沉積率隨著顆粒物的密度而增大，原因在于較大密度顆粒物的單位質量較大，慣性較大，造成流場跟隨性弱。從圖11還可以看出，顆粒物的沉積率隨著通道當量直徑的減小而增大。原因在于管道的當量直徑越大，管道主流流域的速度分布較均勻，湍流度較小，顆粒物在流體中的跟隨性較好，沉積率更低。這與朱亮宇等［20］得到的顆粒物沉積率隨管道直徑的變化趨勢一致。

3.3.5 顆粒物粒徑對沉積率的影響

選取流體進口速度2m·s-1，進口溫度533K，加熱段熱流密度為80000 W·m-2的工況，顆粒物種類選取PbO，粒徑選取0.5μm、1μm、5μm 和10μm。通過模擬計算，不同通道內顆粒物沉積率如圖12所示。

圖12 顆粒物沉積率隨顆粒物粒徑的變化關系Fig.12 Relationship between particle deposition rate and particle size

從圖12可以看出，不同通道內沉積率變化趨勢相同，顆粒物沉積率隨著顆粒物的粒徑增大而增大。原因在于流速一定的情況下，粒徑越大的顆粒物會導致單顆粒的質量增大，流場跟隨性變弱，越容易在進口處發(fā)生沉積。

3.3.6 進口速度對沉積率的影響

選取進口溫度533K，加熱段熱流密度為80000 W·m-2，顆粒物種類為Fe3O4，粒徑2μm 的工況，分別設置流體進口速度為0.8m·s-1、2m·s-1、4m·s-1，通過模擬計算，不同通道內顆粒物沉積率如圖13所示。

圖13 顆粒物沉積率隨顆粒物速度的變化關系Fig.13 Relationship between particle deposition rate and particle velocity

從圖13可以看出，不同通道內沉積率變化趨勢相同，顆粒物沉積率隨著流體流速的增大而降低。原因在于流速增大使得顆粒物受到的切應力增大，導致原來低速條件下沉積在壁面的顆粒物脫離管壁，降低了壁面附近的顆粒物濃度。顆粒物在通道內的沉積容易導致傳熱惡化，上述結果表明：增大鉛鉍的流速，有利于回路的正常安全運行。

3.3.7 各參數對顆粒物沉積的灰色關聯(lián)度分析

利用灰色關聯(lián)度分析在其余條件相同時顆粒物種類、顆粒物濃度、進口速度對顆粒物沉積的影響程度。沉積率為管道中沉積的顆粒物與所有顆粒物之比，以沉積率作為目標參數，顆粒物種類、顆粒物濃度、進口速度作為影響參數。關聯(lián)度分辨系數分別選擇0.5、0.3、0.1進行計算，結果如圖14所示。

圖14 影響顆粒物沉積率的各因素關聯(lián)度Fig.14 Correlation degree of various factors affecting particle deposition rate

從圖14可以看出，選取的分辨率與各因素之間的關聯(lián)度大小呈負相關，同時關聯(lián)度的大小也受分辨率的影響。其中，粒徑跟沉積率的關聯(lián)程度最大，即兩者關系最為密切；其次，顆粒物種類跟沉積率的關聯(lián)程度較??；最后，速度跟沉積率的關聯(lián)程度最小。這是因為隨著顆粒物粒徑的增大，受重力的影響開始增大。而對于較大粒徑的顆粒物而言，重力是最主要的沉積影響因素，因而關聯(lián)程度最大。顆粒物種類主要受顆粒物自身密度的影響，即相同粒徑下不同顆粒質量不同所受重力不同。但由于顆粒物尺寸較小，不同顆粒物間的質量差別不大，因而關聯(lián)度較小。流體速度會對顆粒物有一個攜帶作用，形成切應力，在速度相差不是十分大時，相比重力對顆粒物的沉降影響來說較小，因而關聯(lián)度最小。

4 結語

通過對鉛鉍快堆棒束通道中細顆粒物運動沉積的計算，得到鉛鉍快堆棒束通道內顆粒物的沉積運動情況。并分析了不同顆粒物種類、顆粒物粒徑、顆粒物速度對顆粒物沉積運動的影響，并利用灰色關聯(lián)度分析法分析比較各參數對顆粒物沉積的影響程度大小。

1）反應堆堆芯流道內顆粒物的沉積主要發(fā)生在入口階段，具體表現為進口段大面積附著沉積，中段及后段以點狀沉積為主；隨著軸向距離的加大，湍動能大小是影響顆粒物徑向分布的主要因素。

2）顆粒物密度的增大會加強子通道內顆粒物的沉積情況。顆粒物的密度增大使得顆粒物的單位質量增大，造成流場跟隨性變弱，更容易沉積。

3）顆粒物粒徑的增大會加強顆粒物的沉積情況。顆粒物粒徑的增大使得單顆粒的質量增大，流場跟隨性變弱，越容易在進口處發(fā)生沉積。

4）初始顆粒物進口速度的增大會降低顆粒物的沉積情況。顆粒物進口速度增大使得顆粒物受到的切應力增大，導致原來低速條件下沉積在壁面的顆粒物脫離管壁，降低了壁面附近的顆粒物濃度。增大鉛鉍的流速有利于回路的正常安全運行。

5）對顆粒物沉積的灰色關聯(lián)度大小為粒徑＞種類＞速度。其中，顆粒物速度對顆粒物沉積影響程度最小，顆粒粒徑影響最大。

作者貢獻聲明周濤負責提出設計和研究思路，實施研究方案，對文章的知識性內容作批評性審閱，獲取研究經費，指導；唐劍宇負責實施研究，采集數據，分析計算／解釋數據，文章部分起草和修改；蔣屹負責實施研究，采集數據，計算／解釋數據，文章部分起草。