¹中国原子能科学研究院，北京

²哈尔滨工程大学核科学与技术学院，黑龙江 哈尔滨

³核工业学院，北京

收稿日期：2023年2月22日；录用日期：2023年2月28日；发布日期：2023年4月26日

摘要

铅铋堆中的液态铅铋金属流体与常规介质物性差异较大，其流致振动特性需要明确。针对反应堆中的圆柱横向流致振动，本文对单棒在不同铅铋来流工况下流致振动特性进行了数值研究。首先通过悬臂棒的理论固有频率验证了模态分析的正确性，通过模态分析获得单棒在空气和铅铋中的固有属性。然后采用流固耦合仿真计算开展不同工况下的铅铋流质中单棒流致振动特性研究。结果表明铅铋绕流固定单棒时漩涡脱落频率对应的斯特劳哈尔数为0.2，与经典单圆柱绕流斯特劳哈尔数一致。对于与流体有耦合的弹性棒，当入口流速和温度恒定时，受铅铋流体冲击，棒先发生脉动幅值较小的静变形，然后逐渐稳定地脱落漩涡，使得振幅逐渐增大。

关键词

铅铋介质，流致振动，漩涡脱落

Numerical Study on Vibration Characteristics of Single Rod Suffered Liquid Lead-Bismuth Cross-Flow

Dongdong Wang¹, Hai Zhao², Zhiwei Zhou¹, Duan Wang³, Lin Sun¹, Jiong Zhang¹, Hongyi Yang^1*

¹China Institute of Atomic Energy, Beijing

²College of Nuclear Science and Technology, Harbin Engineering University, Harbin Heilongjiang

³Institute of Nuclear Industry, Beijing

Received: Feb. 22^nd, 2023; accepted: Feb. 28^th, 2023; published: Apr. 26^th, 2023

ABSTRACT

The physical properties of liquid lead-bismuth metal fluid in lead-bismuth reactor are quite different from those of conventional medium, and its fluid-induced vibration characteristics need to be clarified. For the cross-flow induced vibration in reactor, the vibration characteristics of a single rod suffered different lead-bismuth cross-flow conditions are studied numerically in this paper. Firstly, the theoretical natural frequency of the cantilever rod is used to verify the correctness of the modal analysis, and the natural properties of the single rod in air and lead bismuth are obtained by modal analysis. Then the single rod fluid-induced vibration in lead-bismuth fluid under different working conditions was studied by fluid-structure coupling simulation. Results show that the Strouhal number corresponding to the vortex shedding frequency is 0.2, which is consistent with that of classical flow around the circular cylinder. For the elastic rod coupled with the fluid, when the inlet velocity and temperature are constant, the rod is impacted by the lead-bismuth fluid, and the static deformation with small fluctuation amplitude occurs first, and then the vibration amplitude increases gradually because of the steady vortex shedding.

Keywords:Lead Bismuth Medium, Flow Induced Vibration, Vortex Shedding

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

铅铋堆作为国际上最具潜力的第四代反应堆堆型之一，具有自然循环流量大，固有安全性好、热效率高的优点，但是由于铅铋流质密度大，流致振动问题较为突出。流致振动指固体与流体的物理场之间的相互作用所导致的流固耦合现象，流体力作用在固体上使之发生的形变会反过来影响流动特性，是一种复杂的非线性动态关系。因流致振动影响反应堆安全的事故并不少见，英国安格赛核电站由于锅炉管振动停工近3年时间后才得以恢复运营，日本文殊堆也发生过因流致振动导致二回路管道破裂致使钠泄漏的事故 [1] [2] [3] 。2012年美国San Onofre核电站因蒸汽发生器传热管流致振动导致部分传热管磨穿被迫永久关闭 [4] 。

大量学者对管束及棒束的流致振动进行了研究。Erfan Taheri [5] 等人研究了不同雷诺数和KC数下弹性支承圆柱体的振动问题，得到了圆柱体的振动速度随频率比的变化关系。Mingrui Zheng [6] 等人研究了绕流结构的涡激振动性能，发现结构的固有频率和截面形状引起的漩涡脱落位置是影响涡激振动响应的重要因素。Dolfen H. [7] 等采用LES模型模拟了带绕丝六边形燃料组件的铅铋流固耦合，结果预测了燃料棒之间的接触现象。Graham Kennedy [8] 等人完成了MYRRHA的一个全尺度的127棒模型实验，初步结果表明流致振动在带绕丝棒束中非常有限，并且在运行过程中不会引起燃料棒疲劳损伤。

上述研究多为水等常规介质，铅铋与常规介质物性存在较大差异，相关流致振动机理及振动参数的响应特性需进一步研究分析。综上，基于上述研究的不足，本文针对铅铋横向流在高温液态下悬臂棒的流致振动特性进行计算分析。

2. 流致振动仿真计算

由于燃料棒通道内或铅铋堆螺旋管蒸汽发生器管外部流动中均存在一定的横向流导致的涡激振动 [9] [10] ，故本研究采用悬臂式结构(即一端固定一端悬空)针对单棒结构模型在不同的流动工况下进行仿真计算，探究横向流流致振动现象的影响机理。计算的各个工况参数信息如下表1所示。

表1. 计算工况

针对模拟反应堆启动流动工况，棒的流致振动响应设置为入口流速随时间线性增加、入口温度随时间线性增加，获得结构振动响应随时间的变化；针对模拟铅铋流量和温度的脉动工况，棒的流致振动响应设置为入口流速随时间脉动、入口温度脉动，获得结构振动响应随时间的变化；针对模拟稳态运行的流动工况，获得结构振动响应随时间的变化关系。

2.1. 几何模型

建立单棒流致振动三维流体域与固体域模型，几何模型如图1所示，结构物性参数设置为：棒结构采用合金类材料，密度为21,000 kg/m³，杨氏模量为1.5 GPa，泊松比为0.41，结构阻尼比 $\zeta =0.0617$ 。采用该弹性模量以降低结构的固有频率，使结构更易振动。

图1. 几何结构

2.2. 网格划分与无关性验证

流体域和固体域的网格划分如图2、图3所示。流体域采用六面体结构化网格，固体域采用Solid 186单元。Solid 186是高阶3维20节点固体结构单元，具有中间节点，可以精确模拟本文单棒结构。

图2. 流体域网格

图3. 固体域网格

对所划分的网格进行无关性验证，图4为铅铋2 m/s绕流单棒的阻力系数及斯特劳哈尔数随网格数的变化，在网格数增加到约13万之后，阻力系数随网格数基本不变，但斯特劳哈尔数还在上升。在网格数增加到为约30万时，此时斯特劳哈尔数基本不变约为0.21，与经典单圆柱绕流中空气和水介质的斯特劳哈尔数一致。因此，最终采用32.2万的网格划分方案进行铅铋流场计算。

图4. 阻力系数Cd与斯特劳哈尔数Sr随网格数的变化

2.3. 计算求解

铅铋可视为不可压的牛顿流体，对于本研究的流致振动来说，仅考虑其流动，其流动控制方程如式所示：

$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_j{u}_i\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\mu {\nabla }^2{u}_i+\rho f_i$ (1)

针对流体域计算，采用Transition SST湍流模型以实现对边界层分离的精确预测，Transient SST四方程模型是由SST k-ω两方程湍流模型与关于间歇因子 $\gamma$ 和当地边界层动量厚度雷诺数 $\tilde{R}{e}_{\theta t}$ 的输运方程及相关经验公式构成。该方程相对于SST k-ω两方程湍流模型修正了产生项、耗散率项和和混合函数项。当地边界层动量厚度雷诺数 $\tilde{R}{e}_{\theta t}$ 影响到SST k-ω模型的修正，而间歇因子 $\gamma$ 指出了边界层中层流、过渡层和湍流的位置,即层流区 $\gamma =0$ ，湍流区 $\gamma =1$ ，过渡区 $\gamma$ 取值为0到1之间。其中关于流动间歇因子 $\gamma$ 和当地边界层动量厚度雷诺数 $\tilde{R}{e}_{\theta t}$ 的输运方程如式(2)所示：

$\{\begin{array}{l}\frac{\partial \left(\rho \gamma \right)}{\partial t}+\nabla \cdot \left(\rho U\gamma \right)=\nabla \cdot \left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\gamma }}\right)\nabla \gamma \right)+P_{\gamma }-D_{\gamma }\\ \frac{\partial \left(\rho \tilde{R}{e}_{\theta ,t}\right)}{\partial t}+\nabla \cdot \left(\rho U\tilde{R}{e}_{\theta ,t}\right)=\nabla \cdot \left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\theta ,t}}\right)\nabla \tilde{R}{e}_{\theta ,t}\right]+P_{\theta ,t}\end{array}$ (2)

式中 $P_{\gamma }$ 、 $P_{\theta ,t}$ 为产生项， $D_{\gamma }$ 为扩散项。

采用SIPLEC算法求解流场。流场动网格的运动由扩散方程控制：

$\nabla \cdot \left(\gamma \nabla \stackrel{\to }{u}\right)=0$ (3)

式中 $\stackrel{\to }{u}$ 是网格位移速度。在变形边界上，边界条件是网格运动与边界相切(即法向速度分量消失)。方程(3)描述了如何利用扩散系数 $\gamma$ 控制规定的边界运动扩散影响到变形网格内部运动，常数意味着边界运动在整个网格中均匀扩散。对于非均匀扩散系数，高扩散率区域中的网格节点倾向于一起移动(即相对运动较少)。

针对固体域计算，对于固体部分的结构动力学方程，使用有限元方法将其离散为有限元半离散方程，如下：

$M\left\{\ddot{u}\left(t\right)\right\}+C\left\{\stackrel{\dot{}}{u}\left(t\right)\right\}+F^i\left(t\right)=F^a\left(t\right)$ (4)

式中，M为结构质量矩阵，C为结构阻尼矩阵， $\left\{\ddot{u}\left(t\right)\right\}$ 为节点加速度矢量， $\left\{\stackrel{\dot{}}{u}\left(t\right)\right\}$ 为节点速度矢量， $\left\{u\left(t\right)\right\}$ 为节点位移矢量， $F^i\left(t\right)$ 为内部载荷矢量， $F^a\left(t\right)$ 为外部载荷矢量。利用Newmark迭代求解上述方程组。通过将结构物性参数的输入方程，可实现对不同结构物性参数下的瞬态结构动力学计算。结构域采用瞬态动力学模块求解。

流固耦合交界面运动位移与力满足如下方程：

$n_l\cdot {\psi }_l=n_p\cdot {\psi }_p，d_l=d_p$ (5)

式中，n_l、n_p分别为流体域和结构域交界面的方向矢量；ψ_l，ψ_p分别为流体域和结构域交界面上的应力；d_l，d_p为流体域和结构域交界面上的位移。

3. 计算结果分析

3.1. 计算结果模态分析

固有频率的确定为流致振动计算时的时间步长确定和结构动力学响应特性分析奠定基础。同时，当铅铋绕流棒时的激励频率等于该固有频率时会发生共振，所以固有频率也为绕流流速的选定提供依据。对悬臂梁结构进行模态分析以确定棒在真空中(干模态)、铅铋液体中(湿模态)的固有频率，由于悬臂棒关于z轴对称，在x、y方向自由度上的频率相等，因此每两阶的频率相等。下面仅展示奇数阶的振型特征，偶数阶的振型和频率与此相同。理论上，悬臂梁结构的固有频率计算式如式所示：

$f_n=\frac{1}{2\pi }{\left({\beta }_{n}l\right)}^2\sqrt{\frac{EI}{\rho S{l}^4}},\left(n=1,2,\cdots \right)$ (6)

式中n为模态阶数，当n = 1，2，3时， ${\beta }_{1}l=1.875$ ， ${\beta }_{3}l=7.855$ ，E为杨氏模量，单位为Pa，I为截面惯量，单位为m⁴，ρ为密度，单位为kg/m³，S为横截面积单位为m²，l为长度，单位为m。

由仿真计算得到本研究结构的干模态第一阶固有频率为16.63 Hz，振型图如图5所示。计算结果与式(6)理论计算的一阶模态频率结果误差为0.06%，说明结构有限元模态计算结果可靠。图6为棒的第三阶振动振型，其第二阶振型为二阶弯曲振型，对应的频率为97.7 Hz。

图5. 干模态第一阶振型

图6. 干模态第三阶振型

棒在铅铋中的第一阶振动振型如图7所示，可知其第一阶振型为一阶弯曲振型，对应的频率为13.8 Hz。图8、图9分别为棒在铅铋中的第三阶振动振型和第五阶振动振型，其第三阶振型为二阶弯曲振型，第五阶振型为三阶弯曲振型，对应的频率分别为85.7 Hz、221.6 Hz。

图7. 棒在铅铋中的第一阶振型

图8. 棒在铅铋中的第三阶振型

图9. 棒在铅铋中的第五阶振型

3.2. 流致振动分析

工况1.1棒自由端的振动位移如图10所示，应变如图11所示。初始时刻的振动响应具有一定的脉动，后期的脉动不明显，而静变形较明显。这是因为初始时刻的流速较低，结构受力静变形相对较小，随着流速的增大，结构静变形增大，相对振动振幅较小。随流速的线性增加，流向(x方向)位移呈逐渐增大的变化，而横向(y方向)和轴向(z方向)的位移变化不明显，流向位移明显大于横向和轴向的位移变化。这是因为棒迎流受流速不断增大的铅铋流体冲击，且铅铋的温度呈线性增加，与恒定流速相比，脱落的涡不规律使得涡激励效应较小，所以横向脉动幅值较小。而轴向刚度更大，所以轴向变化同样小。总体来说振动位移偏向一侧呈增大的趋势。图11中应变变化具有相似的趋势。

针对工况1.2，棒自由端的振动位移如图12所示。棒受流速和温度线性增大的铅铋流体作用，流向位移变化趋势与工况1.1一样，均不断地增大，而横向振动位移在超过3 s后具有一定幅值的脉动，且脉动幅值不断增大，表明在超过3 s后棒后漩涡脱落效应开始出现，但漩涡脱落的流向激励幅值与流向静变形相比较小，因此导致流向的脉动幅值不明显。

图10. 工况1.1棒自由端振动位移随时间变化

图11. 工况1.1棒中点应变随时间变化

图12. 工况1.2棒自由端振动位移随时间变化

工况二为模拟铅铋流量和温度的脉动，探究棒的流致振动响应，此时入口流速、入口温度均随时间脉动，棒自由端的振动位移如图13所示，棒流向、轴向位移随时间变化均具有一定的脉动，横向虽有变化但并不明显。

图13. 工况2棒自由端振动位移随时间变化

工况三为模拟稳态运行，流场参数保持稳定。针对工况3.1，棒中点应变的变化如图14所示。从应力、应变随时间变化看，流向振动周期约为14 Hz，横向振动周期为流向振动周期一半，对应的斯特劳哈尔数约为0.2，即此时结构振动频率与流场激励频率一致。同时该频率也与3.1节中棒在铅铋内的湿模态频率13.8 Hz相近，表明该工况位于共振区。随时间增加，振动幅值逐渐降低，这是由于棒变形到静平衡位置。所以可以预测在后面的时间内棒的振动幅值会进一步增加。从图15可以看出，横向振动幅值在约t = 1.5 s后有增加的趋势。

图14. 工况3.1棒中点应变随时间变化

图15. 工况3.1棒自由端振动位移随时间变化

针对工况3.2，棒自由端的振动位移变化如图16所示。从图中可以看到随着时间的推移，棒的振幅先减小，后逐渐增加到基本稳定振动。这是由于棒先到达静变形最大的位置处，然后涡开始稳定地脱落，这时横向振幅开始增大且随铅铋湍流和涡激励而振动，流向和轴向振动周期约为横向振动周期的2倍。图17为棒中点应力变化，其中流向应力具有与流向位移相同的变化趋势，这是因为悬臂梁棒结构的位移与应力存在一定的线性换算关系。

针对工况3.3，棒自由端位移随时间变化如图18所示，同样随时间增加，棒的振动振幅先减小，后增大，后期增大是因为铅铋绕流棒的尾流具有较稳定的漩涡脱落激振。

图16. 工况3.2棒自由端振动位移随时间变化

图17. 工况3.2棒中点应力变化

图18. 工况3.3自由端振动位移随时间变化

4. 结论

本文为探究铅铋横向流在高温液态下流致悬臂棒的振动特性，采用计算流体动力学与瞬态结构动力学进行单棒流致振动双向流固耦合三维仿真计算，获得悬臂棒结构在不同铅铋横向流激励条件下的振动响应数据，基于不同边界条件的流致振动响应分析，得到主要结论如下：

1) 由于铅铋附加质量效应，棒在铅铋中的固有频率小于在真空中的固有频率；

2) 悬臂棒在铅铋中的湿模态第一阶模态振型与在真空中的干模态第一阶振型一样，均为一阶弯曲振型；

3) 当入口流速和温度随时间线性增大时，初始一段时间内棒的位移变化以静变形为主，振幅不明显。后期铅铋绕流棒的尾流存在漩涡脱落，使得棒横向和流向振幅逐渐增大，而流向振幅相对于流向的静变形较小。当入口流速和温度脉动时，棒的流向存在一定的脉动效应，而横向不明显。当入口流速和温度恒定时，受铅铋流体冲击，棒先发生脉动幅值较小的静变形，然后逐渐稳定地脱落漩涡，脉动振幅逐渐增大。

本文主要分析了典型工况下铅铋绕流单棒及其流致振动特性，未来的工作需要详细对比研究铅铋与其他介质流致振动的差异，并增加模拟的棒数，拓展模拟工况致流弹失稳等流致振动机理区域。

文章引用

王东东,赵 海,周志伟,王 端,孙 林,张 泂,杨红义. 液态铅铋横向流作用下单棒流致振动特性数值研究
Numerical Study on Vibration Characteristics of Single Rod Suffered Liquid Lead-Bismuth Cross-Flow[J]. 核科学与技术, 2023, 11(02): 155-167. https://doi.org/10.12677/NST.2023.112017

参考文献