(中國核動(dòng)力研究設(shè)計(jì)院，四川 成都 610213)

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

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

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

1 流體仿真模型

1.1 網(wǎng)格模型

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

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

圖1 網(wǎng)格模型Fig.1 Mesh model

1.2 計(jì)算模型

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

固體壁面均設(shè)置為固壁邊界，其邊界設(shè)置為無滑移條件，近壁面應(yīng)用標(biāo)準(zhǔn)壁面函數(shù)。

壓力插值格式選擇標(biāo)準(zhǔn)格式，動(dòng)量選擇二階迎風(fēng)格式，湍動(dòng)能及湍流耗散率皆選取一階迎風(fēng)格式。

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

2 分流域耦合方法

由于網(wǎng)格量和軸對稱條件的限制，因此將控制棒驅(qū)動(dòng)線的流場分為控制棒單棒和驅(qū)動(dòng)桿兩個(gè)部分。兩個(gè)流場迭代計(jì)算同時(shí)進(jìn)行，在每個(gè)時(shí)間步內(nèi)進(jìn)行數(shù)據(jù)交換以保持同步。

基本計(jì)算流程如圖2所示，計(jì)算過程如下：

1)控制棒單棒流體域計(jì)算為主計(jì)算流程，在每個(gè)時(shí)間步長內(nèi)首先計(jì)算當(dāng)前時(shí)間步長Δt，通過DEFINE_DELTAT宏修改時(shí)間步長，并輸出Δt到文件time.txt，提供給驅(qū)動(dòng)桿流域計(jì)算，驅(qū)動(dòng)桿流域讀入time.txt后同樣通過DEFINE_DELTAT宏修改時(shí)間步長；

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

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

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

(1)

式中：M——控制棒和驅(qū)動(dòng)桿的總重;

g——重力加速度。

同時(shí)更新速度，并寫入文件vel.txt，根據(jù)更新速度值利用DEFINE_GRID_MOTION宏更新動(dòng)網(wǎng)格區(qū)域網(wǎng)格。

v=v+Δv

(2)

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

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

6)數(shù)據(jù)交互完成并更新網(wǎng)格之后，兩個(gè)流域同時(shí)開始迭代計(jì)算，迭代完成后進(jìn)入下一時(shí)間步直到計(jì)算結(jié)束。

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

3 方法拓展性

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

4 仿真結(jié)果及分析

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

表1 網(wǎng)格敏感性分析

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

圖3 流體阻力Fig. 3 Fluid resistance

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

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

圖4 控制棒驅(qū)動(dòng)線落棒速度和位移變化曲線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則展示了控制棒前端壓力場分布情況。可以看到隨著控制棒的快速下插，控制棒前端形成較高的壓力，導(dǎo)向管中的流體通過控制棒與導(dǎo)向管間的環(huán)形間隙向上流動(dòng)，形成速度較大的間隙射流，這與前文描述的變化趨勢相一致。

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

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

5 結(jié)論

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