堵流工況下棒束子通道流場實驗與模擬研究

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

(中國核動力研究設計院 核反應堆熱工水力技術重點實驗室，四川 成都 610213)

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

1 實驗裝置

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

圖1 棒束組件內流場測量實驗裝置流程圖Fig.1 Schematic diagram of rod bundle flow field measurement facility

圖2a所示為實驗本體結構示意圖，實驗本體主體為四面可視的矩形通道結構，進出口有穩流裝置，內部放置5×5棒束組件。組件由支撐段和測量段(圖2b)組成，支撐段為硬質材料，測量段為設計外徑9.5 mm、中心距12.6 mm的薄壁FEP管，以最大程度接近水的折射率。堵流件安裝在FEP管外，堵流區域為3×3棒束區域(圖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以上，以最大程度減弱入口和出口效應，測量區域約為20D。圖2c中綠色虛線所示位置為拍攝豎直流場所在截面位置示意，具體為棒中心線和通道中心線。

a——本體結構；b——測量段；c——棒束截面結構圖2 實驗本體及棒束組件結構示意圖Fig.2 Schematic diagram of test facility and rod bundle

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

2 數值模擬方法

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

(1)

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

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

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

本次模擬針對實驗工況，所涉及的模型為冷態流動，模擬為穩態模擬，主要求解x、y、z方向的速度以及壓力。模擬邊界條件為流量入口，流量、壓力、溫度根據實驗測量結果設置，模擬中迭代計算到殘差小于10-5或不再減小為止。

3 結果及分析

3.1 流速測量結果

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

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

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

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

(2)

圖7為根據實驗測量結果得到的通道C1、C2、C3、C4在不同高度處的流速不均勻指數，其中虛線所在位置為堵流件，軸向位置0為堵流件軸向中心位置。可看出，通道C1和C4的流速不均勻指數變化較小，僅在虛線所示堵流件區域出現了較小的增加，并很快恢復正常，表明其受堵流件的影響較小，冷卻劑可正常流動。而通道C2和C3流速不均勻指數在堵流件處出現了較大的增加，而后緩慢減小；且流速不均勻指數在堵流件上游1D處已開始增加，表明其對流場的影響在上游1D處已開始。雖然在堵流件下游約3D后流速不均勻指數下降較緩，但仍未恢復到堵流前水平。

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

3.2 子通道流量結果

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

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

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

(3)

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

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

3.3 模擬結果及分析

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

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

a——堵流件下游1D；b——堵流件下游4D圖11 不同距離處流速的實驗與模擬結果對比Fig.11 Comparison of experimental and simulated velocities at different locations

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

(4)

式中：vi,Exp為流速測量值；vi,Simu為對應位置流速的模擬值。

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

RMSE越小，表明模擬值與實驗測量值偏差越小，模擬不確定度越小。麒麟KILI V1.0和商用CFD軟件盡量保證網格總數和最小網格尺寸相同。可看出，整體上1D處的模擬值不確定度大于4D處。1D處回流較強，對湍流模擬準確度要求較高，因此采用RANS模型難免對局部細節模擬準確度有限。模擬不確定度整體隨網格總數增加而降低，由于網格總數增加主要是對堵流件附近區域網格進行加密，因此1D處不確定度下降較快。在1D處RMSE隨著網格總數由30萬逐漸增加到220萬不斷降低，在4D處網格總數達到140萬后RMSE下降很少，表明在1D處需要更精細的局部網格才能實現網格無關，在4D處因為湍流強度降低，網格大小降低到一定程度后即實現了網格無關。總體而言，麒麟KILI V1.0在1D處與商業數值分析軟件不確定度相當，在4D處部分網格總數下麒麟KILI V1.0優于商業軟件。網格總數較小的20萬和40萬，商業CFD軟件模擬準確度絕對值上更低，主要原因是商業CFD軟件采用貼體網格和局部邊界層加密，更有利于湍流模擬。麒麟KILI V1.0采用IST技術的網格，網格劃分不區分流固區域，網格不貼體，且采用局部整體網格加密，因此其湍流模擬準確度會降低。在網格加密到一定程度后，網格尺寸對模擬準確度影響較小，所以麒麟KILI V1.0在網格總數增加后與商業CFD軟件差別減小，甚至更優。

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

4 結論

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