Analysis of Influence of Expansion on Control Rod Drop in Nuclear Reactor at High Temperature
-
摘要: 核反应堆内的高温环境会导致控制棒组件产生热膨胀现象，进而影响控制棒组件的下落时间。为研究热膨胀对控制棒组件下落的影响，建立考虑膨胀变形的控制棒下落流固耦合模型，求解控制棒组件下落过程的动力学方程。采用有限元对导向管和控制棒进行稳态热力耦合分析，得到高温条件下控制棒和导向管的膨胀变形。基于所建立模型比较考虑热膨胀和不考虑热膨胀条件下控制棒的下落动力学过程。分析结果表明，热膨胀现象延迟了控制棒整体下落时间，但对控制棒进入缓冲段入口的时间影响较小，因此在一般情况下热膨胀对落棒时间的影响可以忽略，但在工程中仍需要重视控制棒膨胀引起的后果。本文研究结果对控制棒结构设计以及控制棒堆内落棒时间分析有重要意义。Abstract: The high temperature environment in the nuclear reactor will cause the thermal expansion of the control rod assembly, which will affect the falling time of the control rod assembly. In order to study the effect of thermal expansion on the drop of control rod assemblies, a fluid-structure coupling model for the drop of control rod assemblies was established considering expansion deformation, and the dynamic equations for the drop process of control rod assemblies were solved. The steady state thermodynamic coupling analysis of the guide tube and control rod is performed using finite element method to obtain the expansion deformation of the control rod and guide tube under high temperature conditions. Based on the established analytical model, the dynamic processes of control rod drop with and without expansion is compared. The analysis results show that the thermal expansion phenomenon delays the total drop time of the control rod, but has little impact on the drop time before the buffer section. Therefore, in general, the impact of thermal expansion on the drop time can be ignored, but the consequences caused by the control rod expansion still need to be paid attention to in the engineering design. The conclusion of this paper is important for the design of control rod and the analysis of in-pile control rod drop time.
-
Key words:
- Control rod drop /
- Thermal expansion /
- Finite element /
- Fluid-structure coupling
-
0. 引 言
核反应堆安全性是核能发展的重要前提，而控制棒的快速下落是保证核反应堆安全停堆的关键。如果控制棒在下落过程中出现卡棒行为或者下落时间过长的情况，会导致非常严重的核反应堆事故，因此研究控制棒下落非常重要。关于控制棒下落这一问题的研究，国内外学者做了许多工作。文献[1]采用有限元方法对压水堆中控制棒下落过程进行了分析，解决了下落过程中控制棒的流固耦合问题。文献[2]采用多点动力学方法研究了单根控制棒跌落事故，建立了瞬态分析耦合动力学方程。文献[3-4]研究了铅铋快堆以及超临界水冷堆控制棒的下落过程，利用计算流体动力学（CFD）瞬态动网格数值分析方法得到了控制棒下落位移、速度以及加速度与时间的关系。
然而包括上述在内的大多数研究只是分析了控制棒原始结构的下落，没有考虑控制棒结构在反应堆环境下发生变化（如热膨胀）的情况。在实际工程中，反应堆内的高温高压、辐射环境可能会使控制棒出现热膨胀现象。若膨胀显著影响了控制棒的下落时间，则核反应堆安全性难以保证，因此研究膨胀对控制棒下落的影响十分重要。文献[5]统计了控制棒和驱动机构组件的异常事件，其中由于控制棒膨胀或变形造成的事故占比高达17%。文献[6-7]研究表明控制棒在反应堆内部的高温高压、辐射环境下可能会出现组织结构和热稳定性改变以及高温蠕变的现象。文献[8]利用ANSYS软件对控制棒驱动机构进行了稳态热力学分析，计算了在交变温度条件下产生蠕变疲劳损伤的热应力。文献[9]使用增量连续塑性原理开发本构方程对反应堆堆芯材料的辐射膨胀进行研究，并对堆芯的膨胀变形和干扰问题进行了工程评估。文献[10]在钠冷快堆系统分析中开发了能够处理堆芯内热膨胀效应的模块，发现热膨胀效应在瞬态演变中起主要作用。文献[11]建立了金属燃料钠冷快堆堆芯径向膨胀的实验模型。
上述文献从不同角度研究了控制棒热膨胀现象，然而针对热膨胀现象对控制棒下落时间的具体影响却鲜有研究。因此，本文主要研究热膨胀现象对控制棒下落过程的影响程度。首先通过ANSYS软件分析高温下控制棒组件的热膨胀应变，然后将应变结果代入到本文建立的控制棒下落动力学方程中，得到热膨胀后控制棒的下落行为。
1. 理论模型
1.1 热膨胀模型
反应堆内，控制棒与导向管在高温条件下长期运行会因为高温作用引起热膨胀现象。本文只考虑控制棒与导向管的径向膨胀。由于导向管的外部膨胀对落棒时间影响较小，本文只考虑高温条件下导向管沿着内径向内膨胀，以及控制棒沿径向向外膨胀的影响。不同高度的堆芯温度不同，控制棒与导向管的热膨胀大小也会与高度相关，从而可以得到变直径的控制棒和导向管。由此，建立如图1所示的控制棒下落热膨胀模型，导向管向内膨胀未画出。以图1导向管底部中心为原点o，向上为正建立oz坐标系。图1a中，控制棒底部在流水孔上方，导向管分为导向管底部位置Z0至流水孔位置Z1、流水孔位置至控制棒底部位置Z、控制棒底部位置至导向管顶部位置Z2三部分，区域分别为$ [{Z_0},{Z_1}] $、$[{Z_1},Z]$和$ [Z,{Z_2}] $。在该位置下，流体${Q_0}$从底部孔流入，经流水孔后分为两部分。一部分从流水孔中流出，记为${Q_1}$；另一部分继续向上流动，记为${Q_{{\text{1\_Z}}}}$；当流动到控制棒底部后，流体则在控制棒与导向管组成的环形管道中流动，记为${Q_2}$(z)。图1b中，控制棒底部在流水孔下方，这时导向管分为导向管底部至控制棒底部、控制棒底部至流水孔、流水孔至导向管顶部三部分，区域分别为$[{Z_0},Z]$、$[Z,{Z_1}]$和$[{Z_1},{Z_2}]$。该图中，流体${Q_0}$从底部流入之后，会先进入控制棒与导向管组成的环形管道，记为${Q'_{{\text{1\_Z}}}}$(z)，再经过流水孔后分为两部分，一部分${Q_1}$从流水孔流出，另一部分环形管道中的流体${Q_2}$(z)继续向上流动。导向管内的流体满足不可压缩条件。
图 1 控制棒下落热膨胀模型黄色部分—控制棒向外膨胀部分；P—压强；v—控制棒下落速度Figure 1. Thermal Expansion Model of Control Rod Drop为方便计算分析，将上述模型简化为圆管和圆环模型。当流体在导向管内流动时，可以看作圆管模型；当流体在导向管与控制棒环形通道中流动时，可以看作圆环模型。
1.2 模型简化
1.2.1 圆管模型
导向管内流体流动如图2所示，流体从导向管左侧流向右侧，流动满足N-S方程[12]，见式（1）。取导向管内的两个点分别为$ {Z_{\text{a}}} $和$ {Z_{\text{b}}} $，将式（1）在$ [{Z_{\text{a}}},{Z_{\text{b}}}] $上积分，可得压降方程式（2）。
$$ \begin{aligned}[b]\frac{1}{S_{\text{1}}}\frac{\partial Q_{\text{t}}}{\partial t}-\frac{Q_{\text{t}}^{\text{2}}}{S_{\text{1}}^3}\frac{\partial S_{\text{1}}}{\partial\text{z}}= & -\frac{1}{\rho}\frac{\partial P}{\partial\text{z}}-\mathrm{sgn}\left(Q_{\text{t}}\right)\frac{C_{\text{GT}}l_{\text{GT}}}{2}\frac{Q_{\text{t}}^{\text{2}}}{S_{\text{1}}^3}- \\ & \mathrm{sgn}\left(v_{\text{f}}\right)\times\frac{1}{2}K_{\text{c}}\left(v_{\text{f}}S_1\right)^2\end{aligned} $$ (1) 式中，$ {S_{\text{1}}} $为导向管横截面面积；$\rho $为流体密度；$ {l_{{\text{GT}}}} $为管壁周长；${{\mathrm{sgn}}} $为符号判断函数；$ {C_{{\text{GT}}}} $为管壁摩擦系数；$ {v_{\text{f}}} $为液体流速；$ {K_{\text{c}}} $为形阻系数；$ Q_{\text{t}} $为流量；t为时间。
$$ \begin{aligned}[b] {P_{{{\text{Z}}_{\text{b}}}}} - {P_{{{\text{Z}}_{\text{a}}}}} = & - \rho \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {\frac{{{{\dot Q}_{\text{t}}}}}{{{S_{\text{1}}}}}} {\text{d}}{\textit{z}} - Q_{\text{t}}^{\text{2}}\frac{\rho }{2}\left[ {\frac{1}{{S_{\text{1}}^2\left( {{Z_{\text{b}}}} \right)}} - \frac{1}{{S_{\text{1}}^2\left( {{Z_{\text{a}}}} \right)}}} \right] -\\ & \frac{\rho }{2}{{\mathrm{sgn}}} \left( {{Q_{\text{t}}}} \right)\int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{C_{{\text{GT}}}}{{\left( {\frac{{{Q_{\text{t}}}}}{{{S_{\text{1}}}}}} \right)}^2}\frac{{{l_{{\text{GT}}}}}}{{{S_{\text{1}}}}}{\text{d}}{\textit{z}}} -\\ & Q_{\mathrm{t}}^2\frac{\rho }{2}\int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{{\mathrm{sgn}}} \left( {\frac{{{Q_{\mathrm{t}}}}}{{{S_1}}}} \right){K_{\text{c}}}{\text{dz}}} \end{aligned} $$ (2) 令：
$$ \left\{ \begin{gathered} {I_1}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \rho \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {\frac{{{\text{dz}}}}{{{S_{\text{1}}}}}} \\ {I_2}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \frac{\rho }{2}\left[ {\frac{{Q_{\text{t}}^{\text{2}}\left( {{Z_{\text{b}}}} \right)}}{{S_{\text{1}}^2\left( {{Z_{\text{b}}}} \right)}} - \frac{{Q_{\text{t}}^{\text{2}}\left( {{Z_{\text{a}}}} \right)}}{{S_{\text{1}}^2\left( {{Z_{\text{a}}}} \right)}}} \right] \\ {I_3}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \frac{\rho }{2}{{\mathrm{sgn}}} \left( {{Q_{\text{t}}}} \right)\int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{C_{{\text{GT}}}}{{\left( {\frac{{{Q_{\text{t}}}}}{{{S_{\text{1}}}}}} \right)}^2}\frac{{{l_{{\text{GT}}}}}}{{{S_{\text{1}}}}}{\text{dz}}} \\ {I_5}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = {Q_{\mathrm{t}}^2}\frac{\rho }{2}\int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{{\mathrm{sgn}}} \left( {{v_{\text{f}}}} \right){K_{\text{c}}}{\text{dz}}} \\ \end{gathered} \right. $$ (3) 式中，$ {I_1} $(Za, Zb)为由于流体加速度导致的压降系数（惯性压降系数）；$ {I_2} $(Za, Zb)为流体对流产生的压降（对流压降）；$ {I_3} $(Za, Zb)为导向管壁摩擦压降；$ {I_5} $(Za, Zb)为导向管直径改变导致的压降。
则式（2）简化为：
$$ {P_{{{\text{Z}}_{\text{b}}}}} - {P_{{{\text{Z}}_{\text{a}}}}} = - {\dot Q_{\text{l}}}{I_1}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) - {I_2}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) - {I_3}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) - {I_5}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) $$ (4) 1.2.2 圆环模型
当导向管内存在控制棒时，流体在导向管与控制棒之间的环形截面内流动，如图3所示。同样地，流体流动满足N-S方程[12]，即：
$$ \begin{aligned}[b] \frac{{\partial {v_{\text{f}}}\left( {\textit{z}} \right)}}{{\partial t}} +& {v_{\text{f}}}\left( {\textit{z}} \right)\frac{{\partial {v_{\text{f}}}\left( {\textit{z}} \right)}}{{\partial {\textit{z}}}} =- \frac{1}{\rho }\frac{{\partial P}}{{\partial {\textit{z}}}} -\\ & {{\mathrm{sgn}}} \left[ {{v_{\text{f}}}\left( {\textit{z}} \right)} \right]{C_{{\text{CR}}}}\frac{{v_{_{\text{f}}}^{\text{2}}\left( {\textit{z}} \right)}}{2}\frac{{{l_{{\text{CR}}}}}}{{{S_2}}} -\\ & {{\mathrm{sgn}}} \left[ {{v_{\text{f}}}\left( {\textit{z}} \right) - v} \right]{C_2}\frac{{{{\left[ {{v_{\text{f}}}\left( {\textit{z}} \right) - v} \right]}^2}}}{2}\frac{{{l_{{\text{CR}}}}}}{{{S_2}}} -\\ & {{\mathrm{sgn}}} \left[ {{v_{\text{f}}}\left( {\textit{z}} \right)} \right] \times \frac{1}{2}{K_{\text{c}}}{\left[ {{v_{\text{f}}}\left( {\textit{z}} \right){S_2}} \right]^2} \end{aligned} $$ (5) 式中，$ {C_{{\text{CR}}}} $为控制棒摩擦系数；$ {l_{{\text{CR}}}} $为控制棒周长；$ {S_2} $为导向管和控制棒组成的环形面积；vf(z)为$ {\textit{z}} $位置处流体流速；C2为外圆管壁的摩擦系数。
根据不可压缩条件，可知导向管和控制棒之间的流量$ {Q_{\text{n}}}\left( {\textit{z}} \right) $满足如下关系：
$$ {Q_{\text{n}}}\left( {\textit{z}} \right) = {Q_{\text{m}}} - v{S_{{\text{CR}}}}\left( {\textit{z}} \right) $$ (6) 式中，$ {S_{{\text{CR}}}}\left({\textit{z}} \right) $为$ {\textit{z}} $位置处控制棒横截面积。
同时可知流量和流速满足关系式$ {v_{\text{f}}}\left( {\textit{z}} \right) = \dfrac{{{Q_{\text{n}}}\left( {\textit{z}} \right)}}{{{S_2}\left( {\textit{z}} \right)}} $，式（5）可以变形为：
$$ \begin{aligned}[b] & \frac{1}{S_2}\left[\dot{Q}_{\text{m}}-S_{\text{CR}}\dot{v}\right]+\frac{Q_{\text{n}}\left( {\textit{z}} \right)}{S_2}\frac{\partial}{\partial {\textit{z}} }\left[\frac{Q_{\text{n}}\left( {\textit{z}} \right)}{S_2}\right] \\ & =-\frac{1}{\rho}\frac{\partial P}{\partial {\textit{z}}}-\mathrm{sgn}\left[Q_{\text{n}}\left( {\textit{z}} \right)\right]\frac{C_{\text{GT}}}{2}\frac{l_{\text{GT}}}{S_2}\left[\frac{Q_{\text{n}}\left( {\textit{z}} \right)}{S_2}\right]^2- \\ & \mathrm{sgn}\left[\frac{Q_{\text{n}}\left( {\textit{z}} \right)}{S_2}-v\right]\frac{C_{\text{CR}}}{2}\frac{l_{\text{CR}}}{S_2}\left[\frac{Q_{\text{n}}\left( {\textit{z}} \right)}{S_2}-v\right]^2\end{aligned} $$ (7) 取圆环中的两点${Z_{\text{a}}}$和${Z_{\text{b}}}$，将式（7）在$\left[ {{Z_{\text{a}}},{Z_{\text{b}}}} \right]$上积分可得压降方程：
$$ \begin{aligned}[b] & P\left( {{Z_{\text{a}}}} \right) - P\left( {{Z_{\text{b}}}} \right) = - \ddot {\textit{z}} {I_0}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) +{\dot Q_{\text{m}}}{I_1}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + \\[-2pt] & Q_{\text{m}}^2{I_2}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) +{I_3}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {I_4}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {I_5}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) \end{aligned} $$ (8) $$ \left\{ \begin{gathered} {I_0}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \rho \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {\frac{{{S_{{\text{CR}}}}}}{{{S_2}}}{\text{d}}{\textit{z}}} \; \\[-2pt] {I_4}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \frac{\rho }{2}\int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{{\mathrm{sgn}}} \left[ {\frac{{{Q_{\text{n}}\left( {\textit{z}} \right)}}}{{{S_2}}} - v} \right]{C_{{\text{CR}}}}{{\left[ {\frac{{{Q_{\text{n}}\left( {\textit{z}} \right)}}}{{{S_2}}} - v} \right]}^2}\frac{{{l_{{\text{CR}}}}}}{{{S_2}}}{\text{d}}{\textit{z}}} \\ \end{gathered} \right. $$ (9) 式中，${I_0}$(Za, Zb)为控制棒直径改变导致的压降系数；${I_4}$(Za, Zb)为控制棒壁面摩擦导致的压降。${I_1}$(Za, Zb)、${I_2}$(Za, Zb)、${I_3}$(Za, Zb)和${I_5}$(Za, Zb)表达式见式（3），只需将式中Ql换成Qn(z)，${S_1}$换成${S_2}$。
1.3 控制棒下落动力学方程
控制棒运动受牛顿第二定律控制，其控制方程可以写为：
$$ {F_{{\text{CR}}}} = {M_{{\text{CR}}}} \times \ddot {\textit{z}} $$ (10) 式中，MCR为控制棒的质量；${F_{{\text{CR}}}}$为控制棒受到的合力。
$$ {F_{{\text{CR}}}} = {F_{\text{G}}} + {F_{\text{b}}} + {F_{{\text{f\_CR}}}} + {F_{{\text{r\_CR}}}} $$ (11) 式中，FG、Fb分别为重力、浮力（已知条件）；${F_{{\text{f\_CR}}}}$、${F_{{\text{r\_CR}}}}$分别为控制棒受到的摩擦阻力和压差阻力（需要对其求解）。
根据式（9）中的${I_4}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right)$，可以得到在区间$\left[ {{Z_{\text{a}}},{Z_{\text{b}}}} \right]$内控制棒受到的摩擦阻力为：
$$ {F_{{\text{f\_CR}}}} = \frac{\rho }{2}\int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{{\mathrm{sgn}}} \left[ {\frac{{{Q_{\text{n}}\left( {\textit{z}} \right)}}}{{{S_2}}} - v} \right]{C_{{\text{CR}}}}{{\left[ {\frac{{{Q_{\text{n}}\left( {\textit{z}} \right)}}}{{{S_2}}} - v} \right]}^2}{l_{{\text{CR}}}}} {\text{d}}{\textit{z}} $$ (12) 区间$\left[ {{Z_{\text{a}}},{Z_{\text{b}}}} \right]$内控制棒受到的压差阻力为：
$$ {F_{{\text{r\_CR}}}} = \left( {{P_{{{\text{Z}}_{\text{a}}}}} - {P_{{{\text{Z}}_{\text{b}}}}}} \right){S_{{\text{CR}}}}\left( {{Z_{\text{a}}}} \right) + \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {\left[ {P\left( {\textit{z}} \right) - {P_{{{\text{Z}}_{\text{b}}}}}} \right]} {\text{d}}{S_{{\text{CR}}}}\left( {\textit{z}} \right) $$ (13) 式中，${S_{{\text{CR}}}}\left( {{Z_{\text{a}}}} \right)$为控制棒在${Z_{\text{a}}}$处的横截面积。
结合式（8）与式（13）可得：
$$ \begin{aligned}[b] {F_{{\text{r\_CR}}}} = & \big[ - \dot v{I_0}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {{\dot Q}_{\text{m}}}{I_1}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {I_2}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + \\ & {I_3}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {I_4}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {I_5}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) \big]{S_{{\text{CR}}}}\left( {{Z_{\text{a}}}} \right) - \\ &\dot vI{I_0}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + {\dot Q_{\text{m}}}I{I_1}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + I{I_2}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) +\\ &I{I_3}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) + I{I_4}\left( {{Z_{\text{a}}},{Z_b}} \right) + I{I_5}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) \end{aligned} $$ (14) $$ \left\{ \begin{gathered} I{I_0}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{I_0}\left( {{\textit{z}},{Z_{\text{b}}}} \right)} {\text{d}}{{{S}}_{{\text{CR}}}}\left( {\textit{z}} \right) \\ I{I_1}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{I_1}\left( {{\textit{z}},{Z_{\text{b}}}} \right)} {\text{d}}{{{S}}_{{\text{CR}}}}\left( {\textit{z}} \right) \\ I{I_2}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{I_2}\left( {{\textit{z}},{Z_{\text{b}}}} \right)} {\text{d}}{{{S}}_{{\text{CR}}}}\left( {\textit{z}} \right) \\ I{I_3}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{I_3}\left( {{\textit{z}},{Z_{\text{b}}}} \right)} {\text{d}}{{{S}}_{{\text{CR}}}}\left( {\textit{z}} \right) \\ I{I_4}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{I_4}\left( {{\textit{z}},{Z_{\text{b}}}} \right)} {\text{d}}{{{S}}_{{\text{CR}}}}\left( {\textit{z}} \right) \\ I{I_5}\left( {{Z_{\text{a}}},{Z_{\text{b}}}} \right) = \int_{{Z_{\text{a}}}}^{{Z_{\text{b}}}} {{I_5}\left( {{\textit{z}},{Z_{\text{b}}}} \right)} {\text{d}}{{{S}}_{{\text{CR}}}}\left( {\textit{z}} \right) \\ \end{gathered} \right. $$ (15) 由此，根据简化模型得到的受力公式进一步得到图1下落过程中控制棒与导向管的受力。
压降平衡方程（4）、（8）和受力方程（10）共同组成了控制棒下落和流体流动的流固耦合模型。此外，为使方程组封闭，需要添加流量代数方程组。当控制棒底部位于流水孔上方时，流量方程为：
$$ \left\{ \begin{gathered} {Q_0} = {Q_1} + {Q_{{\text{1\_Z}}}} \\ {Q_2}\left( {\textit{z}} \right) = {Q_{{\text{1\_Z}}}} - v{S_{{\text{CR}}}}\left( {\textit{z}} \right),\;\;Z \geqslant {Z_1} \\ \end{gathered} \right. $$ (16) 当控制棒底部位于流水孔下方时，流量方程为：
$$ \left\{ \begin{gathered} {{Q'}_{{\text{1\_Z}}}}\left( {\textit{z}} \right) = {Q_0} - v{S_{{\text{CR}}}}\left( {\textit{z}} \right) \\ {Q_2}\left( {\textit{z}} \right) = {Q_1} + {{Q'}_{{\text{1\_Z}}}}\left( {\textit{z}} \right),\;\;Z < {Z_1} \\ \end{gathered} \right. $$ (17) 将流量代数方程组和微分方程组合并，得到控制棒下落控制方程，并采用辛Euler中点格式[13]对控制棒下落进行求解，可以得到控制棒下落过程中每一个时间步的位移、速度和加速度。
2. 稳态热力耦合分析
根据压水堆内结构相对位置可以得到导向管与控制棒的温度区间，如图4所示。此外，在使用ANSYS软件进行热膨胀分析时，还需要知道导向管与控制棒的材料及材料参数，如密度、弹性模量、线膨胀系数和热导率等。
图 4 控制棒和导向管的轴向温度分布图Figure 4. Axial Temperature Distribution of Control Rod and Guide Tube导向管材料选用压水堆常用的Zr-4合金，密度为6560 kg/m3，弹性模量为95.2 GPa，泊松比为0.33[14]。根据文献[15-16]可得到20℃~800℃范围内Zr-4合金相对热膨胀的经验公式和298~1073 K范围内Zr-4合金热导率的计算公式。
控制棒包壳材料选用常用的316钢，密度为7980 kg/m3，弹性模量为210.3 GPa，泊松比为0.3。在20℃~400℃的平均热膨胀系数为1.75×10−5℃−1，平均热导率为1.84 W/(m·K)。
堆芯结构和流量相关参数如表1所示。本文将上述物理参数输入到ANSYS软件中，根据控制棒和导向管真实参数建立控制棒和导向管仿真模型，并划分网格，如图5和图6所示（原模型结构细长，图中只展示部分长度）。控制棒由于是规整的圆柱体结构，采用六面体单元，网格大小设置为$1.6 \times {10^{ - 3}}{\text{ m}}$。同理，导向管也采用六面体单元，网格大小设置为$7 \times {10^{ - 4}}{\text{ m}}$。
表 1 堆芯结构和流量相关参数Table 1. Parameters of Core Structure and Flow Rate参数名 数值 控制棒长度/m 4 控制棒直径/mm 9.7 导向管底部长度/m 0.5 导向管底部内径/mm 10 导向管顶部长度/m 3.5 导向管顶部内径/mm 11.4 侧孔距底部距离/m 0.52 反应堆内每个回路内的流量/(m3·s−1) 1.8×10−3 按照上述假定的温度区间，在ANSYS软件稳态热力学模块中给控制棒和导向管施加温度。为了能够求解该温度区间内的热膨胀位移，出于保守分析考虑，对导向管外侧施加固定约束，使热膨胀后管道尽可能缩小。通过计算得到了控制棒和导向管的径向变形云图，如图7所示。图7a从整体来看控制棒径向变形是轴对称的，且沿径向向外膨胀，当离圆心越远时，径向膨胀变形越大。图7b中，由于温度沿长度方向变化，因此轴截面的径向变形也沿长度方向变化。图7c中显示，变形满足轴对称，但由于外壁固定，整体沿径向向内膨胀，且越靠近圆心膨胀位移越大。
由图7可以得到控制棒和导向管沿长度方向变化的热膨胀位移，进而描述出控制棒和导向管热膨胀后的直径变化与轴向坐标的关系，如图8所示。其中，导向管有缓冲段以及非缓冲段部分，由于这两部分的直径与厚度均不相同，在研究热膨胀后的直径变化与轴向坐标关系时需要分别单独分析。由图8a可知，由于控制棒向外膨胀，所以直径变化与轴向坐标呈正相关。由图8b和图8c可知，导向管直径变化与轴向坐标呈负相关，是因为沿径向是向内膨胀的。
图 8 不同位置直径变化与轴向坐标关系图Figure 8. Relationship between Diameter Change at Different Positions and Axial Coordinate3. 热膨胀对控制棒下落的影响
根据上述计算得到控制棒和导向管热膨胀后尺寸，将其代入到控制棒下落动力学方程中，求解得到控制棒下落时间以及每个时间步下控制棒的位移、速度，如图9所示。由图9可知，热膨胀后的下落时间大于热膨胀前的，说明热膨胀会延长控制棒下落时间。同时热膨胀后的速度曲线在前半段略低于热膨胀前，说明热膨胀增加了控制棒下落时的阻力，从而使控制棒花费更长的时间下落到底部。
图 9 热膨胀前后控制棒下落过程对比图Figure 9. Comparison of Rod Drop Processes before and after Expansion将膨胀前后计算得到的2个关键参数：控制棒进入缓冲段的时间和控制棒整体落棒时间进行对比，结果如表2所示。由表2可知，热膨胀后控制棒进入缓冲段的时间增加了0.74%，整体落棒时间增加了3.28%。因此，热膨胀对控制棒进入缓冲段时间的影响可以忽略，对整体落棒时间有一定影响。
表 2 热膨胀前后落棒时间相对偏差Table 2. Relative Deviation of Control Rod Time before and after Thermal Expansion对比项 热膨胀前 热膨胀后 相对偏差 控制棒进入缓冲段时间 1.4270 s 1.4375 s 0.74% 控制棒整体落棒时间 1.8587 s 1.9197 s 3.28% 4. 结 论
针对热膨胀对落棒过程的影响，本文建立了热膨胀条件下控制棒下落流固耦合模型，重点分析了控制棒受到的摩擦阻力、压降阻力，采用有限元方法对导向管和控制棒进行稳态热力耦合分析。得到结论如下：
（1）热膨胀会增加控制棒下落时的阻力，延长下落时间。最终热膨胀导致控制棒进入缓冲段的时间增加了0.74%，整体落棒时间增加了3.28%。
（2）落棒时间是核反应堆事故安全分析的重要输入之一，控制棒进入缓冲段时间的增加与堆芯负反应性引入密切相关。由于热膨胀对控制棒进入缓冲段时间的影响较小，因此在一般情况下热膨胀对落棒时间的影响可以忽略。但是在实际工程中，控制棒组件会受到高压、辐射等因素影响，其膨胀比例可能会更大甚至存在异常肿胀，对下落时间造成更为严重的影响，因此在工程中仍需要重视控制棒膨胀引起的后果。
-
图 1 控制棒下落热膨胀模型
黄色部分—控制棒向外膨胀部分；P—压强；v—控制棒下落速度
Figure 1. Thermal Expansion Model of Control Rod Drop
图 4 控制棒和导向管的轴向温度分布图
Figure 4. Axial Temperature Distribution of Control Rod and Guide Tube
图 8 不同位置直径变化与轴向坐标关系图
Figure 8. Relationship between Diameter Change at Different Positions and Axial Coordinate
图 9 热膨胀前后控制棒下落过程对比图
Figure 9. Comparison of Rod Drop Processes before and after Expansion
表 1 堆芯结构和流量相关参数
Table 1. Parameters of Core Structure and Flow Rate
参数名 数值 控制棒长度/m 4 控制棒直径/mm 9.7 导向管底部长度/m 0.5 导向管底部内径/mm 10 导向管顶部长度/m 3.5 导向管顶部内径/mm 11.4 侧孔距底部距离/m 0.52 反应堆内每个回路内的流量/(m3·s−1) 1.8×10−3 
下载: 导出CSV
表 2 热膨胀前后落棒时间相对偏差
Table 2. Relative Deviation of Control Rod Time before and after Thermal Expansion
对比项 热膨胀前 热膨胀后 相对偏差 控制棒进入缓冲段时间 1.4270 s 1.4375 s 0.74% 控制棒整体落棒时间 1.8587 s 1.9197 s 3.28%
下载: 导出CSV
-
[1] YOON K H, KIM J Y, LEE K H, et al. Control rod drop analysis by finite element method using fluid–structure interaction for a pressurized water reactor power plant[J]. Nuclear Engineering and Design, 2009, 239(10): 1857-1861. doi: 10.1016/j.nucengdes.2009.05.023 [2] PARHIZKARI H, AGHAIE M, ZOLFAGHARI A, et al. Analysis of control rod drop accident in PWRs with multipoint kinetics method[J]. Annals of Nuclear Energy, 2016, 88: 194-203. doi: 10.1016/j.anucene.2015.10.039 [3] KENNEDY G, LAMBERTS D, VAN TICHELEN K, et al. Experimental and numerical study of the MYRRHA control rod system dynamics[C]//Proceedings of 2017 International Congress on Advances in Nuclear Power Plants. Shimbashi, Minato: Atomic Energy Society, 2017: 24-28. [4] 张吉斌,高希龙,何航行,等. 控制棒落棒动力学数值计算[J]. 核动力工程,2020, 41(6): 218-223. [5] 周红,肖志,陶书生,等. 运行核电厂控制棒组件及其驱动机构异常事件的经验反馈[J]. 核安全,2013(1): 19-22,35, doi: 10.16432/j.cnki.1672-5360.2013.01.012. [6] 单秉昆. 核电站控制棒组件(RCCA)肿胀机理分析[J]. 科技视界,2014(35): 320-321. [7] 蔡家藩,陈增武,乔维,等. 核电站反应堆控制棒束组件检查技术应用[C]//中国核学会. 中国核科学技术进展报告——中国核学会2009年学术年会论文集(第一卷·第3册). 北京: 原子能出版社,2009: 463-470. [8] 周建钦,梁尚明,赵军. 高温下控制棒驱动机构结构完整性分析[J]. 机械,2018, 45(7): 13-17. [9] HALL M M. Stress state dependence of in-reactor creep and swelling[J]. Journal of Nuclear Materials, 2010, 396(1): 112-118. doi: 10.1016/j.jnucmat.2009.10.063 [10] NIKITIN E, FRIDMAN E. Extension of the reactor dynamics code DYN3D to SFR applications – Part I: thermal expansion models[J]. Annals of Nuclear Energy, 2018, 119: 382-389. doi: 10.1016/j.anucene.2018.05.015 [11] YUN S, MIKHAILOV G, KIM S J, et al. A study on the radial expansion reactivity of a metal-fueled sodium-cooled fast reactor via a physics experiment[J]. Annals of Nuclear Energy, 2021, 153: 108072. doi: 10.1016/j.anucene.2020.108072 [12] 孔珑. 工程流体力学[M]. 北京: 中国电力出版社,2001: 226-231. [13] 赵珂,陈昌义,席炎炎,等. 控制棒下落与流体流动的耦合状态方程及其保辛算法[J]. 应用数学和力学,2022, 43(9): 935-943. [14] 王瑞红,丁向东,肖林,等. 渗氢前后Zr-4合金板的力学性能及有限元模拟[J]. 中国有色金属学报,2002, 12(3): 544-549. [15] 薛淑娟,王云惠,赵文金,等. 锆、锆-4合金及新锆合金的热膨胀[J]. 核动力工程,2004, 25(3): 236-240. [16] 张莉丽. Zr-4合金堆外热物理性能的测量研究[J]. 中国原子能科学研究院年报,1990: 280-282. -
下载:
计量
- 文章访问数: 108
- HTML全文浏览量: 34
- PDF下载量: 25
- 被引次数: 0
