控制棒下落与流体流动的耦合状态方程及其保辛算法*

赵 珂，陈昌义，席炎炎，黄东威，吴 锋，钟万勰

（1.大连理工大学 工程力学系 工业装备结构分析国家重点实验室，辽宁 大连 116024；2.中广核研究院有限公司，广东 深圳 518000）

（我刊编委钟万勰来稿）

引 言

控制棒是核反应堆在发生紧急情况时控制停堆的重要部件，控制棒的下落时间是保证核反应安全的关键指标.近年来，国内外许多学者对控制棒下落时间的问题做了大量研究：文献[1]通过商业软件ADINA 建立了单根控制棒二维流-固仿真模型，并将ADINA 得到的计算结果与未公开的企业内部代码得到的结果进行比较，落棒时间误差在7%以内；文献[2]基于FLUENT 动网格算法，建立了单根控制棒三维流体仿真模型，将得到的控制棒落棒时间与试验结果进行比较，误差小于6%.采用二维、三维动力学模型的缺点是自由度多、计算量大，因此计算效率较低，往往只能分析单根控制棒，不能分析整个控制棒组件.考虑到控制棒和导向管均具有长细比大的特点，且控制棒的直径与导向管的直径相近，文献[3]基于一维Navier-Stokes 方程建立了控制棒下落的数学模型，得到的结果能够很好地匹配试验结果；文献[4-5]构建了基于一维Navier-Stokes 方程的落棒动力学模型，通过编译C++程序进行计算.这些研究中常将流体与固体分开进行分析[4-5]：首先，假定控制棒的运动状态恒定不变，计算流体流动的流量和压强分布，得到流体阻力；然后，假定流体为拟稳态流动，其速度和压强状态不变，根据流体阻力计算控制棒的速度等变量.这种算法本质上是动力学求解算法中的算子分裂法[6]，因此只有一阶精度，计算精度较低，需要取非常小的时间步长[3]才能保证“拟稳态流动”近似的精度.此外，控制棒下落过程中，流场中的流速、压强等会发生非线性突变，采用单一的时间步长，算子分裂算法还存在不收敛问题，本文的研究将致力于解决这些问题.

在控制棒下落过程中，控制棒的下落与流体流动是同时发生的，因此应该采用流固耦合方法，以真实地描述控制棒的下落过程.基于流固耦合思想，本文首先建立了控制棒运动方程和导向管内流体压降平衡的非线性微分方程组；然后引入关键恒等变换=v，得到非线性微分方程组的状态方程.为精确求解该非线性状态方程，本文建立了一种时间步长自适应保辛算法.该方法可以根据流动的实际状态实时调整时间步长，精确捕捉流动状态的突变，避免数值发散，并可以在流动变化平稳的时间段采用较大的时间步长，提高计算效率.

1 理论模型

1.1 物理模型

控制棒和导向管的纵剖面图如图1所示，其中阴影部分表示控制棒，控制棒下面是导向管，控制棒的直径小于导向管的直径，且控制棒的长度大于导向管的长度.导向管和控制棒全部竖直放置，导向管管壁和底部中心处均有流水孔，分别记为侧孔和底部流水孔.导向管和控制棒完全浸在水中，底部流水孔保持进水状态，而侧孔存在流进和流出两种状态的突变.在控制棒下落过程中，导向管固定不动.初始时刻，控制棒在导向管侧孔上方，如图1(a)所示.控制棒由静止释放后，在重力、流体作用力等不同力的作用下开始下落，经过一段时间，控制棒底部下落至导向管侧孔下方，如图1(b)所示.

对图1中的结构建立坐标.以导向管底部截面圆心为坐标原点o，以导向管中心轴为z轴，向上为正.图1中，导向管底部Z0的坐标为0，对应的压强为PZ0；导向管侧孔的坐标为Z1，对应的压强为PZ1；控制棒底部的坐标为Z，对应的压强为PZ，Z在下落过程中会随时间而变化；导向管顶部的坐标为Z2，对应的压强为PZ2.

图1(a)中，控制棒底部在导向管侧孔上方.根据导向管内流量的不同可将导向管分成三个区间，分别为[Z,Z2]，[Z1,Z]，[Z0,Z1]，对应的流量分别记为Q2，Q1_Z，Q0.从导向管侧孔流进或流出的流量记为Q1，以流出为正.同样地，图1(b)中也根据导向管内流量的不同，将导向管分成三个区间，分别为[Z1,Z2]，[Z,Z1]，[Z0,Z]，对应的流量记为Q2，Q1_Z，Q0，而导向管侧孔处的流量仍记为Q1，以流出为正.

图1 控制棒在导向管内的下落示意图Fig.1 The falling diagram of the control rod in the guide tube

为方便论述，导向管顶部以上的控制棒部分所受的力暂不考虑；且假定控制棒在下落过程中保持竖直状态，即控制棒没有横向偏移和横向振动；反应堆内的流体满足不可压缩条件.

1.2 控制棒的运动方程

控制棒的运动受Newton 第二定律控制，其控制方程可以写成

其中，M是控制棒的质量；w是控制棒的位移；表示控制棒的加速度；F是控制棒受到的合力，可表示为

式(2)中，Fg为控制棒的重力；Fb为控制棒受到的浮力；frod为控制棒受到的流体摩擦阻力；Fr为控制棒受到的流体压差阻力；Fc为控制棒受到的机械摩擦力，由实验拟合得到；Fs为控制棒受到的弹簧阻力，此力在控制棒快要到导向管底部时产生，防止由于控制棒速度过大从而与导向管产生强烈碰撞，可表示为

其中，Fprl表示弹簧预紧力，Ks表示弹簧刚度，hs表示控制棒接触弹簧时的位移.

在某区间段[Za,Zb]内，摩擦力frod可表示为[7]

其中，ρ为流体密度，vf为流体流速，v为控制棒的下落速度，C2为控制棒的摩擦因数，l2为控制棒的湿周长，sgn为符号函数.

在式(2)所示控制棒受到的几种合力中，重力和浮力很容易计算，摩擦阻力可通过式(4)来计算，而压差阻力Fr的计算涉及到流体的非定常流动，计算相对复杂，将在下一小节详细阐述.

1.3 压降和压差阻力计算

导向管内的流体流动可分为两种情况，一种为图1所示Q0对应的区间，流体在变截面的圆形导向管内流动；第二种为图1所示Q2对应的区间，流体在变截面的导向管和控制棒之间的环形通道内流动.由于通道的长细比很大，流动可模拟为一维N-S 方程[1]：

其中，t为时间，P为压强，C1为导向管的摩擦因数，l1为导向管的湿周长，S1为控制棒与导向管之间的环形面积，Kc为局部形阻系数，由实验测定，令环形面积之间的流量为Q=vfS1.将上式在区间 [Za,Zb]上积分，可得压降方程：

引入如下变量：

则压降方程(6)化简为

式(7)中的I1～I5均有特定的物理意义：I1表示由于局部加速度导致的压降变化系数，称为惯性压降系数；I2表示对流加速度产生的压降系数，称为对流压降系数；I3表示导向管受到的流体摩擦力产生的压降；I4表示控制棒受到的流体摩擦力产生的压降；I5表示由于导向管形状改变产生的形阻压降.根据求得的压降，可进一步求出控制棒在区间[Za,Zb]受到的压差阻力[5]：

其中，Srod表示控制棒的横截面积.

以上分析讨论的是导向管内存在控制棒时的压降和压差阻力计算.当导向管内没有控制棒时，式(8)压降的计算要减去控制棒摩擦力产生的压降部分，即减去I4(Za,Zb)，则压降方程为

2 控制棒与流体耦合运动状态方程

当控制棒底部在导向管侧孔上部或者下部时，流场分布是不同的，压降方程的建立以及控制棒所受压差阻力的计算也是不同的，因此需分别讨论.

2.1 控制棒底部在导向管侧孔之上

首先讨论控制棒底部在导向管侧孔上部的情况.如1.1 小节论述，当控制棒底部在导向管侧孔之上时，将导向管分成三个区间，分别为[Z,Z2]，[Z1,Z]，[Z0,Z1].在区间 [Z,Z2]内，由于流体在圆环内流动，其两端压降计算采用式(8)，可表示为

式(11)中的最后一项表示由于局部损失产生的压降.为方便表示，记I11=I1(Z,Z2)，I11的第一个下标表示压降类型，此处为惯性压降系数；第二个下标表示导向管的区间，此处为区间[Z,Z2].其他压降类型和其他导向管区间也类似表示，所以上式变为

在[Z1,Z]和 [Z0,Z1]两段，流体在圆管内流动，因此可利用式(10)计算两端压强，分别表示为

式(14)中最后一项表示流体从导向管底部流水孔流进导向管时产生的局部突变压降，其中表示局部压降损失系数，SK表示导向管底部流水孔的面积.在导向管侧孔内外存在局部突变的压降[4]可表示为

其中，C3为局部压降损失系数[7]；ζC3是开孔的局部水头损失系数，一般取1；nC3为导向管侧孔的数目；SC3为侧孔的面积.根据式(12)、(13)、(15)可得

根据式(14)、(15)可得

式(16)中的PZ1-PZ2表示导向管侧孔外与导向管顶部外侧之间的压降；式(17)中的PZ0-PZ1为导向管底部外侧与导向管侧孔外之间的压降.这两个压降的计算通常通过实验求得.

式(16)与式(17)可以组成一个非线性方程组，其中的未知量有Q0，Q1，Q2.观察图1(a)中的四个流量及控制棒的运动速度v，根据不可压缩条件，可以得出如下关系：

式(16)、(17)与式(1)可以组成三个方程和三个未知量的非线性微分方程组.此时可以使用将控制棒的运动和流体的流动解耦的方法求解.首先，假定控制棒运动状态不变，计算流体流动的流量和压强，从而得到流体阻力；然后，假定流体的流动状态不变，根据流体阻力计算控制棒的速度和位移.计算控制棒的运动状态时，文献[3]采用了向前差分的方法，文献[4-5]采用了Taylor 展开的方法.然而，这些传统方法违背了控制棒下落运动的真实情况，只能通过不断缩小迭代时间步长来提高精度，从而导致计算效率较低.

本文将控制棒的运动和流体的流动同时分析，即通过流-固耦合对控制棒下落进行求解.此时引入关键恒等变换=v，将微分方程的个数升级为四个，非线性微分方程组变为

其中

将式(19)写成状态方程的形式：

式(21)即为描述控制棒运动与导向管内流体流动的流-固耦合状态方程，其中U为状态向量，是待求解的未知量.当给定初始时刻的状态向量，即可直接对上述非线性微分状态方程进行时程积分，得到各个时间节点的状态向量，同时得到控制棒的位移、速度、流体的流量等所有未知量.对于该方程有许多成熟算法，其中保辛算法[8-11]具有精度高和稳定性好的优点，在解决动力学[12]问题上有独特的优势，并在动力学的最优控制问题[13-14]、水波问题[15-16]等领域得到广泛运用.本论文选择辛Euler 中点格式，此部分将在后续介绍.

2.2 控制棒底部在导向管侧孔之下

同理，如1.1 小节中所述，当控制棒底部在导向管侧孔下方时，将导向管分成三个区间，分别为[Z1,Z2]，[Z,Z1]，[Z0,Z].

类似于上一小节的方法，可以得到如下非线性微分方程组：

其中

本小节与上一小节的不同主要体现在计算导向管内压降和控制棒所受压差阻力上.本小节中控制棒所在区间为[Z1,Z2]，[Z,Z1]，在计算区间[Z,Z1]内导向管的压降时，需考虑控制棒所受摩擦阻力产生的压降；同时，本小节中控制棒受到的压差阻力需将这两个区间上控制棒受到的压差阻力相加.

同理，将式(24)写成与式(21)完全相同的状态方程，其中f1，f2和f3的表达式见式(25)，而H(U)的表达式如下：

至此，我们给出了描述控制棒下落与流体流动统一的耦合状态方程，无论控制棒在导向管侧孔上方还是下方，耦合方程的形式是相同的，即为f(U)=H(U)，不同之处在于H(U)和f(U)的表达形式.如果考虑控制棒上部分驱动杆结构和流场分析，也可表示为如式(21)所示的统一的状态方程求解.

3 时间步长自适应保辛算法

求解式(21)所示状态方程，本文选择辛Euler 中点格式，在此详细介绍.假设状态方程(21)中t0时刻U0已知，要计算U1.将状态方程变形为

则根据辛Euler 中点格式[9]有

其中

式(28)是一组非线性微分方程，可以迭代求解，首先假设=U0，然后按n=0,1,2,··· 开始循环：

直到与之间的相对误差小于允许值，则停止计算，此时U1=.再以U1为初值，计算U2，U3,···.计算到临界边界条件的值时，结束计算.

在整个落棒过程中，导向管内的流场是非定常的，甚至出现突变，因此状态向量的时程曲线并非光滑变化，将导致数值计算的不稳定.具体体现为迭代格式(30)在实际计算时不收敛.一般需要将时间步长取得足够小才能保证数值计算收敛.然而在大部分时程积分区间，状态向量的时程曲线是光滑的，时间步长并不需要取很小就可以保证收敛.综合以上两方面考虑，在保证整个时程分析收敛的同时，尽可能地提高计算效率，本文引入时间步长自适应算法.该方法将采用一个稳定的、较大的时间步长，记为Δt0（比如0.01 s）.因为在大部分时间段内，流场没有突变，状态向量变化是光滑的，因此取较大的时间步长 Δt0，可以保证迭代的收敛.当流场发生突变时，仍然取 Δt0会出现迭代格式不收敛的情况，据此可以根据迭代步数来判断流场是否发生突变.在实际计算时，以给定的允许迭代步数Ni为依据，如果某次计算时，迭代步数大于Ni，则判断该时间步长内，流场发生了突变，这时将时间步长减半为0.5Δt0，重新迭代.重新迭代后，如果迭代步数仍然超过Ni，则将时间步长再次减半为(0.5)2Δt0，重新迭代；如果迭代步数没有超过Ni就收敛，则结束本次迭代，进入下个时间步计算.当给定Δt=Δt0时，在每个时间步的具体计算流程为：

1）令=U0，n=0.

2）判断n＜Ni是否成立，如果是，计算式(30)，得到，进入步骤3）；如果否，Δt=进入步骤1）.

4）U1=，进入下一个时间步进行计算.

4 算 例

以本文建立的落棒分析模型分析某小型反应堆的控制棒下落过程，并将计算结果与在核工业长期使用的某商业软件计算结果进行对比.整个控制棒下落组件分为上下两个部分，这两部分由星型架连接.其中下部分主要是24 根控制棒插入到24 根导向管中，上部分是驱动杆和驱动机构，驱动机构包括：热套管、调节器、钩爪、行程套筒等.小堆的主要输入参数如表1所示.

表1 小堆的主要输入参数Table 1 Main input parameters for a small reactor

图2所示为采用本文模型和商业软件计算得到的控制棒在下落过程中的位移、速度和加速度的时程曲线对比图，其中实线为本文方法的计算结果，虚线为商业软件的计算结果.本文模型采用自适应时间步长计算，时间步长为Δt0=0.01 s ；商业软件的时间步长为0.001 s.从图2可以看出，本文模型采用时间步长 0.01 s计算得到的位移、速度和加速度时程曲线与商业软件采用时间步长为0.001 s计算得到的位移、速度和加速度时程曲线高度吻合.计算结果不仅验证了本文建立的状态方程的正确性，且说明了本文方法的高效性.

图2 小型反应堆下本文模型与商业软件关于位移、速度、加速度随时间的变化对比Fig.2 The comparison of time-varying displacements,velocities and accelerations,between the proposed model and the commercial software for the case of a small reactor

以小堆为例，本文模型的初始时间步长 Δt0分别取为0.001 s，0.005 s，0.01 s，0.015 s，商业软件的初始时间步长取为0.001 s.将本文模型得到两个关键参数T5和T5+T6（T5为进入缓冲段时间，T5+T6为控制棒整体落棒时间）与商业软件得到T5和T5+T6进行比较，结果如表2所示.

表2中，括号内的百分数表示本文模型的计算结果与商业软件的计算结果的相对误差.从表中可以看到，本文使用的4个时间步长所得结果与商业软件对比，相对误差均小于1%.本文使用的最大时间步长Δt0=0.015 s，是商业软件时间步长的15 倍.而当商业软件的时间步长取为0.002 s 时，计算结果发散.

表2 本文模型在不同时间步长下与商业软件关于T5 和T5+T6的比较Table 2 The comparison of T5 and T5+T6 between the commercial software and the proposed model for different initial time steps

5 总 结

针对反应堆内控制棒下落问题，本文建立了导向管内控制棒下落和流体流动的理论模型，重点分析了控制棒受到的摩擦阻力、压降阻力等作用力.在此基础上，建立了描述控制棒下落与流体流动的耦合非线性方程组，再通过引入关键恒等变换v=，得到了流固耦合的非线性状态方程，进一步通过高精度的辛Euler 中点格式对得到的状态方程进行了时程分析.考虑到在不同时程分析段内导向管中流动状态的突变，本文将时间步长自适应算法引入到时程分析中.最后以某小型反应堆控制棒下落为例进行模型和算法的验证.算例中本文模型选取的最大时间步长 Δt0是商业软件时间步长的15 倍，关于落棒时间T5和T5+T6的相对误差小于1%，从而证明了本文模型和算法的正确性和高效性.本文建立的流固耦合非线性状态方程和时间步长自适应算法不仅为控制棒下落分析提供更符合实际的计算思路，还为更加复杂的圆管内流动模型或者圆环内流动模型提供了高精度、高效率的解决方案.

致谢本文作者衷心感谢大连市“青年科技之星”项目（2018RQ06）对本文的资助.