莫錦濤

（中國核動力研究設計院核反應堆系統(tǒng)設計技術(shù)重點實驗室 成都 610213）

控制棒驅(qū)動線的落棒行為是指該設備在反應堆的壽命周期內(nèi)，在高輻照、高溫、高壓并伴隨有流致振動的條件下，控制棒組件（這里還包括驅(qū)動桿、可拆接頭等部件）在驅(qū)動機構(gòu)得到釋放指令之后，在反應堆內(nèi)的運動情況。由于控制棒驅(qū)動線落棒行為關(guān)系到事故工況下的安全停堆。因此落棒時間是反應堆結(jié)構(gòu)設計時首要考慮的參數(shù)之一。

目前反應堆控制棒驅(qū)動線的落棒時間主要靠驅(qū)動線的相關(guān)綜合試驗來進行驗證，但與此同時，控制棒驅(qū)動線的落棒時間的仿真研究也在開展。目前主要有兩種方法用于反應堆控制棒驅(qū)動線的落棒時間的計算，一種方法是將控制棒驅(qū)動線簡化為一維水力模型，對這種方法而言，需要求解的落棒動力學方程中無法通過理論計算得到水力阻力和機械摩擦阻力系數(shù)，這些數(shù)據(jù)只能依賴相關(guān)的試驗獲得，且無法考慮流道形狀變化對仿真結(jié)果的影響。Donis等［1］首次通過求解一維動力學方程來計算得到控制棒落棒時間，在此基礎上建立了數(shù)學模型用于分析單根控制棒落棒行為。美國、法國等為了對核電站壓水堆控制棒落棒時間進行計算，均研制開發(fā)了相應的專用軟件［2-4］，但由于這些軟件均采用了一維動力學模型，引入了大量與試驗相關(guān)的修正系數(shù)，因此在對不同堆型進行計算時需要對程序進行大量修改，拓展性較差。

另外一種方法是計算流體力學（Computational Fluid Dynamics，CFD）動網(wǎng)格計算方法。例如，肖聰?shù)龋?-6］基于計算流體力學動網(wǎng)格技術(shù)，對某反應堆的控制棒驅(qū)動線導向組件以及單根控制棒建立了相應的三維流體仿真模型，并利用FLUENT動網(wǎng)格技術(shù)進行瞬態(tài)動力學計算。但是由于三維模型計算量過大，因此該方法只對驅(qū)動線的局部流域進行了仿真建模分析，且由于計算量過大，計算時間過長，這樣的方法距離實際工程應用還很遙遠。到目前為止還沒有針對完整驅(qū)動線的有限元計算分析。

隨著CFD技術(shù)手段的進步，近年來，重疊網(wǎng)格技術(shù)取得了長足發(fā)展。動網(wǎng)格和滑移網(wǎng)格已經(jīng)不再是模擬運動部件在流體中相對運動的唯一選擇。重疊網(wǎng)格主要用來處理流體域中的部件運動，該技術(shù)的核心在于采用了背景網(wǎng)格和部件網(wǎng)格。部件網(wǎng)格用于模擬運動部件，通過部件網(wǎng)格和背景網(wǎng)格在不同相對位置處相互重合區(qū)域的網(wǎng)格進行數(shù)據(jù)交換，重疊網(wǎng)格較容易地解決了復雜相對運動下的流固耦合計算問題［7-15］。

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

1 流體仿真模型

1.1 網(wǎng)格模型

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

驅(qū)動線流場分為控制棒流域和驅(qū)動桿流域兩個部分，忽略星形架等結(jié)構(gòu)，并對流體域流道作了適當簡化。以圖1所示的控制棒流域為例，導向管內(nèi)的流域被劃分為背景網(wǎng)格，而圍繞控制棒的部分流域被劃分為部件網(wǎng)格，部件網(wǎng)格外邊界定義為Overset邊界。背景網(wǎng)格在對應控制棒部件網(wǎng)格徑向間隙的部位進行了局部加密（本例中設置了10層網(wǎng)格），需要說明的是，控制棒部件網(wǎng)格徑向間隙尺寸是由控制棒與導向管之間的最小間隙決定的（通常為緩沖段），而背景網(wǎng)格需要在控制棒部件網(wǎng)格徑向間隙加密區(qū)掃過的區(qū)域進行相應的網(wǎng)格加密（如圖1中10層網(wǎng)格深色區(qū)域所示），保證背景網(wǎng)格和部件網(wǎng)格的重疊區(qū)域至少有5層以上的網(wǎng)格。

圖1 分流域耦合網(wǎng)格模型Fig.1 Multi-field coupling mesh model

1.2 計算模型

由于采用二維軸對稱模型，因此計算時求解的是二維軸對稱形式下的動量方程和連續(xù)性方程，湍流模型采用k-ε湍流模型。

本文計算條件選取反應堆冷態(tài)工況，常溫常壓，并假設控制棒從最高棒位由靜止狀態(tài)開始下落。計算流體域所有出口均設置為全壓邊界，流體域出口包括緩沖段流水孔、底部端塞排水孔等。固體壁面采用了無滑移無滲透的壁面條件并且采用壁面函數(shù)模擬邊界層。模型求解器選用Coupled算法進行求解。為了兼顧計算穩(wěn)定性和計算求解速度，除了壓力和動量方程采用二階離散格式，其他方程均選取一階迎風格式。

2 分流域耦合方法

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

基本的耦合計算流程如圖2所示，完整的計算過程如下：

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

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

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

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

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

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

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

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

為保證計算結(jié)果的準確性，對模型進行了網(wǎng)格無關(guān)性研究，主要考慮控制棒部件網(wǎng)格和背景網(wǎng)格的重疊區(qū)域網(wǎng)格層數(shù)（r方向）以及流體域軸線方向（z方向）網(wǎng)格層數(shù)對計算的影響。如表1所示，為控制棒（低位）在導向管流速2 m·s-1條件下計算得到的穩(wěn)態(tài)流體阻力值。表1中選取不同的徑向（r）和軸向（z）網(wǎng)格數(shù)時的計算結(jié)果。從表1中可以看出，控制棒流場計算對徑向網(wǎng)格數(shù)和軸向網(wǎng)格數(shù)均比較敏感，隨著徑向網(wǎng)格層數(shù)和軸向網(wǎng)格層數(shù)的增加，流體阻力值呈增大趨勢，但計算誤差的減小幅度在降低。綜合考慮計算精度和穩(wěn)定性，本文最終選取重疊區(qū)域網(wǎng)格層數(shù)（r方向）10層，軸線方向（z方向）網(wǎng)格層數(shù)6 000層作為網(wǎng)格劃分方案。

一般而言，控制棒下落的過程中受到的力包括機械摩擦力、浮力、重力、流體阻力，在這些力的共同作用下完成落棒運動，為了簡化，假設控制棒在對中條件下落棒，因此機械摩擦力可以忽略。控制棒在初始落棒位置時，由于落棒速度小，流體阻力很小，因此控制棒在重力作用下加速運動。在控制棒插入導向管過程中，不斷對導向管中的流體進行擠壓。受擠壓的流體一部分被壓出底部端塞排水孔和導向管側(cè)壁流水孔，另外一部分則沿著控制棒與導向管之間的環(huán)形空間流出導向管?？刂瓢羲艿牧黧w阻力與控制棒端部的壓差阻力以及控制棒與導向管之間的間隙大小有關(guān)，當間隙較小時，控制棒與導線管之間的速度梯度較大，流體作用在控制棒上的剪切力也較大。因此當控制棒下落到緩沖段時，帶來流通截面積的瞬間變化，因此導向管內(nèi)產(chǎn)生巨大壓差。同時隨著導向管與控制棒間的環(huán)腔體積瞬間減小，流體剪切力迅速增大，控制棒阻力也相應變大。在這些阻力的共同作用下，控制棒運動速度迅速減小，最終趨于平穩(wěn)達到最終位置。

圖3為控制棒驅(qū)動線落棒速度和位移變化曲線，結(jié)合圖4中的控制棒和驅(qū)動桿流體阻力變化。可以看到與實驗得到的控制棒落棒運動規(guī)律相一致，在進入緩沖段之前，控制棒在重力作用下不斷加速，但是由于控制棒速度變快，相應的流體阻力變大，因此速度增加的趨勢逐漸變緩。當控制棒進入緩沖段之后，由于控制棒底部形成的壓差阻力以及緩沖段環(huán)形小間隙帶來的摩擦阻力急劇增加，使得控制棒運動速度呈現(xiàn)斷崖式減小的趨勢，最終下降到一定速度后，穩(wěn)定下移直到最終位置。從計算得到的控制棒落棒運動規(guī)律，可以看到與實際落棒運動特征非常符合，但是實驗得到的落棒速度比計算值偏小，這可能是由于本文模型假設控制棒沿著軸線方向在對中條件下落，并未考慮實際落棒過程中由于錯對中等因素帶來的偏心、機械摩擦等，導致計算得到的落棒阻力偏小，因此落棒速度也就偏大。圖5給出了控制棒進入緩沖段前后流水孔附近的壓力分布變化，從圖5中可以非常直觀地看到，在進入緩沖段之前，流水孔附近的流體由于控制棒落棒運動的擠壓而壓力增加，部分流體由于高壓而沿流水孔流出。而當控制棒進入緩沖段，流水孔附近的流體失去了運動擠壓而壓力迅速降低。

表1 網(wǎng)格敏感性分析Table 1 Mesh-density sensitivity analysis

圖3 控制棒驅(qū)動線落棒速度和位移變化曲線Fig.3 The velocity and displacement curves of the control rod driving line

圖4 控制棒驅(qū)動線落棒流體阻力變化曲線Fig.4 The fluid resistance curve of the control rod driving line

為了進一步驗證耦合方法的正確，同時也對緩沖段流水孔的作用做一些研究。圖6和圖7分別給出了在流水孔直徑減小一半、甚至不開孔的條件下，控制棒的落棒速度、流體阻力的對比?？梢钥吹?，在不開孔的條件下，由于缺少流水孔對導向管中受擠壓流體的泄流作用，導致控制棒所受的壓差阻力變大，控制棒下落速度整體變低（相比開孔條件下最高速度偏小0.7 m·s-1左右），進而使得控制棒落棒時間延長（相比開孔條件下到達緩沖段的時間變慢0.5～0.6 s）。此外，由于不開孔條件下，控制棒進入緩沖段的速度較小，因此緩沖段帶來的阻力增加幅度也較小。

驅(qū)動桿所受流體阻力主要與落棒速度有關(guān)，驅(qū)動桿下降過程中由于耐壓殼內(nèi)液體體積增加，造成了負壓，并且導致流體沿著驅(qū)動桿和耐壓殼之間的環(huán)形間隙向上流動，彌補驅(qū)動桿上端移動所形成的空腔。間隙的流體流動同樣給驅(qū)動桿帶來了流體阻力。因此，在不開孔條件下，驅(qū)動桿流體阻力由于整體落棒速度的降低而偏小。

從圖7對比還可以看到，落棒時間和流體阻力對孔徑變化的影響敏感程度較小，但是對是否開孔的影響非常敏感。因此，流水孔的設置對控制棒落棒來說不可或缺，但是其開孔大小則可以進行一定的優(yōu)化。

圖5 緩沖段流水孔壓力分布變化(a)控制棒進入緩沖段前，(b)控制棒進入緩沖段后Fig.5 The pressure distribution of the holes in buffer section of control rod drive line(a)Before the buffer,(b)After the buffer

圖6 控制棒驅(qū)動線落棒速度變化對比Fig.6 The comparison of velocity curves of the control rod driving line

4 結(jié)語

本文提出了一種基于重疊網(wǎng)格的反應堆控制棒落棒行為分流域耦合仿真方法，該方法分別建立了控制棒單棒和驅(qū)動桿的二維軸對稱重疊網(wǎng)格模型，保證控制棒和驅(qū)動桿的部件網(wǎng)格和背景網(wǎng)格能夠根據(jù)落棒運動規(guī)律自適應地變化。兩個流域在每個時間步長內(nèi)交換流體阻力計算結(jié)果，并根據(jù)運動方程求解得到的速度來更新控制棒和驅(qū)動桿的運動狀態(tài)，實現(xiàn)分流域耦合。本方法建模簡單，方法通用性好，計算中考慮了驅(qū)動線流道形狀的影響，且在計算時間和求解精度之間取得了良好的折中，此外，本方法在計算條件允許的情況下，還能較容易地擴展到三維模型。本文采用該方法模擬了某反應堆驅(qū)動線控制棒的落棒行為，仿真結(jié)果得到的速度、位移隨時間變化與實驗值符合良好，計算得到的速度場分布情況等與一般理解上的控制棒驅(qū)動線落棒規(guī)律特性一致，說明該方法計算結(jié)果可信。此外，還對不同緩沖段流水孔大小情況下控制棒的落棒速度、流體阻力變化情況進行了對比分析，發(fā)現(xiàn)落棒時間和流體阻力對孔徑變化的影響敏感程度較小，但是對是否開孔的影響非常敏感。流水孔的設置對控制棒落棒來說不可或缺，但是其開孔大小則可以進行一定的優(yōu)化。

圖7 控制棒驅(qū)動線落棒流體阻力變化對比Fig.7 The comparison of fluid resistance curves of the control rod driving line