劉 新 陳先龍 張 修 于曉雷

（深圳中廣核工程設(shè)計有限公司 深圳 518172）

控制棒下落時間計算模型

劉 新 陳先龍 張 修 于曉雷

（深圳中廣核工程設(shè)計有限公司 深圳 518172）

反應(yīng)堆驅(qū)動線是反應(yīng)堆控制的“生命線”，執(zhí)行重要的安全功能?？刂瓢粝侣鋾r間不僅是驅(qū)動線設(shè)計考核的重要指標，也是核電站安全分析的關(guān)鍵參數(shù)。基于壓水堆驅(qū)動線的結(jié)構(gòu)特性，對控制棒組件在下落過程進行受力分析，提出落棒時間計算的數(shù)學(xué)模型，采用有限差分法求解運動方程，并且在此基礎(chǔ)上編制計算程序。計算結(jié)果最終表明：程序計算值與嶺澳核電站二期機組熱態(tài)落棒試驗值符合較好，基于此計算方法編制的落棒時間計算程序可以用于壓水堆停堆時的安全評價。

控制棒組件，驅(qū)動線，落棒時間

控制棒組件在反應(yīng)堆中的落棒時間對于反應(yīng)堆的安全運行至關(guān)重要。反應(yīng)堆驅(qū)動線由運動組件和固定組件兩部分組成。運動組件包括控制棒棒束組件(Rod Cluster Control Assembly, RCCA)和驅(qū)動桿，固定組件是運動組件的下降通道，主要包括：控制棒驅(qū)動機構(gòu)(Control Rod Drive Mechanism, CRDM)、熱套管、控制棒導(dǎo)向筒組件和燃料組件導(dǎo)向管。反應(yīng)堆驅(qū)動線各部件結(jié)構(gòu)復(fù)雜，落棒時間計算需要考慮反應(yīng)堆內(nèi)冷卻劑的流速、溫度和壓力的影響，也要考慮外界因素（如地震）的影響。

控制棒在下落過程中，由于通道流體流速的變化，導(dǎo)致受到的水力學(xué)力（水力摩擦力和水力壓力）會不斷的變化，而且變化值與下落位置是非線性關(guān)系，從而增加了計算落棒時間的難度[1]。竇一康等[2]提出了地震工況下用有限差分法對振動方程作離散，對時間歷程則用Wilson-θ法進行逐步分析，對由碰撞引起的非線性采用修正的Newton法作迭代處理的計算方法。但是，常微分方程的 Wilson-θ求解方法比較復(fù)雜，因此，我們采用有限差分方法，將下落時間進行離散，同時利用圓管中流體壓差計算的原理，避開了復(fù)雜的常微分方程計算[3]。此外，目前國內(nèi)外有關(guān)控制棒下落時間計算的軟件大多采用Fortran面向過程的編程語言，不利于代碼的閱讀以及程序的擴展，本文在總結(jié)出更精確、更簡單計算模型的基礎(chǔ)上，采用面向?qū)ο缶幊陶Z言C++編寫出了控制棒落棒時間計算程序[4]。

1 理論模型

控制棒運動組件不僅受到重力還受到阿基米德浮力、水力壓力、水力摩擦力等，其中有些力隨時間進行變化，為了計算落棒時間，我們需要對計算模型進行簡化假設(shè)，并且要把運動過程按照時間步長進行離散，假設(shè)在單個時間步長過程中力是恒定的，然后對時間進行積分計算。

1.1 簡化假設(shè)

(1) 在控制棒導(dǎo)向筒內(nèi)，作用在星型架和控制棒的水力學(xué)力被忽略；

(2) 導(dǎo)向管內(nèi)的水溫為常數(shù)且等于堆芯中冷卻劑的平均溫度；

(3) 在計算的任何時刻，流體應(yīng)為單向的。

1.2 控制棒受力分析

控制棒下落過程受到的力包括恒力和變力（如圖1所示），其中，恒力有重力F1、浮力F2、額定機械摩擦力F3、控制棒導(dǎo)向筒在連續(xù)導(dǎo)向段的水力-機械摩擦力F4，變力包括控制棒在導(dǎo)向管中受到的水力學(xué)力F5、驅(qū)動桿在驅(qū)動機構(gòu)中受到的水力學(xué)力F6、地震工況下的附加力F7。

圖1 控制棒下落過程受力分析Fig.1 Force analysis of control rod in dropping.

根據(jù)牛頓第二定律，有如下動力學(xué)方程：

式中，M為控制棒運動組件的質(zhì)量，kg；a(t)為t時刻控制棒運動組件的加速度，m·s-2。

由式(1)可以看出，需要計算下落時間的關(guān)鍵因素就是變力的計算，因為控制棒運動組件是在流體中運動，我們可以根據(jù)流體在圓管中的流動模型以及圓柱體插入圓管中的流動模型分別進行分析，詳細模型說明見§1.3和§1.4。

1.3 流體在圓管中流動時的模型

假定流體不可壓縮，不考慮重力的影響，流體沿管道流動，如圖2所示。流動截面上的靜壓和壁面摩擦力均勻分布，則動量方程為：

式中， Cf為壁面摩擦系數(shù)；U為流體的流速，m·s-1；Pe為圓管的濕周，m；S為管道的截面積，m2。

圖2 流體在圓管中流動模型Fig.2 Fluid model of flows in a pipe.

則動量方程[5]可寫成：

進行空間積分：

1.4 圓柱體在圓管中運動時的流動模型

1.4.1 壓降計算

圓柱體在管道中運動，如圖3所示。假設(shè)圓柱體在管道中的下落速度為V，根據(jù)式(2)得出流體此時的動量方程：

式中，Cf′為圓柱體的壁面摩擦系數(shù)；Pe'為圓柱體的濕周，m；Sa為圓柱體與壁面環(huán)腔之間的截面積，m2。

通過對動量方程積分得：

圖3 圓柱體在圓管中流動模型Fig.3 Fluid model of flow in a pipe with a mobile cylinder.

1.4.2 圓柱體受力分析

流體和圓柱體間的水力學(xué)摩擦力通過壁面摩擦系數(shù)Cf′來確定：

由于流體流動產(chǎn)生的壓降而作用在圓柱體頂部和底部的壓差阻力根據(jù)式(8)進行計算：式中，SCR為圓柱體的截面積，m2；zΔ為控制體的步長，m。

1.5 控制棒組件下落過程中變力分析

控制棒組件在下落通道中受到的流體阻力計算主要的難點為驅(qū)動機構(gòu)中驅(qū)動桿受到的流體阻力和控制棒在燃料組件導(dǎo)向管中受到的流體阻力的計算，根據(jù)§1.3和§1.4的物理模型，我們采取下面的思路進行計算。

1.5.1 驅(qū)動機構(gòu)區(qū)域

重點在于計算驅(qū)動桿在驅(qū)動機構(gòu)區(qū)域流體流速的計算，把驅(qū)動結(jié)構(gòu)分為三段求解，分別為行程套管區(qū)、鉤爪區(qū)、熱套管區(qū)。在三個作用區(qū)域分別利用流體的質(zhì)量守恒方程和§1.4給出的壓差和受力計算方程聯(lián)立得出非線性方程組，然后采用牛頓迭代法進行求解，最后計算出驅(qū)動桿在行程套管區(qū)、鉤爪區(qū)、熱套管區(qū)三個區(qū)域受到的流體水力作用力。

1.5.2 燃料組件導(dǎo)向管區(qū)域

重點在于計算燃料組件導(dǎo)向管中控制棒下方流體的流速以及控制棒與導(dǎo)向管環(huán)形通道中流體流速。根據(jù)§1.3和§1.4壓差的計算方程，得出導(dǎo)向管內(nèi)流體的壓差的方程，然后根據(jù)燃料組件外部的壓差與內(nèi)部壓差的關(guān)系，計算出導(dǎo)向管內(nèi)各處流體的流速，最后計算出控制棒在導(dǎo)向管內(nèi)受到的流體水力作用力。

2 算法分析

從初始時刻t(0)開始計算，按時間步長迭代計算。初始狀態(tài)工況的描述如下：

(1) 移動組件的初速度為零；

(2) 導(dǎo)向管內(nèi)的流體流量是恒定的；

(3) 在管座/熱套管/驅(qū)動機構(gòu)部位的流體流量為零；

(4) 控制棒被提出或被部分提出。

在 t(i) 時刻，我們計算每個產(chǎn)生作用的力的值。我們從中得出加速度，并且假設(shè)在單個時間步長中力是恒定的，得出 t(i+1) 時刻的速度和行程?；谶@些值，重新計算變力。

由上可以看出，在每個時間步長之后，需要重新計算的力有導(dǎo)向管內(nèi)控制棒的水力學(xué)力F5、驅(qū)動機構(gòu)部位驅(qū)動桿的水力學(xué)力F6、附加力F7（地震工況下才需要計算）。我們計算控制棒從任意給定位置依靠重力下落到導(dǎo)向管緩沖段入口的時間以及星型架底端彈簧剛好接觸上管座上表面的時間。

在地震工況下，需借助試驗數(shù)據(jù)通過插值的方法得出t時刻控制棒組件受到的地震附加力[6]，但尚未進行相關(guān)的模擬振動試驗，所以暫時預(yù)留接口。

利用面向?qū)ο蟮乃枷?，對控制棒下落過程進行模型抽象化，創(chuàng)建阻力模型管理器，把下落過程分解出的阻力計算模型逐個添加到阻力模型管理器中。計算流程如圖4所示。

3 程序架構(gòu)

根據(jù)理論模型的建立和算法的分析，設(shè)計控制棒下落時間計算軟件的架構(gòu)如圖5所示。

核心計算程序采用C++編寫，運行環(huán)境為Linux系統(tǒng)。

4 算例驗證

為了驗證理論算法的正確性，本文把編制程序的計算結(jié)果與嶺澳二期核電站4號機組熱態(tài)落棒試驗實測數(shù)據(jù)進行對比，如表1所示。程序計算得出的控制棒下落過程加速度、速度、位移隨時間的變化曲線如圖6所示，該程序能夠較準確地計算壓水反應(yīng)堆控制棒的下落時間，而且計算值較試驗值偏大，說明程序計算值偏保守。

圖4 程序流程圖Fig.4 Program flow chart.

圖5 程序架構(gòu)圖Fig.5 Architecture diagram.

表1 落棒時間計算結(jié)果和試驗結(jié)果對比Table 1 Drop time comparison between calculation results and test results.

圖6 控制棒下落過程加速度(a)、速度(b)和位移(c)的時程曲線Fig.6 Acceleration (a), velocity (b) and displacement (c) time history curve of control rod falling.

5 結(jié)語

根據(jù)對落棒時間下落過程的機理分析，利用流體力學(xué)的動量守恒方程，并且采用有限差分法，將復(fù)雜的動力學(xué)方程轉(zhuǎn)化成容易求解的偏微分方程，利用時間的一階向后差分方式進行求解計算。此外，采用面向?qū)ο蟮母呒壵Z言，對下落過程的各個受力模型進行抽象化處理，使得編制的程序更具有閱讀性和擴展性。

通過程序計算結(jié)果與嶺澳核電站二期機組熱態(tài)落棒試驗值進行比較，表明本文總結(jié)的落棒時間計算模型和計算方法具有可靠性，該程序可用于壓水堆停堆時的安全評價。

1 孫磊, 于建華, 魏永濤，等. 控制棒組件落棒時間與歷程計算[J]. 核動力工程, 2003, 24(1): 59-62 SUN Lei, YU Jianhua, WEI Yongtao, et al. Analysis of drop-time and course of control rod assembly[J]. Nuclear Power Engineering, 2003, 24(1): 59-62

2 竇一康, 姚偉達, 楊仁安, 等. 反應(yīng)堆控制棒在事故工況下的落棒時間分析[F]. 中國核科技報告, 1998 DOU Yikang, YAO Weida, YANG Ren’an, et al. Computation of reactor control rod drop time under accident conditions[F]. China Nuclear Science and Technology Report, 1998

3 李慶揚, 王能超, 易大義. 數(shù)值分析[M]. 北京: 清華大學(xué)出版社, 2008 LI Qingyang, WANG Nengchao, YI Dayi. Numerical analysis[M]. Beijing: Tsinghua University Press, 2008

4 Lippman S B, Lajoie J, Moo B M, et al. C++ primer[M]. Beijing: Electronic Industry Press, 2013

5 陳卓如. 工程流體力學(xué)[M]. 北京: 高等教育出版社, 2004 CHEN Zhuoru. Engineering fluid mechanics[M]. Beijing: Higher Education Press, 2004

6 屈鐵軍, 陳厚群. 強地震作用下核電站控制棒下落時間分析[J]. 計算力學(xué)學(xué)報, 2001, 18(3): 364-370 QU Tiejun, CHEN Houqun. Control rod drop time analysis for the nuclear power plant under action of strong earthquake[J]. Chinese Journal of Computational Mechanics, 2001, 18(3): 364-370

CLCTL364

Calculation model of controlling rod drop time

LIU Xin CHEN Xianlong ZHANG Xiu YU Xiaolei
(China Nuclear Power Design Company Ltd., Shenzhen 518172, China)

Background:Reactor drive line is the "lifeline" of reactor control, which performs important safety functions.Purpose:To provide the basis for the drive line design assessment and nuclear power plant safety analysis.Methods:Through analyzing the stress of pressurized water reactor control rod assembly and control rod drive mechanism, established the mathematical model and the motion equation, and then solved by the finite difference method.Results:The program calculated value accords well with the test value of Ling’ao nuclear power plant phaseⅡ unit hot rod drop.Conclusion:The rod drop time of calculation program based on this method can be used to the safety evaluation when pressurized water reactor is shut down.

Controlling rod assembly, Drive line, Rod drop time

TL364

10.11889/j.0253-3219.2014.hjs.37.110604

劉新，女，1986年出生，2009年畢業(yè)于南華大學(xué)，現(xiàn)主要從事反應(yīng)堆非堆芯熱工水力計算軟件研發(fā)，E-mail: liuxin_067@163.com

2014-04-04，

2014-05-08