一種基于動(dòng)網(wǎng)格的反應(yīng)堆控制棒落棒行為分流域耦合仿真方法

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

反應(yīng)堆控制棒驅(qū)動(dòng)線作為反應(yīng)堆內(nèi)具有相對(duì)運(yùn)動(dòng)的設(shè)備單元。一般由燃料組件、導(dǎo)向組件、控制棒組件、控制棒驅(qū)動(dòng)機(jī)構(gòu)等部分組成。由于控制棒組件對(duì)堆芯反應(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)線簡(jiǎn)化為一維水力模型，取得了大量的成果。這些成果作為控制棒驅(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ì)算，因此對(duì)實(shí)驗(yàn)依賴性較高。隨著CFD計(jì)算水平的提高和計(jì)算機(jī)技術(shù)的進(jìn)步，一些學(xué)者選擇采用CFD動(dòng)網(wǎng)格方法對(duì)控制棒落棒行為進(jìn)行仿真分析，例如，肖聰?shù)然?FLUENT 動(dòng)網(wǎng)格技術(shù)對(duì)單根圓柱形控制棒建立了三維流體仿真模型來進(jìn)行落棒行為仿真分析。同樣，基于流體動(dòng)力學(xué)(CFD)的動(dòng)網(wǎng)格技術(shù)，在新的堆型上模擬了十字形控制棒組件在控制棒導(dǎo)向組件內(nèi)的運(yùn)動(dòng)行為[8-9]。但是由于三維模型計(jì)算量過大，該方法只對(duì)驅(qū)動(dòng)線的局部流域進(jìn)行了仿真建模分析。到目前為止還沒有針對(duì)完整驅(qū)動(dòng)線的有限元計(jì)算分析。

本文提出了一種基于動(dòng)網(wǎng)格的反應(yīng)堆控制棒落棒行為分流域耦合仿真方法，該方法分別建立了控制棒單棒和驅(qū)動(dòng)桿的二維軸對(duì)稱模型，保證控制棒和驅(qū)動(dòng)桿對(duì)應(yīng)流域的網(wǎng)格能夠根據(jù)落棒運(yùn)動(dòng)規(guī)律自適應(yīng)地變化。兩個(gè)流域在每個(gè)時(shí)間步長(zhǎng)內(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)線長(zhǎng)度較長(zhǎng)，而控制棒和驅(qū)動(dòng)桿周圍的流體間隙相對(duì)小很多，要獲得與軸向長(zhǎng)度和徑向間隙長(zhǎng)度相適應(yīng)的網(wǎng)格尺寸，網(wǎng)格量將相當(dāng)大。因此，為了減少網(wǎng)格模型大小，本文選用二維軸對(duì)稱模型來對(duì)驅(qū)動(dòng)線流場(chǎng)進(jìn)行模擬，且為了保證計(jì)算精度和動(dòng)網(wǎng)格更新，整個(gè)流域采用四邊形網(wǎng)格進(jìn)行網(wǎng)格劃分。

驅(qū)動(dòng)線流場(chǎng)分為控制棒流域和驅(qū)動(dòng)桿流域兩個(gè)部分，忽略星形架等結(jié)構(gòu)，并對(duì)流體域流道作了適當(dāng)簡(jiǎn)化。以圖1所示的控制棒流域?yàn)槔刂瓢襞c導(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)而不會(huì)產(chǎn)生負(fù)網(wǎng)格。

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

1.2 計(jì)算模型

本文采用FLUENT 14.5版本進(jìn)行仿真計(jì)算。由于采用二維軸對(duì)稱模型，因此計(jì)算時(shí)求解的是二維軸對(duì)稱的連續(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)格量和軸對(duì)稱條件的限制，因此將控制棒驅(qū)動(dòng)線的流場(chǎng)分為控制棒單棒和驅(qū)動(dòng)桿兩個(gè)部分。兩個(gè)流場(chǎng)迭代計(jì)算同時(shí)進(jìn)行，在每個(gè)時(shí)間步內(nèi)進(jìn)行數(shù)據(jù)交換以保持同步。

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

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

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 方法拓展性

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

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

為保證計(jì)算結(jié)果的準(zhǔn)確性，對(duì)模型進(jìn)行了網(wǎng)格無關(guān)性研究，主要考慮控制棒間隙流域網(wǎng)格層數(shù)(r方向)以及流體域軸線方向(z方向)網(wǎng)格層數(shù)對(duì)計(jì)算的影響。如表1所示，為控制棒(低位)在導(dǎo)向管流速2 m/s條件下計(jì)算得到的穩(wěn)態(tài)流體阻力值。從表中可以看出，控制棒流場(chǎng)計(jì)算對(duì)徑向網(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的變化趨勢(shì)一致。

圖3 流體阻力Fig. 3 Fluid resistance

對(duì)于控制棒而言，控制棒開始下落后，導(dǎo)向管中一部分流體通過導(dǎo)向管側(cè)壁流水孔和底部端塞排水孔向外排出，另外一部分流體通過控制棒與導(dǎo)向管間的環(huán)形間隙向上流動(dòng)，形成間隙射流；同時(shí)導(dǎo)向管內(nèi)壓強(qiáng)隨著控制棒速度的不斷增加而不斷增大，這導(dǎo)致控制棒的流體阻力也在不斷增加，但單根控制棒的流體阻力相對(duì)較小。當(dāng)控制棒最終運(yùn)動(dòng)到緩沖段后，導(dǎo)向管側(cè)壁流水孔被迅速遮蔽，排水作用減弱。控制棒與導(dǎo)向管間的環(huán)形間隙由于緩沖段縮口的存在瞬間減小。同時(shí)流通截面積突變和驅(qū)動(dòng)棒對(duì)導(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)入緩沖段之前，控制棒在重力作用下不斷加速，但是由于速度增加，流體阻力變大，因此速度增加的趨勢(shì)逐漸變緩。當(dāng)控制棒進(jìn)入緩沖段之后，由于控制棒底部形成的壓差阻力及緩沖段環(huán)形小間隙帶來的水力摩擦阻力急劇增加，使得控制棒運(yùn)動(dòng)速度呈現(xiàn)斷崖式減小的趨勢(shì)，最終下降到一定速度后，穩(wěn)定下移直到最終位置。但是實(shí)驗(yàn)得到的落棒速度比計(jì)算值偏小，這可能是由于本文模型假設(shè)控制棒沿著軸線方向在對(duì)中條件下落，并未考慮實(shí)際落棒過程中由于錯(cuò)對(duì)中等因素帶來的偏心、機(jī)械摩擦等，導(dǎo)致計(jì)算得到的落棒阻力偏小，因此落棒速度也就偏大。

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

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

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

圖5展示了控制棒下落過程中，控制棒前端速度場(chǎng)分布情況，圖6則展示了控制棒前端壓力場(chǎng)分布情況。可以看到隨著控制棒的快速下插，控制棒前端形成較高的壓力，導(dǎo)向管中的流體通過控制棒與導(dǎo)向管間的環(huán)形間隙向上流動(dòng)，形成速度較大的間隙射流，這與前文描述的變化趨勢(shì)相一致。

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

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

5 結(jié)論

本文提出了一種基于動(dòng)網(wǎng)格的反應(yīng)堆控制棒落棒行為分流域耦合仿真方法，該方法分別建立了控制棒單棒和驅(qū)動(dòng)桿的二維軸對(duì)稱模型，保證控制棒和驅(qū)動(dòng)桿對(duì)應(yīng)流域的網(wǎng)格能夠根據(jù)落棒運(yùn)動(dòng)規(guī)律自適應(yīng)地變化。兩個(gè)流域在每個(gè)時(shí)間步長(zhǎng)內(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ò)展到三維模型。本文采用該方法對(duì)某反應(yīng)堆驅(qū)動(dòng)線落棒行為進(jìn)行了仿真，仿真結(jié)果表明，無論是速度、水力阻力、位移隨時(shí)間變化，以及速度場(chǎng)、壓力場(chǎng)分布情況等，均與一般理解的控制棒驅(qū)動(dòng)線落棒規(guī)律一致，說明該方法計(jì)算結(jié)果可信。