(中國核動力研究設計院，四川 成都 610213)

反應堆控制棒驅動線作為反應堆內具有相對運動的設備單元。一般由燃料組件、導向組件、控制棒組件、控制棒驅動機構等部分組成。由于控制棒組件對堆芯反應性控制起決定性作用，因此一個合格的驅動線設計必須保證在緊急情況下，控制棒能夠在規(guī)定時間內從反應堆內任意位置迅速插入堆芯，迫使反應堆停堆，防止事態(tài)惡化。因此，一直以來落棒時間也是控制棒驅動線設計首要考慮的指標，同時也是核電站安全分析的重要參考。由于物理樣機造價昂貴，因此大批學者為準確計算落棒時間進行了大量研究工作[1-7]，這些研究工作將控制棒驅動線簡化為一維水力模型，取得了大量的成果。這些成果作為控制棒驅動線設計的重要參考已經在核行業(yè)取得廣泛應用，并形成了許多專用軟件，例如CIGAL 軟件(法國)、DROP 軟件(美國)等。

但是，由于這些程序采用面向過程的編程思路進行編制，當驅動線結構發(fā)生變化時往往需要修改大量代碼，擴展性較差。且由于采用一維水力模型，很多細節(jié)，包括流道形狀變化均依靠引入實驗修正系數來進行計算，因此對實驗依賴性較高。隨著CFD計算水平的提高和計算機技術的進步，一些學者選擇采用CFD動網格方法對控制棒落棒行為進行仿真分析，例如，肖聰等基于 FLUENT 動網格技術對單根圓柱形控制棒建立了三維流體仿真模型來進行落棒行為仿真分析。同樣，基于流體動力學(CFD)的動網格技術，在新的堆型上模擬了十字形控制棒組件在控制棒導向組件內的運動行為[8-9]。但是由于三維模型計算量過大，該方法只對驅動線的局部流域進行了仿真建模分析。到目前為止還沒有針對完整驅動線的有限元計算分析。

本文提出了一種基于動網格的反應堆控制棒落棒行為分流域耦合仿真方法，該方法分別建立了控制棒單棒和驅動桿的二維軸對稱模型，保證控制棒和驅動桿對應流域的網格能夠根據落棒運動規(guī)律自適應地變化。兩個流域在每個時間步長內交換流體阻力計算結果，并根據運動方程求解得到的速度來更新控制棒和驅動桿的運動狀態(tài)，實現分流域耦合。本方法通用性好，計算中考慮了驅動線流道形狀的影響，且在計算時間和求解精度之間取得了良好的折中，此外，本方法在計算條件允許的情況下，還能較容易地擴展到三維模型。

1 流體仿真模型

1.1 網格模型

由于反應堆控制棒驅動線長度較長，而控制棒和驅動桿周圍的流體間隙相對小很多，要獲得與軸向長度和徑向間隙長度相適應的網格尺寸，網格量將相當大。因此，為了減少網格模型大小，本文選用二維軸對稱模型來對驅動線流場進行模擬，且為了保證計算精度和動網格更新，整個流域采用四邊形網格進行網格劃分。

驅動線流場分為控制棒流域和驅動桿流域兩個部分，忽略星形架等結構，并對流體域流道作了適當簡化。以圖1所示的控制棒流域為例，控制棒與導向管之間的流體域由Interface交界面分成了靜止區(qū)域和運動區(qū)域兩個部分，其中在運動區(qū)域，控制棒前端和后端分別定義了兩個剛性面，其邊界類型為Interior，主要作用是與Interface交界面將控制棒包圍起來，所包圍區(qū)域的網格可以隨著控制棒進行剛性移動，而剛性面之外的運動區(qū)域網格則采用動態(tài)網格層變的方法進行網格更新，分割因子設置為 0.4，合并因子設為 0.2。這樣就保證了控制棒能夠在導向管中沿軸向自由移動而不會產生負網格。

圖1 網格模型Fig.1 Mesh model

1.2 計算模型

本文采用FLUENT 14.5版本進行仿真計算。由于采用二維軸對稱模型，因此計算時求解的是二維軸對稱的連續(xù)性方程和動量方程，湍流模型采用k-epsilon湍流模型，使用SIMPLE 算法進行求解。

固體壁面均設置為固壁邊界，其邊界設置為無滑移條件，近壁面應用標準壁面函數。

壓力插值格式選擇標準格式，動量選擇二階迎風格式，湍動能及湍流耗散率皆選取一階迎風格式。

反應堆冷態(tài)工況下，溫度為 20 ℃，壓力為101 kPa。整個流域與外界導通的部分，包括緩沖段的流水孔等，邊界條件均設置為壓力出口，表壓為0。

2 分流域耦合方法

由于網格量和軸對稱條件的限制，因此將控制棒驅動線的流場分為控制棒單棒和驅動桿兩個部分。兩個流場迭代計算同時進行，在每個時間步內進行數據交換以保持同步。

基本計算流程如圖2所示，計算過程如下：

1)控制棒單棒流體域計算為主計算流程，在每個時間步長內首先計算當前時間步長Δt，通過DEFINE_DELTAT宏修改時間步長，并輸出Δt到文件time.txt，提供給驅動桿流域計算，驅動桿流域讀入time.txt后同樣通過DEFINE_DELTAT宏修改時間步長；

2)驅動桿流域積分計算驅動桿流體阻力F_d，輸出輸出F_d到文件Force.txt文件，并等待更新速度文件；

3)控制棒單棒流體域讀入文件Force.txt獲得F_d，并積分計算控制棒流體阻力F_c(該值為單棒流體阻力乘以控制棒數)，通過如下方程求解速度改變量：

Δv=(M×g+F_d+F_c)×Δt/M

(1)

式中：M——控制棒和驅動桿的總重;

g——重力加速度。

同時更新速度，并寫入文件vel.txt，根據更新速度值利用DEFINE_GRID_MOTION宏更新動網格區(qū)域網格。

v=v+Δv

(2)

4)驅動桿流域讀入vel.txt，得到速度值，利用DEFINE_GRID_MOTION宏更新動網格區(qū)域網格。

5)驅動桿流域和控制棒單棒流體域利用DEFINE_CG_MOTION宏更新值更新Interior邊界面(剛性面)的速度。

6)數據交互完成并更新網格之后，兩個流域同時開始迭代計算，迭代完成后進入下一時間步直到計算結束。

圖2 分流域耦合方法流程圖Fig. 2 Flow chart of multi-field coupling method

3 方法拓展性

本分流域耦合方法可以很容易由二維軸對稱模型拓展到二維全流域模型甚至三維模型。甚至可以實現三維流域與二維流域之間的跨維度耦合。以三維模型為例，UDF代碼部分幾乎不需要進行修改。圖1中的Interior剛性面在三維模型中替換為圓面，Interface交界面則為圓柱面。但是為了保證動網格動態(tài)網格層變方法順利更新網格，需要對運動區(qū)域的網格沿軸線方向劃分規(guī)整的分層網格。

4 仿真結果及分析

為保證計算結果的準確性，對模型進行了網格無關性研究，主要考慮控制棒間隙流域網格層數(r方向)以及流體域軸線方向(z方向)網格層數對計算的影響。如表1所示，為控制棒(低位)在導向管流速2 m/s條件下計算得到的穩(wěn)態(tài)流體阻力值。從表中可以看出，控制棒流場計算對徑向網格數和軸向網格數均比較敏感，隨著徑向網格層數和軸向網格層數的增加，計算誤差在降低。綜合考慮計算精度和穩(wěn)定性，本文最終選取間隙流域網格層數(r方向)5層，軸線方向(z方向)網格層數6 000層作為網格劃分方案。

表1 網格敏感性分析

一般而言，當驅動桿和控制棒從一定高度由靜止狀態(tài)在重力作用下沿著驅動線流道下落，初始階段，由于流體阻力較小，落棒速度不斷增加。此時驅動桿上端由于耐壓殼內由于液體體積增加，造成了負壓，并且導致流體沿著驅動桿和耐壓殼之間的環(huán)形間隙向上流動，彌補驅動桿上端移動所形成的空腔。隨著驅動桿速度的不斷增加，驅動桿上端負壓和環(huán)形間隙中的流體速度也不斷增加，導致驅動桿流體阻力不斷變大。這與圖3(b)中0～0.7 s的變化趨勢一致。

圖3 流體阻力Fig. 3 Fluid resistance

對于控制棒而言，控制棒開始下落后，導向管中一部分流體通過導向管側壁流水孔和底部端塞排水孔向外排出，另外一部分流體通過控制棒與導向管間的環(huán)形間隙向上流動，形成間隙射流；同時導向管內壓強隨著控制棒速度的不斷增加而不斷增大，這導致控制棒的流體阻力也在不斷增加，但單根控制棒的流體阻力相對較小。當控制棒最終運動到緩沖段后，導向管側壁流水孔被迅速遮蔽，排水作用減弱。控制棒與導向管間的環(huán)形間隙由于緩沖段縮口的存在瞬間減小。同時流通截面積突變和驅動棒對導向管底部流體的迅速擠壓導致導向管內壓強急劇增加，在控制棒底部形成很大的壓差阻力，如圖3(a)所示t=0.8 s左右形成的水力阻力脈沖所示。在此壓差阻力的作用下，控制棒組件速度迅速減小。隨著控制棒運動速度的減小，驅動桿和導向管內壓差也相應不斷減小，因此流體阻力也減小，并在低棒位附近最終趨于平緩。

如圖4所示，為控制棒驅動線落棒速度和位移變化曲線，可以看到與控制棒落棒一般運動規(guī)律相一致，在進入緩沖段之前，控制棒在重力作用下不斷加速，但是由于速度增加，流體阻力變大，因此速度增加的趨勢逐漸變緩。當控制棒進入緩沖段之后，由于控制棒底部形成的壓差阻力及緩沖段環(huán)形小間隙帶來的水力摩擦阻力急劇增加，使得控制棒運動速度呈現斷崖式減小的趨勢，最終下降到一定速度后，穩(wěn)定下移直到最終位置。但是實驗得到的落棒速度比計算值偏小，這可能是由于本文模型假設控制棒沿著軸線方向在對中條件下落，并未考慮實際落棒過程中由于錯對中等因素帶來的偏心、機械摩擦等，導致計算得到的落棒阻力偏小，因此落棒速度也就偏大。

圖4 控制棒驅動線落棒速度和位移變化曲線Fig. 4 Variation curve of falling speed anddisplacement of the control rod drive line

圖5 控制棒前端速度場Fig. 5 Velocity field at the front end of the control rod

圖6 控制棒前端壓力場Fig. 6 Pressure field at the front of the control rod

圖5展示了控制棒下落過程中，控制棒前端速度場分布情況，圖6則展示了控制棒前端壓力場分布情況?？梢钥吹诫S著控制棒的快速下插，控制棒前端形成較高的壓力，導向管中的流體通過控制棒與導向管間的環(huán)形間隙向上流動，形成速度較大的間隙射流，這與前文描述的變化趨勢相一致。

圖7 壓力分布曲線Fig. 7 Pressure distribution curve

圖7給出了控制棒和驅動桿運動到某一位置時的表面壓力，橫坐標代表控制棒和驅動桿長度方向，其中0 m處表示棒和桿的上端，可以看到控制棒和驅動桿的上端均形成了負壓，而在前段則形成正壓，這與前文描述的規(guī)律一致。其中，驅動桿上端由于與耐壓殼形成相對密閉的空間，因此負壓作用明顯，控制棒則由于處于開闊空間，因此負壓作用較弱。

5 結論

本文提出了一種基于動網格的反應堆控制棒落棒行為分流域耦合仿真方法，該方法分別建立了控制棒單棒和驅動桿的二維軸對稱模型，保證控制棒和驅動桿對應流域的網格能夠根據落棒運動規(guī)律自適應地變化。兩個流域在每個時間步長內交換流體阻力計算結果，并根據運動方程求解得到的速度來更新控制棒和驅動桿的運動狀態(tài)，實現分流域耦合。該方法通用性好，計算中考慮了驅動線流道形狀的影響，且在計算時間和求解精度之間取得了良好的折中，此外，該方法在計算條件允許的情況下，還能較容易地擴展到三維模型。本文采用該方法對某反應堆驅動線落棒行為進行了仿真，仿真結果表明，無論是速度、水力阻力、位移隨時間變化，以及速度場、壓力場分布情況等，均與一般理解的控制棒驅動線落棒規(guī)律一致，說明該方法計算結果可信。