控制棒落棒時間分析軟件CRAC的研發(fā)與驗證

馬 超，周 毅，盧嘉川，蘆 韡，李 云，趙 波，陳 浩，馬 強，秦 勉，郭曉明

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

壓水堆緊急停堆時，控制棒組件依靠重力下落控制反應性，實現(xiàn)安全停堆。落棒時間計算是反應堆安全分析和燃料組件及相關組件設計的重要部分，關系到導向管及控制棒設計是否滿足落棒時間準則。針對落棒開展理論研究，研發(fā)落棒時間的分析工具，獲得落棒位移-速度-時間關系，對反應堆安全分析及燃料組件設計具有重要意義。

針對壓水堆落棒時間分析技術，國內(nèi)外研究者開展大量研究工作，通過試驗與理論分析，認為落棒過程中控制棒驅(qū)動線受到的載荷主要為水力阻力和機械阻力。針對機械阻力，研究者利用有限差分或有限元法，求解驅(qū)動線各部件的振動方程，通過設置表面接觸條件，計算控制棒受到的非線性碰撞力與摩擦力[1-2]。水力阻力具體包括浮力、慣性力、局部形阻力和沿程阻力。針對水力阻力，研究者利用流體力學理論和經(jīng)驗公式相結合的方法計算落棒過程中流體阻力[3]，還有研究者利用CFD軟件動網(wǎng)格技術，開展落棒過程中水力緩沖模擬[4-7]。上述研究工作考慮的落棒模型較為簡單，僅圍繞某一種載荷建模，如控制棒所受水力阻力或機械阻力，忽略或簡化考慮其他載荷因素，距離工程應用較遠。動網(wǎng)格技術需要耗費大量計算資源，只能構建簡單控制棒結構的仿真模型，如單根控制棒和水力緩沖部件，無法開展全尺寸控制棒驅(qū)動線模擬。

國內(nèi)外核電企業(yè)采用理論和經(jīng)驗關系式相結合的方式求解落棒過程中的水力阻力和機械阻力[8-10]，在試驗數(shù)據(jù)基礎上，形成針對特定堆型的落棒專用程序，開展落棒時間計算。本文介紹適用于華龍一號工程設計的落棒時間分析程序CRAC，滿足華龍一號落棒時間分析需求。

1 CRAC功能

CRAC通過合理簡化控制棒驅(qū)動線結構與載荷形式，建立動力學方程，基于流體力學壓降平衡關系、彈簧胡克定律及經(jīng)驗機械阻力關系，求解落棒過程中的流體阻力與機械阻力，獲取反應堆正常工況及地震停堆工況，落棒總體時間及隨時間變化的控制棒位移、速度信息，適用于華龍一號控制棒驅(qū)動線的落棒時間計算。

2 理論模型

2.1 總體理論模型

控制棒驅(qū)動線結構如圖1所示，導向通道在軸向由下至上分為3部分：堆芯燃料組件部分、導向筒部件部分和控制棒驅(qū)動機構部分。反應堆正常運行狀態(tài)下，上部的驅(qū)動機構鉤爪提升驅(qū)動桿，將控制棒組件大部分停留在堆芯外部，控制棒棒束保持在導向筒部件中；需要停堆時，鉤爪部件斷電，松開驅(qū)動桿，控制棒驅(qū)動線依靠重力下落到燃料組件中，實現(xiàn)停堆。驅(qū)動桿、控制棒為細長結構(長徑比大于100∶1)，控制棒驅(qū)動線與通道之間間隙較小，限制了徑向運動，因此在分析落棒過程中，近似為一維軸向運動。

圖1 反應堆控制棒驅(qū)動線結構Fig.1 Structure of reactor control rod drive line

圖2示出控制棒驅(qū)動線受力示意圖，針對控制棒軸向運動建模，所有作用力都作用在重心，建立運動方程如下：

圖2 控制棒驅(qū)動線受力示意圖Fig.2 Force analysis schematic of control rod drive line

(1)

2.2 流體阻力計算模型

控制棒在驅(qū)動機構、導向筒及燃料組件3部分通道中運動。由于在導向筒中流場較為寬闊，水力阻力較小，故忽略導向筒內(nèi)作用在行星架與控制棒上的水力阻力。運動部件所受流體阻力，主要考慮控制棒在燃料組件中及驅(qū)動桿在驅(qū)動機構中受到的流體阻力。引入如下假設：1) 流體不可壓縮；2) 落棒過程中在燃料組件導向管和緩沖段內(nèi)的水溫為常數(shù)；3) 在每個計算時間步內(nèi)，固定壁面摩擦系數(shù)與運動壁面摩擦系數(shù)為定值(兩者選取不同參考速度)。

2.2.1控制棒在燃料組件內(nèi)所受到的水力阻力 圖3示出控制棒插入到導向管和緩沖段示意圖，采用集總參數(shù)法描述流道速度信息，如表1所列。由于堆芯流量可能在較大范圍內(nèi)變化及控制棒組件機械緩沖回彈現(xiàn)象存在，各流道速度方向可能發(fā)生變化。在水力阻力建模前，需定義速度標志位，保證對流體動力學方程求解的正確。

圖3 控制棒插入到導向管和緩沖段示意圖Fig.3 Schematic of control rod dropin thimble guide and dashpot

表1 各通道流體速度定義Table 1 Definition of fluid velocity direction in different flow channels

1) 控制棒在導向管內(nèi)運動

在燃料組件中所受水力阻力，由壓差引起的作用力和沿程阻力組成：

FFA=NCR(SCR(pB-p2)+lABCRτCR)

(2)

式中：FFA為控制棒組件在燃料組件內(nèi)受到的水力阻力；NCR為控制棒組件中棒總數(shù)；SCR為單根控制棒橫截面積；pB為控制棒下端界面處壓強；p2為導向管頂端外部壓強；lA為控制棒插入導向管軸向長度；BCR為單根控制棒濕潤周長；τCR為控制棒在導向管內(nèi)受到的平均流體切應力。

連續(xù)性方程為：

vSCR-vASA-4vDHSDH=0

(3)

4vDHSDH=vTGSTG

(4)

式中：STG為導向管圓形通道過流面積；SA為導向管環(huán)形通道過流面積；SDH為單個流水孔過流面積(周向分布共4個)。流水孔外側和控制棒下端之間的流體壓差為：

(5)

式中：p1為流水孔外側的流體壓強；ξDH為流水孔局部阻力系數(shù)；λTG為導向管圓形通道的沿程阻力系數(shù)；DTG為導向管通道水力直徑；lTG為導向管總長度。

在導向管環(huán)形通道內(nèi)，控制棒下端與導向管上端之間壓差表示如下：

(6)

λ按Colebrook公式計算[11]：

(7)

式中：Ks為壁面絕對粗糙度；Re為雷諾數(shù)，Re=vrefDρ/η，vref為參考速度，η為流體黏度；D為水力直徑。

從流水孔外側到導向管頂端壓降Δp12=p1-p2，由燃料組件特性、進入燃料組件流量及流體溫度共同確定，因此Δp12視為已知量，根據(jù)壓降平衡關系：

Δp12=(p1-pB)+(pB-p2)

(8)

在確定每一時間步v(t)條件下，聯(lián)立式(3)～(6)及(8)可確定vA(t)，從而確定當前時刻控制棒壓差(pB-p2)和切向應力分量：

(9)

根據(jù)式(2)獲得控制棒受到的水力阻力FFA。

2) 控制棒在緩沖段內(nèi)運動

當控制棒插入緩沖段時，控制棒組件受到的水力阻力FFA由式(10)計算，此時導向管內(nèi)環(huán)形流道長度為lA=lTG。

(10)

式中：lB為控制棒插入緩沖段內(nèi)長度；CCRB為緩沖段環(huán)形流道中控制棒壁面摩擦系數(shù)。此時，連續(xù)性方程如下：

4vDHSDH+vASA-vBSB=0

(11)

vSCR+vDPSDP+vBSB=0

(12)

式中：SB為緩沖段環(huán)形流道過流面積；SDP為緩沖段圓形流道過流面積?？刂瓢粝露伺c導向管頂端之間壓降分為兩段：

(13)

(14)

式中：pDH為流水孔內(nèi)側流體壓強；CDP為環(huán)形流道緩沖段壁面?zhèn)饶Σ料禂?shù)；ξ1BS為環(huán)形流道截面突然擴張的局部壓力損失系數(shù)，即從SB突擴至SA的壓力損失系數(shù)。由導向管緩沖段對應的外部流道壓降平衡關系為：

pB-pDH=(pB-p0)+

(p0-p1)+(p1-pDH)

(15)

式中，p0為軸肩螺釘阻尼孔底端壓力。

(16)

(17)

式中：lDP為緩沖段總長；λDP為緩沖段圓形通道沿程損失系數(shù)；λSH為軸肩螺釘阻尼孔中沿程損失系數(shù)；DDP為緩沖段圓形通道水力直徑；DSH為軸肩螺釘孔水力直徑；lSH為軸肩螺釘孔深；ξSH為軸肩螺釘阻尼孔局部壓力損失系數(shù)；SSH為軸肩螺釘孔過流面積。因此式(15)可轉(zhuǎn)化為：

(18)

式中，Δp01=p0-p1視為已知。

在已知v(t)情況下，聯(lián)立式(13)、(14)、(16)和(18)可確定vA(t)、vB(t)，進而確定壓差(pTG-pA2)，通過式(10)計算控制棒插入緩沖段時受到的水力阻力FFA。

2.2.2驅(qū)動桿在驅(qū)動機構內(nèi)受到的水力阻力 驅(qū)動桿在驅(qū)動機構中受到的水力阻力建模過程與2.2.1節(jié)內(nèi)容相似，通過求解連續(xù)性方程和動力學方程，可得驅(qū)動桿受到的水力阻力FCRDM。

2.3 機械阻力計算模型

在控制棒驅(qū)動線落棒過程中，驅(qū)動桿與控制棒會發(fā)生變形，與導向通道發(fā)生隨機碰撞接觸，機械運動行為復雜，直接針對驅(qū)動線相關設備結構形式建模，計算模型將非常龐大且求解困難。因此在CRAC研發(fā)過程中，考慮基于控制棒驅(qū)動線相關試驗結果與工程經(jīng)驗，建立機械阻力相對于錯對中量及插入深度變化的計算模型：

Ff=F(d,z)+F0

(19)

式中：函數(shù)F(d，z)輸入變量包括上部導向組件與隔熱套之間錯對中量d及落棒位移z，F(xiàn)(d，z)是在參考試驗數(shù)據(jù)條件下，考慮工程保守性建立的插值函數(shù)；F0為對中條件下的機械摩擦力。

3 CRAC計算流程

圖4示出CRAC計算流程，通過輸入整個系統(tǒng)初始參數(shù)，包括結構尺寸參數(shù)、冷卻劑物性參數(shù)、控制棒初始位置與速度、各通道中流體初始速度和壓強，確定控制棒驅(qū)動線初始加速度。將時間步長設為Δt，第i時間步對應時刻ti=ti-1+Δt，該時刻所有載荷與流體狀態(tài)將根據(jù)時刻ti-1的特性確定。每一時間步求出控制棒驅(qū)動線所受合力，求得加速度及其他下落特性參數(shù)。重復上述過程，直至控制棒組件完全插入燃料組件，實現(xiàn)對落棒過程的模擬計算。

圖4 CRAC計算流程圖Fig.4 Calculation flow chart of CRAC

4 CRAC驗證

4.1 落棒試驗簡介

CF2燃料組件是由中國核動力研究設計院自主設計研發(fā)，是第三代核電華龍一號核心部件。為滿足第三代核電燃料組件剛度要求，CF2在導向管緩沖結構方面采用管中管新型設計。該燃料組件將應用于華龍一號國外首堆巴基斯坦卡拉奇2號機組的首循環(huán)。本次落棒試驗中回路參數(shù)為：壓力15.5 MPa，溫度315 ℃，工況列于表2[12]。試驗測量得到：控制棒下落至緩沖段入口的落棒時間(T5)；控制棒下落至全插入位置的落棒時間(T5+T6)；控制棒下落速度、位移與時間關系。落棒時間測量誤差為±0.01 ms，落棒速度測量系統(tǒng)相對誤差為±0.28%。

4.2 計算結果與試驗值對比

根據(jù)試驗條件及控制棒驅(qū)動線結構參數(shù)，利用CRAC進行表2所列5種工況下的落棒模擬計算。圖5示出G1和G5工況下速度-時間和位移-時間曲線。由圖5可知：CRAC計算結果與試驗結果變化趨勢一致，但計算最大速度值要比試驗值略低；CRAC計算的位移曲線趨勢較為平緩，落棒時間結果要比試驗值大。因為CRAC中機械摩擦阻力及水力局部壓降系數(shù)均為經(jīng)驗函數(shù)或經(jīng)驗系數(shù)，充分考慮落棒時間計算保守性，保證軟件計算值大于真實值。

圖5 工況G1(a)、G5(b)條件下的落棒結果對比Fig.5 Comparison of control rod drop result for G1 (a), G5 (b) conditions

表2 落棒試驗工況Table 2 Condition of drop experiment

落棒試驗獲得的速度曲線在控制棒進入緩沖段之前達到最大，此時曲線出現(xiàn)明顯的突然上升、然后下降過程。產(chǎn)生該現(xiàn)象的可能原因是，控制棒下落的過程中，子彈頭部位會產(chǎn)生局部過壓，在導向管中運動過程中，距離流水孔較遠時，過壓產(chǎn)生的水力阻力壓差增大較為平緩連續(xù)，而在控制棒子彈頭接近流水孔時，流水孔處內(nèi)外壓差較大，流水孔起到泄壓作用，使得控制棒子彈頭部位過壓在流水孔附近產(chǎn)生先降低(泄壓)再升高(進入緩沖段)的變化，因此速度曲線在此處先迅速升高再迅速下降。由于CRAC采用一維集總參數(shù)模型，將控制棒簡化為圓柱體，無法考慮子彈頭部位的壓力場與速度場急劇變化現(xiàn)象，因此獲得的速度-時間曲線是逐漸變化，沒有突然上升再下降的過程。

表3列出試驗數(shù)據(jù)(T5，T5+T6)與計算結果對比，計算值均大于試驗值，實現(xiàn)對試驗數(shù)據(jù)的包絡，滿足反應堆安全分析保守性要求。落棒時間計算值同試驗值符合性較好，最大相對誤差為15.45%，證明CRAC對華龍一號堆型落棒時間計算的適用性以及準確性。

表3 落棒時間計算與試驗數(shù)據(jù)對比Table 3 Comparison of calculated and measured values for control rod drop time

5 結論

本文介紹了CRAC的功能、理論模型、計算流程，并利用試驗數(shù)據(jù)對CRAC的計算準確性進行驗證，計算結果與試驗數(shù)據(jù)相對誤差在15.45%以內(nèi)，兩者符合較好，且計算結果均大于試驗結果，表明軟件計算精度與保守性能滿足華龍一號堆型安全停堆時間準則分析需求。