壓力容器環(huán)形下降段內(nèi)含硼溶液流動混合特性的數(shù)值模擬分析

杜為安，張 琦，楊繼行，閆超星，張 翼，劉成洋

(1.中國艦船研究設(shè)計中心，湖北 武漢 430061；2.上海交通大學(xué) 核科學(xué)與工程學(xué)院，上海 200240)

典型壓水反應(yīng)堆在一回路冷卻劑內(nèi)添加硼酸來補償剩余反應(yīng)性，硼酸溶液在反應(yīng)堆壓力容器以及一回路管道內(nèi)的分布特性直接影響堆芯功率[1-3]。當(dāng)堆芯供冷不足、堆芯內(nèi)產(chǎn)生大量水蒸氣時，由于硼酸在蒸汽中溶解度較小，幾乎不含硼的水蒸氣經(jīng)過蒸汽發(fā)生器后將形成不含硼的冷凝液(稱為稀釋水團(tuán))。稀釋水團(tuán)將聚集在壓力容器入口上游的相關(guān)管路中[4]，這部分稀釋水團(tuán)進(jìn)入壓力容器后將會和含硼冷卻劑混合，并最終進(jìn)入堆芯。堆芯應(yīng)急冷卻系統(tǒng)啟動后，安全注射水箱中的高濃度硼溶液被注入壓力容器內(nèi)與一回路冷卻劑混合后進(jìn)入堆芯。安注冷卻劑和一回路冷卻劑之間還可能導(dǎo)致熱分層[5-6]。上述含硼溶液的混合過程直接關(guān)系到堆芯反應(yīng)性的變化，攪混不均勻可能導(dǎo)致堆芯反應(yīng)性激增或功率分布偏移，進(jìn)而威脅堆芯安全。

壓力容器的環(huán)形下降段位于堆芯入口上游，研究含硼溶液在環(huán)形下降段內(nèi)的混合過程對反應(yīng)堆安全具有重要意義[1,7]。國內(nèi)外已開展過含硼溶液混合相關(guān)的實驗[2,8-9]，受測量技術(shù)限制，現(xiàn)階段的相關(guān)實驗只能獲得環(huán)形下降段內(nèi)局部位置的平均硼濃度分布。計算流體力學(xué)(CFD)軟件成本較低、可重復(fù)性高，并可對三維場的流動與混合過程進(jìn)行分析，被廣泛應(yīng)用于反應(yīng)堆熱工水力的相關(guān)研究。本文基于計算流體力學(xué)軟件ANSYS Fluent 18.0對環(huán)形下降段內(nèi)的流動與混合特性進(jìn)行分析。

1 數(shù)值模型

1.1 控制方程

本文基于不可壓縮流體假設(shè)，忽略溫度變化對濃度擴(kuò)散的影響，控制方程如下。

1) 連續(xù)方程

連續(xù)方程即質(zhì)量守恒方程，其微分形式為：

(1)

式中：u、v、w為x、y、z方向速度分量，m/s；ρ為密度，kg/m3；t為時間，s。

2) 動量方程

動量方程實際上是牛頓第二定律，x、y、z方向的動量方程為：

(2)

(3)

(4)

式中：p為流體微元體上的壓力，Pa；τxx、τxy、τxz等為黏性應(yīng)力τ的分量，Pa；fx、fy、fz為x、y、z方向的單位質(zhì)量力，m/s2，若質(zhì)量力只有重力，z軸豎直向上，則fx=fy=0，fz=-g。

3) 能量方程

包含有熱交換的系統(tǒng)必須滿足能量守恒定律，其本質(zhì)為熱力學(xué)第一定律，表達(dá)式為：

(5)

式中：E為流體微團(tuán)總焓值，J/kg，由內(nèi)能、勢能和動能組成，E=h-p/ρ+u2/2；hj為組分j的焓值，J/kg；Jj為組分j的擴(kuò)散通量；k為流體熱傳導(dǎo)系數(shù)，W/(m·K)；Sh為包括化學(xué)反應(yīng)及其他體積熱源項。

4) 組分質(zhì)量守恒方程

在存在質(zhì)的交換或多種化學(xué)組分的特性控制系統(tǒng)中，每一組分都要遵循組分質(zhì)量守恒定律。對于一特定的系統(tǒng)，組分質(zhì)量守恒定律可描述為：系統(tǒng)內(nèi)某組分對時間的變化率等于通過系統(tǒng)界面的凈擴(kuò)散流量和化學(xué)反應(yīng)生成率之和。

由組分守恒定律，得到組分s的組分質(zhì)量守恒方程:

(6)

式中：cs、ρcs為組分s的體積濃度和質(zhì)量濃度；Ds為組分s的擴(kuò)散系數(shù)；Ss為源項。式(6)從左到右分別為時間變化率、對流項、擴(kuò)散項和源項。

1.2 模型建立與網(wǎng)格劃分

本文參照哈爾濱工程大學(xué)激光診斷實驗室的硼濃度混合分析實驗[10]建立相應(yīng)的物理模型，以便于對比實驗結(jié)果與數(shù)值模擬結(jié)果。表1為壓力容器物理模型的相關(guān)參數(shù)。圖1為壓力容器物性模型的俯視圖和正視圖。流體經(jīng)入口管進(jìn)入環(huán)形下降段，流入下腔室后從下腔室上部的7×7矩形通道流出。下腔室并不是本文關(guān)注的重點，在保證流體流動方向不變的情況下，簡化下腔室內(nèi)構(gòu)建，僅保留半球形結(jié)構(gòu)和上部孔板。

表1 壓力容器模型參數(shù)Table 1 Parameter of reactor model

圖1 壓力容器簡化模型Fig.1 Mock RPV model

采用ANSYS ICEM對模型進(jìn)行網(wǎng)格劃分，網(wǎng)格結(jié)構(gòu)如圖2所示。將壓力容器模型分為入口管段、環(huán)形下降段、下腔室和矩形出口通道4個主要區(qū)域，各區(qū)域均采用六面體結(jié)構(gòu)化網(wǎng)格，各區(qū)域交接面的節(jié)點一一對應(yīng)。經(jīng)無關(guān)性驗證后(圖3)，最終網(wǎng)格數(shù)量約為127萬，網(wǎng)格質(zhì)量為0.55。

1.3 湍流模型

Rohde等[11]結(jié)合硼混合實驗數(shù)據(jù)進(jìn)行湍流模型驗證，并指出不考慮浮升力影響的情況下，SST模型可準(zhǔn)確計算壓力容器內(nèi)硼濃度分布。SST模型對壁面附近的低雷諾數(shù)流動區(qū)域采用標(biāo)準(zhǔn)k-ω模型，對充分發(fā)展的自由流動區(qū)域采用標(biāo)準(zhǔn)k-ε模型，因此很大程度上降低了標(biāo)準(zhǔn)k-ω模型對自由剪切流的敏感程度[12-13]。SST模型的輸運方程如下：

圖2 壓力容器模型網(wǎng)格劃分Fig.2 Meshing of RPV model

圖3 網(wǎng)格無關(guān)性驗證Fig.3 Independence validation of mesh number

(7)

式中：Gk為由時均速度梯度引起的湍流動能產(chǎn)生項；Gω為湍流渦頻率產(chǎn)生項；Yk和Yω分別為湍流動能及湍流渦頻率的耗散項；Dω為交叉擴(kuò)散項。

為驗證CFD計算結(jié)果的準(zhǔn)確性，將數(shù)值模擬和實驗得到的入口橫截面平均質(zhì)量分?jǐn)?shù)進(jìn)行對比，結(jié)果如圖4所示?？煽闯?，入口橫截面的平均質(zhì)量分?jǐn)?shù)隨注射時間的增加而增大，在10 s后增長速度明顯減緩。CFD計算結(jié)果與實驗結(jié)果較為一致，平均相對誤差為3.41%。說明本文計算結(jié)果具有較高的準(zhǔn)確性。

圖4 數(shù)值模擬結(jié)果與實驗對比Fig.4 Comparison between CFD simulation and experiment

1.4 邊界條件設(shè)置

開啟混合物模塊，自定義硼酸溶液和水溶液的混合物，通過編寫UDF設(shè)置二者的組分占比，入口流速分別設(shè)置為0.2、0.3、0.4 m/s，硼酸質(zhì)量分?jǐn)?shù)c設(shè)置為：t≤0 s時，c=0；t>0 s時，c=0.1。系統(tǒng)運行壓力為0.1 MPa，混合物的物性設(shè)置為常壓下水的物性。忽略硼示蹤劑對流體的物性影響，流體密度為998.2 kg/m3，運動黏度為0.001 Pa·s。入口邊界條件設(shè)置為速度入口，出口邊界條件設(shè)置為壓力出口，壁面類型選擇無滑移壁面邊界條件。模擬混合過程前首先將壓力容器內(nèi)的全部流體定義為水，t>0 s時含水混合物自入口注入壓力容器。

2 數(shù)值模擬結(jié)果與分析

2.1 環(huán)形下降段內(nèi)流動特性分析

1) 二維速度場分析

圖5 水平方向速度矢量分布Fig.5 Transverse velocity vector

圖5為橫截面的速度矢量分布，對應(yīng)的入口雷諾數(shù)為11 448。可看出，在0 cm橫截面周向位置[10°，30°]和[-30°，-10°]區(qū)域存在較大尺度的渦結(jié)構(gòu)，由于入口流體與內(nèi)壁發(fā)生劇烈碰撞，流體呈放射狀散開，在靠近入口管道的區(qū)域形成回流并與環(huán)腔內(nèi)流體混合。Eltayeb等[14]通過可視化實驗在該區(qū)域觀測到了旋渦的形成與瓦解過程，旋渦的尺寸與位置存在一定的隨機(jī)性。而數(shù)值模擬結(jié)果僅能觀測到1個尺度較大的渦，這是由于兩方程渦黏模型基于平均化假設(shè)，在瞬態(tài)條件下無法計算出脈動值[15]。在單回路流動實驗的初始階段(0～2 s)，示蹤劑在上述區(qū)域的擴(kuò)散過程與速度矢量分布相對應(yīng)，染色劑與內(nèi)壁發(fā)生碰撞后形成逆流向外壁發(fā)展。隨后流體沿環(huán)形區(qū)域周向位置流動，期間不斷與環(huán)腔壁面發(fā)生碰撞，水平方向的速度隨周向位置的延伸而衰減，在[90°，130°]和[-130°，-90°]的區(qū)域最為明顯。在水平位置較低的-15 cm橫截面，[45°，135°]和[-45°，-135°]區(qū)域的流體流動不受上游渦結(jié)構(gòu)的影響，速度分布較為平穩(wěn)。

圖6、7分別為環(huán)形下降段內(nèi)水平方向和豎直方向的速度分布。如前文所述，流體進(jìn)入環(huán)腔后沿周向位置擴(kuò)散，期間不斷與壁面發(fā)生碰撞，水平方向動能不斷損失，水平方向的速度隨周向位置的增加而衰減。注意到-10～-25 cm截面水平方向速度先增加后減小，其中以圖6a最為明顯，這是由于靠近入口管道的區(qū)域受放射狀回流的影響。部分入口流體與內(nèi)壁碰撞后將沿豎直方向進(jìn)入環(huán)腔下部空間；環(huán)腔上部空間的入口流體受質(zhì)量力影響也會逐漸沿豎直方向沉降，可看出水平位置越低的截面，豎直方向的速度分量越大。

2) 三維流線圖分析

圖8為單回路穩(wěn)態(tài)流動條件下環(huán)形下降段和下腔室內(nèi)的流線圖。流體經(jīng)入口管進(jìn)入壓力容器后與環(huán)腔內(nèi)壁發(fā)生劇烈碰撞并形成回流，部分回流受邊界條件的束縛沿內(nèi)壁呈放射狀散開。隨后入口流體分成兩股沿著環(huán)形下降段周向位置流動，豎直方向的速度分量隨周向位置的延伸而增加，到達(dá)入口對側(cè)后環(huán)腔左右兩股流體開始匯合。

圖6 水平方向速度分布Fig.6 Transverse velocity distribution

圖7 豎直方向速度分布Fig.7 Longitudinal velocity distribution

圖8 環(huán)腔與下腔室內(nèi)流線圖Fig.8 Streamlines in down-comer and lower plenum

流線圖可反映微元流動介質(zhì)在流場中的軌跡。從圖8可看出，在入口正下方和環(huán)腔上部空間存在流動停滯區(qū)，該區(qū)域隨入口雷諾數(shù)的增加而收縮。在主泵恢復(fù)工作的工況下，主泵轉(zhuǎn)速由控制器驅(qū)動呈點動式增加而非線性增加，持續(xù)的低速轉(zhuǎn)動可能導(dǎo)致稀釋水團(tuán)滯留在環(huán)腔上部空間，與冷卻劑交混過程緩慢。

2.2 硼輸運特性分析

1) 橫向濃度場分布

將數(shù)值模擬得到的橫截面混合物質(zhì)量分?jǐn)?shù)分布與文獻(xiàn)[10]所報道的示蹤劑濃度分布進(jìn)行比較，結(jié)果如圖9、10所示。從數(shù)值模擬結(jié)果可看出，混合流體進(jìn)入環(huán)腔后與內(nèi)壁發(fā)生碰撞，隨后沿著周向位置發(fā)展并不斷與環(huán)腔內(nèi)流體混合；周向位置越大，混合物擴(kuò)散的速度越慢。在圖9、10中，數(shù)值模擬圖像中的比色卡為混合物份額，而實驗結(jié)果圖像的比色卡為染色劑(RhB)的濃度，二者的擴(kuò)散趨勢高度一致，進(jìn)一步驗證了模型的可靠性。SST模型能得到濃度擴(kuò)散整體趨勢，但無法獲得湍流擴(kuò)散的細(xì)節(jié)。在數(shù)值模擬結(jié)果中可觀測到入口位置內(nèi)壁和外壁之間存在一定的濃度梯度，而實驗結(jié)果顯示該區(qū)域的混合較為充分，濃度梯度較小。

a——數(shù)值模擬質(zhì)量分?jǐn)?shù)分布；b——LIF實驗染色劑濃度分布圖9 0 cm橫截面實驗與數(shù)值模擬結(jié)果對比Fig.9 Contrast between numerical simulation and experimental results at 0 cm cross section

a——數(shù)值模擬質(zhì)量分?jǐn)?shù)分布；b——LIF實驗染色劑濃度分布圖10 -10 cm橫截面實驗與數(shù)值模擬結(jié)果對比Fig.10 Contrast between numerical simulation and experimental results at -10 cm cross section

圖11為周向位置質(zhì)量分?jǐn)?shù)數(shù)值模擬與濃度實驗結(jié)果。從圖11a可觀測到，混合物進(jìn)入環(huán)腔后在入口位置持續(xù)出現(xiàn)1個質(zhì)量分?jǐn)?shù)的峰值，由于入口流體在水平方向通過的距離較短，且尚未與環(huán)形下降段內(nèi)外壁發(fā)生碰撞，混合特性較差。隨后高濃度硼溶液沿周向位置流動并不斷與水發(fā)生混合，環(huán)腔內(nèi)混合物所占區(qū)域不斷增加，如圖11b所示。可注意到，在t>5 s后，周向角為[30°，110°]和[-30°，-110°]區(qū)域的質(zhì)量分?jǐn)?shù)平穩(wěn)增長，說明該區(qū)域的混合物已與水充分混合。從圖11可看出，周向位置質(zhì)量分?jǐn)?shù)不斷展平，質(zhì)量分布區(qū)域均勻。由于壓力容器模型加工精度所限環(huán)腔兩側(cè)濃度分布輕微不對稱，整體上，數(shù)值模擬結(jié)果與實驗結(jié)果較為一致，實際流動過程染色劑發(fā)展的邊緣位置存在一定脈動。

2) 三維濃度場分布

圖12為硼酸溶液注入環(huán)形下降段后的質(zhì)量分?jǐn)?shù)分布。含硼溶液進(jìn)入壓力容器后呈放射狀散開，并繞環(huán)腔壁邊界流動；混合物到達(dá)周向位置較大的區(qū)域后水平方向的速度衰減，其輸運特性過渡為沿豎直方向的沉降，以圖12b中t=24 s時刻最為明顯。從圖12可看出，單回路流動混合時，入口正下方存在1個流動停滯區(qū)，在整個混合過程中模擬硼酸溶液無法進(jìn)入到流動停滯區(qū)?？梢?，同等時間下，入口混合流體雷諾數(shù)越大，混合物所分布的區(qū)域相應(yīng)越大，整個環(huán)腔內(nèi)的流通停滯區(qū)相應(yīng)收縮。文獻(xiàn)[10]所報道的示蹤劑濃度擴(kuò)散特性與該現(xiàn)象一致。

a——周向位置質(zhì)量分?jǐn)?shù)分布數(shù)值模擬；b——周向位置濃度分布LIF實驗結(jié)果圖11 數(shù)值模擬與實驗結(jié)果定量對比Fig.11 Quantitative contrast between numerical simulation and experimental results

a——Re=5 833；b——Re=11 448圖12 壓力容器三維質(zhì)量分?jǐn)?shù)分布Fig.12 Three-dimensional mass fraction distribution in RPV

3 結(jié)論

本文基于SST模型對環(huán)形下降段內(nèi)的硼混合問題進(jìn)行了分析，得到如下結(jié)論。

1) 單回路條件下，環(huán)形下降段上部空間水平方向的速度梯度較大，下部空間流體流動較為穩(wěn)定；環(huán)形下降段內(nèi)存在流動停滯區(qū)，且該區(qū)域隨入口雷諾數(shù)的增加而收縮。

2) 含硼溶液從單一入口進(jìn)入壓力容器后將沿環(huán)形下降段周向輸運，期間水平方向動能不斷損失，在環(huán)腔內(nèi)的輸運特性逐漸過渡為沿豎直方向的沉降。