堵流工況下棒束子通道流場實(shí)驗(yàn)與模擬研究

陳登高，畢景良，黃彥平，袁德文，昝元峰，徐建軍

(中國核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆熱工水力技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

核電作為一種高效的清潔能源受到了能源行業(yè)的青睞，并逐漸成為一種重要的能源生產(chǎn)方式。壓力容器中的堆芯是熱能產(chǎn)生的部位，也是反應(yīng)堆的關(guān)鍵部件之一，其熱工水力特性直接關(guān)系到燃料元件的冷卻效率及安全問題。在極端情況下，燃料組件子通道可能因?yàn)榘羰[脹、異物堆積等發(fā)生通道堵流現(xiàn)象，傳熱惡化，甚至超溫、堆芯融化等危險(xiǎn)工況，因此有必要開展相關(guān)研究。目前，針對(duì)通道堵流的研究主要集中在矩形通道的研究，如郭玉川[1]研究了板元件在堵塞情況下的流場分布特性，Xu等[2]及Yuan等[3]研究了板元件加熱情況下堵流時(shí)流動(dòng)阻力、傳熱等特性。此外，采用模擬方法研究子通道堵塞特性也是常用的方法，如李健全等[4]使用系統(tǒng)程序研究了堵塞情況下的反應(yīng)堆整體熱工水力變化；Lee等[5]、Guo等[6]、Salama等[7-9]則使用CFD分析軟件計(jì)算了堵流情況下矩形通道的流動(dòng)和傳熱特性。目前已有大量研究人員使用激光粒子測速(PIV)方法開展了子通道流場測量[12-15]，但僅有很少工作是堵流情況下棒束子通道流場的PIV測量[10-11]，且尚無堵塞情況下子通道流量的測量及分析。堵流下的子通道流量分布是建立系統(tǒng)程序的重要數(shù)據(jù)，因此有必要開展相關(guān)的研究。本文利用PIV方法測量5×5棒束堵塞情況下子通道的流場分布特性及子通道流量分布特性，并基于自研數(shù)值分析軟件麒麟KILI V1.0和商用CFD軟件進(jìn)行堵流現(xiàn)象的數(shù)值模擬，結(jié)合實(shí)驗(yàn)結(jié)果分析棒束子通道在局部通道堵塞情況流場、子通道流動(dòng)分布等特性。

1 實(shí)驗(yàn)裝置

本文首先采用可視化實(shí)驗(yàn)方法對(duì)堵流情況下5×5棒束組件子通道流量分配進(jìn)行了測量，實(shí)驗(yàn)回路如圖1所示，主要包括實(shí)驗(yàn)本體、水箱、主泵、流量計(jì)、閥門和測控儀表等。主泵從水箱中取水加壓后一路進(jìn)入實(shí)驗(yàn)本體，一路回流至水箱，兩路都設(shè)置調(diào)節(jié)閥以調(diào)節(jié)進(jìn)入實(shí)驗(yàn)本體的水流量。通過文丘里流量計(jì)測量進(jìn)入本體的水流量，通過壓力變送器和熱電偶測量水壓和溫度，所有測量儀表都進(jìn)行了檢定或校準(zhǔn)以滿足相應(yīng)的測量精度要求。

圖1 棒束組件內(nèi)流場測量實(shí)驗(yàn)裝置流程圖Fig.1 Schematic diagram of rod bundle flow field measurement facility

圖2a所示為實(shí)驗(yàn)本體結(jié)構(gòu)示意圖，實(shí)驗(yàn)本體主體為四面可視的矩形通道結(jié)構(gòu)，進(jìn)出口有穩(wěn)流裝置，內(nèi)部放置5×5棒束組件。組件由支撐段和測量段(圖2b)組成，支撐段為硬質(zhì)材料，測量段為設(shè)計(jì)外徑9.5 mm、中心距12.6 mm的薄壁FEP管，以最大程度接近水的折射率。堵流件安裝在FEP管外，堵流區(qū)域?yàn)?×3棒束區(qū)域(圖2c)，中間通道堵流比例最大，堵流中心通道(圖2c中紅色部分，C2-3、C2-4、C3-3、C3-4)堵流比例為72%，堵流邊通道1(圖2c中黃色部分，C1-3、C1-4、C2-2、C3-2、C4-3、C4-4)堵流比例約為20%，堵流通道2(圖2c中綠色部分，C1-2、C4-2)堵流比例約為10%。堵流件安裝位置距入口大于40D，距出口距離也在20D以上，以最大程度減弱入口和出口效應(yīng)，測量區(qū)域約為20D。圖2c中綠色虛線所示位置為拍攝豎直流場所在截面位置示意，具體為棒中心線和通道中心線。

a——本體結(jié)構(gòu)；b——測量段；c——棒束截面結(jié)構(gòu)圖2 實(shí)驗(yàn)本體及棒束組件結(jié)構(gòu)示意圖Fig.2 Schematic diagram of test facility and rod bundle

PIV面激光以xz平面照射測量區(qū)域，高速相機(jī)則沿y方向拍攝。采用跨幀模式進(jìn)行圖像采集，然后使用專業(yè)軟件進(jìn)行流速場計(jì)算，即獲得了測量區(qū)域的流場分布。相機(jī)拍攝幀率為15 Hz，每次針對(duì)同一工況拍攝200組以上照片，然后求取平均流場。本文還基于PIV測量結(jié)果采用數(shù)值積分方法獲得了棒束子通道流量結(jié)果，具體測量和計(jì)算方法參考文獻(xiàn)[16]。

2 數(shù)值模擬方法

數(shù)值模擬基于麒麟KILI V1.0進(jìn)行，該軟件采用浸入界面(IST)方法進(jìn)行網(wǎng)格劃分，基于NS方程求解流場、溫度場、壓力場等。該軟件支持常見CAD軟件三維幾何模型的導(dǎo)入，然后進(jìn)行IST方法的笛卡爾網(wǎng)格劃分，而后進(jìn)行邊界條件設(shè)置、模型設(shè)置、求解設(shè)置等，最后進(jìn)行迭代求解和后處理。為了保證透光性，實(shí)驗(yàn)所用FEP管壁厚僅約0.5 mm，因此易變形導(dǎo)致其定位與設(shè)計(jì)不同。為了保證數(shù)值模擬所建模型最大程度地還原實(shí)驗(yàn)棒束，基于PIV所拍攝的棒束圖進(jìn)行圖像處理并測量了棒束的真實(shí)尺寸，測量原理示意圖如圖3所示，基于PIV拍攝通道截面圖，測量棒束通道參考尺寸LRef和待測量尺寸LB對(duì)應(yīng)的像素?cái)?shù)量PRef和PB，然后進(jìn)行長度計(jì)算：

(1)

圖3 棒束安裝后實(shí)際尺寸測量原理示意圖Fig.3 Schematic diagram of size measurement of installed rod bundle

圖4為采用浸入界面方法時(shí)網(wǎng)格劃分示意圖。圖4a為浸入界面方法的基本原理，即在網(wǎng)格劃分時(shí)不區(qū)分流體與固體域，使用流固界面符號(hào)距離函數(shù)φs進(jìn)行數(shù)值地區(qū)分流體與固體域。在求解NS方程時(shí)，通過引入Heaviside函數(shù)H(φs)自動(dòng)地區(qū)別流體與固體域。該函數(shù)在流體內(nèi)等于1，在固體內(nèi)等于0，在界面上等于0.5。帶堵流件的棒束幾何由三維CAD軟件根據(jù)實(shí)驗(yàn)裝置尺寸測量結(jié)果建立，然后導(dǎo)入麒麟KILI V1.0中進(jìn)行網(wǎng)格劃分，圖4b的固體域即是繪制的棒束三維幾何。劃分網(wǎng)格時(shí)直接根據(jù)實(shí)驗(yàn)裝置矩形流道大小繪制矩形流體域，即圖4c中的背景網(wǎng)格，這樣三維幾何就浸沒在背景網(wǎng)格中(流體域)。然后基于重疊網(wǎng)格技術(shù)在堵流件區(qū)域添加加密的子網(wǎng)格對(duì)堵流塊進(jìn)行加密，以提高堵流件區(qū)域流場解析精度。本次模擬主要針對(duì)堵流件上下游約9D的范圍，為了避免入口效應(yīng)，也繪制了堵流件上游10D的網(wǎng)格。堵流件對(duì)子通道流場有很強(qiáng)的影響，局部網(wǎng)格精細(xì)程度對(duì)堵流件下游臨近流場模擬精細(xì)程度影響很大，因此網(wǎng)格數(shù)量對(duì)流場模擬結(jié)果影響很大，因此本文也探索了不同網(wǎng)格數(shù)量對(duì)模擬結(jié)果的影響。

圖4 采用浸入界面方法時(shí)網(wǎng)格劃分示意圖Fig.4 Grid generation of simulated rod bundle using IST method

本次模擬針對(duì)實(shí)驗(yàn)工況，所涉及的模型為冷態(tài)流動(dòng)，模擬為穩(wěn)態(tài)模擬，主要求解x、y、z方向的速度以及壓力。模擬邊界條件為流量入口，流量、壓力、溫度根據(jù)實(shí)驗(yàn)測量結(jié)果設(shè)置，模擬中迭代計(jì)算到殘差小于10-5或不再減小為止。

3 結(jié)果及分析

3.1 流速測量結(jié)果

本實(shí)驗(yàn)中進(jìn)口流量穩(wěn)定在約14 t/h，對(duì)應(yīng)的水平截面平均質(zhì)量流速約1 566 kg/(m2·s)，平均流速約1.6 m/s。實(shí)驗(yàn)中測量了堵流件上下游約7D范圍內(nèi)棒束中心豎直截面、子通道豎直中心截面的流場。測量得到的典型豎直截面流速使用截面平均流速進(jìn)行了歸一化處理，即實(shí)際流速/截面平均流速，結(jié)果如圖5、6所示。圖5、6所示結(jié)果皆為100組PIV測量結(jié)果的平均值，對(duì)應(yīng)總時(shí)長約6.5 s，橫坐標(biāo)使用棒束間距P進(jìn)行歸一化處理，縱坐標(biāo)使用棒束外徑D進(jìn)行歸一化處理?？煽闯觯跊]有堵流件的正常狀態(tài)下，堵流件所在位置上下游流場基本都以本子通道的豎直流動(dòng)為主，流場未發(fā)生偏斜，不同子通道流速分布相似。在堵流情況下，C1和C2都出現(xiàn)了流場偏斜和流速不均的情況。對(duì)照?qǐng)D2c可知C1子通道最大堵流比例約20%，C1-3和C1-4子通道下游約3D范圍內(nèi)流速小于截面平均流速的50%，可見僅堵塊下游少量流場受影響。對(duì)于C2-3和C2-4，子通道堵塞比例為72%，堵流件下游2D范圍內(nèi)流速整體小于平均流速的20%，超過4D后才恢復(fù)到約50%。C2-3和C2-4子通道在堵流件下游也出現(xiàn)了明顯的局部回流和漩渦，局部流速小于截面平均流速的10%，堵塞現(xiàn)象十分明顯。

a——C1通道；b——C2通道圖5 非堵流狀態(tài)下截面歸一化流速測量結(jié)果Fig.5 Typical normalized velocity from experiment measuring under non-blocking condition

a——C1通道；b——C2通道；c——R2棒束中心截面；d——R3棒束中心截面圖6 堵流工況下歸一化流速測量結(jié)果Fig.6 Typical normalized velocity from experiment measuring under blocking condition

為了直觀地分析堵流對(duì)子通道流速的影響，使用流速不均勻指數(shù)進(jìn)行量化。該指數(shù)等于指定豎直截面高度z處流速絕對(duì)值與截面平均流速的均方根，具體計(jì)算公式為：

(2)

圖7為根據(jù)實(shí)驗(yàn)測量結(jié)果得到的通道C1、C2、C3、C4在不同高度處的流速不均勻指數(shù)，其中虛線所在位置為堵流件，軸向位置0為堵流件軸向中心位置?？煽闯?，通道C1和C4的流速不均勻指數(shù)變化較小，僅在虛線所示堵流件區(qū)域出現(xiàn)了較小的增加，并很快恢復(fù)正常，表明其受堵流件的影響較小，冷卻劑可正常流動(dòng)。而通道C2和C3流速不均勻指數(shù)在堵流件處出現(xiàn)了較大的增加，而后緩慢減?。磺伊魉俨痪鶆蛑笖?shù)在堵流件上游1D處已開始增加，表明其對(duì)流場的影響在上游1D處已開始。雖然在堵流件下游約3D后流速不均勻指數(shù)下降較緩，但仍未恢復(fù)到堵流前水平。

圖7 流速不均勻指數(shù)隨z的變化Fig.7 Velocity uneven factor at different z values

3.2 子通道流量結(jié)果

本研究使用PIV測量獲得了子通道流量結(jié)果，具體測量和計(jì)算方法參見文獻(xiàn)[16]，本次流量測量范圍為堵流件中心上游約2D至下游約5D范圍。圖8為獲得的堵流件中心下游不同距離處截面豎直方向質(zhì)量流速的分布，其中質(zhì)量流速以平均質(zhì)量流速1 500 kg/(m2·s)進(jìn)行了歸一化處理。從圖8可看出，堵流件上游1D處水平截面質(zhì)量流速分布仍較均勻，而在堵流件下游0.5D處，堵流中心通道(C2-3、C2-4、C3-3、C3-4)質(zhì)量流速出現(xiàn)了較大降低，上述子通道質(zhì)量流速普遍在平均值的30%以下。堵流區(qū)域臨近的非堵流區(qū)域則出現(xiàn)了明顯的質(zhì)量流速增加，其最大質(zhì)量流速超過平均值的2倍，表明堵流導(dǎo)致的子通道質(zhì)量流速偏移十分明顯。質(zhì)量流速在堵流中心下游4D后恢復(fù)較明顯，堵流中心質(zhì)量流速已恢復(fù)到平均值的80%左右。圖8a中部分棒束周圍出現(xiàn)了局部較低的流速測量結(jié)果，通過分析原始PIV圖片發(fā)現(xiàn)，導(dǎo)致這一現(xiàn)象的主要原因是部分棒束在加工制造中出現(xiàn)了難以察覺的刮傷，在激光照射下出現(xiàn)局部亮度過高的情況，進(jìn)而導(dǎo)致這些局部區(qū)域的示蹤粒子被掩蓋，PIV相機(jī)拍攝不到清晰的移動(dòng)粒子，進(jìn)而計(jì)算為低流速區(qū)域。在實(shí)驗(yàn)中進(jìn)行了多次更換新棒束，仍在局部出現(xiàn)了上述現(xiàn)象，但對(duì)真實(shí)流速分布沒有影響，也不影響整體流場結(jié)果的判斷。圖8b中C1-3和C1-4子通道質(zhì)量流速明顯高于C4-3和C4-4子通道的現(xiàn)象，原因是安裝后的帶堵流件棒束在R1-X側(cè)與流道壁面的間隙(測量值約4.07 mm)大于R5-X側(cè)(測量值約2.17 mm)，導(dǎo)致在堵流截面處流體經(jīng)C1-3和C1-4側(cè)繞流更多，進(jìn)而導(dǎo)致質(zhì)量流速更多。該現(xiàn)象隨著遠(yuǎn)離堵流區(qū)域而明顯減小，這一差異對(duì)堵流件對(duì)冷卻劑流動(dòng)的改變趨勢判斷影響不大。

a——z=-1D；b——z=0.5D；c——z=1D；d——z=4D圖8 堵流件下游不同距離處截面質(zhì)量流速分布Fig.8 Mass velocity distribution of blockage downstream

為了進(jìn)一步量化比較堵流對(duì)子通道流量的影響，本研究基于測量所得子通道流量計(jì)算了5×5棒束的4×4子通道截面流量不均勻指數(shù)n，其計(jì)算公式為：

(3)

圖9為堵流件下游不同距離處子通道流量分布系數(shù)，其中折線為不同子通道歸一化流量，柱狀圖為子通道截面流量不均勻指數(shù)，0所在位置為堵流件高度中心。可看出，在堵流區(qū)域上游1D處子通道截面流量不均勻指數(shù)明顯增加，表明子通道流量分配已受到影響，堵流件下游0.5D范圍內(nèi)流量不均勻指數(shù)出現(xiàn)了突增，在1D后不均勻程度有了明顯恢復(fù)。圖9中曲線為堵流中心通道(C2-3、C2-4、C3-3、C3-4)、堵流邊通道(C2-2、C3-2)及部分未堵流通道(C2-1、C3-1)的子通道歸一化流量。與流量不均勻指數(shù)類似，3種類型通道都是在堵流件上游1D距離開始出現(xiàn)了流量變化，其中堵流中心通道流量降低，堵流邊通道和未堵流通道則上升。在堵流件下游0.5D范圍內(nèi)，堵流中心通道流量出現(xiàn)陡降，最小流量僅為理論流量的25%左右，而堵流邊通道流量則突增，最大流量約為理論流量的1.7倍，未堵流通道最大則約為1.3倍。由圖9可知，中心堵流對(duì)所在子通道有強(qiáng)烈的限流作用，對(duì)鄰近的子通道則有強(qiáng)烈的增加流量作用，對(duì)鄰近未堵流通道有增流作用。3種子通道的流量變化都在堵流下游1D處出現(xiàn)了明顯的緩解。

圖9 堵流件下游子通道截面流量不均勻指數(shù)Fig.9 Sub-channel flow rate uneven index of blockage downstream

3.3 模擬結(jié)果及分析

使用前述數(shù)值模擬方法進(jìn)行數(shù)值迭代求解可得到5×5棒束在堵流情況下的模擬結(jié)果，本研究也使用商業(yè)CFD軟件(ANSYS FLUENT 15.0)進(jìn)行了模擬，網(wǎng)格繪制時(shí)盡量保證二者網(wǎng)格總數(shù)和最小網(wǎng)格尺寸相同，此外二者湍流模型(標(biāo)準(zhǔn)k-ε模型)、邊界條件、求解算法等相同，以對(duì)比麒麟KILI V1.0與商業(yè)CFD軟件。圖10為模擬所得C1、C2子通道中心豎直截面流場的典型模擬結(jié)果，其中流速為使用截面平均流速歸一化處理的流速，其中白色框部分為堵流件所在位置。C2子通道位于堵流中心通道，堵塞比例較高，模擬結(jié)果中堵流件下游形成了低速區(qū)，這與實(shí)驗(yàn)測量結(jié)果吻合。比較麒麟KILI V1.0和商用CFD軟件的模擬結(jié)果可發(fā)現(xiàn)二者流場分布十分接近，表明麒麟KILI V1.0可很好地模擬流場分布。

a——C1，KILI V1.0；b——C2，KILI V1.0；c——C1，CFD軟件；d——C2，CFD軟件圖10 模擬所得堵流情況下歸一化流速云圖Fig.10 Typical normalized velocity contour from simulation under blocking condition

為了更加直觀地對(duì)比麒麟KILI V1.0和商用CFD軟件差別，提取了堵流件下游不同距離處實(shí)驗(yàn)與模擬的結(jié)果，如圖11所示(C2通道中心截面)。可看出，在兩種距離處模擬結(jié)果可較好地反映實(shí)驗(yàn)結(jié)果，其中1D處堵流中心通道流速較低，實(shí)驗(yàn)測量中有強(qiáng)烈回流；麒麟KILI V1.0與商用CFD軟件在細(xì)節(jié)流速模擬上還原實(shí)驗(yàn)測量值準(zhǔn)確程度相當(dāng)。

a——堵流件下游1D；b——堵流件下游4D圖11 不同距離處流速的實(shí)驗(yàn)與模擬結(jié)果對(duì)比Fig.11 Comparison of experimental and simulated velocities at different locations

圖12為模擬結(jié)果的均方根誤差(RMSE)，其計(jì)算公式為：

(4)

式中：vi,Exp為流速測量值；vi,Simu為對(duì)應(yīng)位置流速的模擬值。

a——z=1D；b——z=4D圖12 模擬值均方根誤差Fig.12 RMSE for different simulation methods

RMSE越小，表明模擬值與實(shí)驗(yàn)測量值偏差越小，模擬不確定度越小。麒麟KILI V1.0和商用CFD軟件盡量保證網(wǎng)格總數(shù)和最小網(wǎng)格尺寸相同?？煽闯?，整體上1D處的模擬值不確定度大于4D處。1D處回流較強(qiáng)，對(duì)湍流模擬準(zhǔn)確度要求較高，因此采用RANS模型難免對(duì)局部細(xì)節(jié)模擬準(zhǔn)確度有限。模擬不確定度整體隨網(wǎng)格總數(shù)增加而降低，由于網(wǎng)格總數(shù)增加主要是對(duì)堵流件附近區(qū)域網(wǎng)格進(jìn)行加密，因此1D處不確定度下降較快。在1D處RMSE隨著網(wǎng)格總數(shù)由30萬逐漸增加到220萬不斷降低，在4D處網(wǎng)格總數(shù)達(dá)到140萬后RMSE下降很少，表明在1D處需要更精細(xì)的局部網(wǎng)格才能實(shí)現(xiàn)網(wǎng)格無關(guān)，在4D處因?yàn)橥牧鲝?qiáng)度降低，網(wǎng)格大小降低到一定程度后即實(shí)現(xiàn)了網(wǎng)格無關(guān)。總體而言，麒麟KILI V1.0在1D處與商業(yè)數(shù)值分析軟件不確定度相當(dāng)，在4D處部分網(wǎng)格總數(shù)下麒麟KILI V1.0優(yōu)于商業(yè)軟件。網(wǎng)格總數(shù)較小的20萬和40萬，商業(yè)CFD軟件模擬準(zhǔn)確度絕對(duì)值上更低，主要原因是商業(yè)CFD軟件采用貼體網(wǎng)格和局部邊界層加密，更有利于湍流模擬。麒麟KILI V1.0采用IST技術(shù)的網(wǎng)格，網(wǎng)格劃分不區(qū)分流固區(qū)域，網(wǎng)格不貼體，且采用局部整體網(wǎng)格加密，因此其湍流模擬準(zhǔn)確度會(huì)降低。在網(wǎng)格加密到一定程度后，網(wǎng)格尺寸對(duì)模擬準(zhǔn)確度影響較小，所以麒麟KILI V1.0在網(wǎng)格總數(shù)增加后與商業(yè)CFD軟件差別減小，甚至更優(yōu)。

圖13為堵流件所在區(qū)域流場特征模擬結(jié)果。流線圖表明堵流件下游區(qū)域產(chǎn)生明顯的回流區(qū)，導(dǎo)致流速突降，這可能導(dǎo)致棒束冷卻不足；矢量圖用x方向流速值進(jìn)行標(biāo)志顏色，紅色表示正向，藍(lán)色表示負(fù)向，可看出，堵流件的阻塞導(dǎo)致所在子通道冷卻劑向周圍子通道流動(dòng)，這是堵流中心子通道流量減少，堵流邊通道流量增加的主要原因。

4 結(jié)論

本文基于棒束子通道流場可視化實(shí)驗(yàn)裝置和5×5棒束子通道堵流實(shí)驗(yàn)本體獲得了子通道堵流72%比例下的流場分布及子通道流量結(jié)果，并使用麒麟KILI V1.0對(duì)堵流情況下的子通道進(jìn)行了數(shù)值模擬分析。結(jié)果表明：1) 在子通道堵流約72%的情況下，堵塞子通道上游約1D范圍流場開始受到明顯影響，堵流通道下游0.5D距離內(nèi)堵流影響最強(qiáng)烈，堵塞中心子通道流速僅為平均流速的約30%，子通道流量僅約為理論流量的25%，堵流影響在1D后出現(xiàn)明顯緩解，堵塞中心子通道流量恢復(fù)到理論流量的60%以上；2) 基于浸入界面方法的麒麟KILI V1.0能對(duì)堵流情況下的棒束子通道流場進(jìn)行高準(zhǔn)度的模擬，其模擬結(jié)果與商用CFD軟件相當(dāng)。模擬結(jié)果表明，堵流件下游的局部回流是導(dǎo)致流速降低的主要原因，堵流通道阻力增大，使得冷卻劑向堵塞通道臨近通道分流，導(dǎo)致堵塞通道流量降低，堵流中心通道相鄰子通道流量增大。