冯志鹏，臧峰刚，张毅雄

(中国核动力研究设计院 核反应堆系统设计技术重点实验室，四川 成都 610041)

在核反应堆结构中，圆柱和圆柱束结构是流体诱发振动的主要部件，例如在横流作用下的蒸汽发生器传热管、燃料棒等。根据国内外核电站中对蒸汽发生器部件失效的多例事故的统计分析可知，传热管的流致振动和与其相关的冲击磨损、微动接触疲劳加上介质腐蚀使得管壁逐渐变薄，是导致传热管承压能力降低而破裂的主要原因。因此，对其流体诱导振动特性进行深入研究是十分必要的。

传热管的流体诱导振动问题得到了许多学者的关注[1]，他们在很大雷诺数范围内对圆柱体的流致振动进行了一系列的研究，包括用有限元方法[2]或计算流体力学(CFD)方法[3]。Gabbai等[4]对过去二十几年中，经受稳态均匀流的大多数圆柱结构的流致振动所取得的研究进展进行了全面综述，详细讨论了目前流体诱导振动问题研究现状的要点与不足。随着工业设备的过程参数越来越高，流体与结构之间存在着强烈的交互作用[5]，因此，需要构建更精确的物理模型来分析结构和流体的相互作用。而现有的研究对象大多是二维弹性支撑刚性柱体，不能反映真实的三维流动状态，且将弹性管简化成弹性支撑的刚性管，忽略了管自身的弹性，不能考虑结构的弹性变形与流体流动的相互影响。

本文基于有限体积法和有限元法，利用动网格控制技术模拟结构的运动，联合CFD和计算结构动力学(CSD)建立传热管的流固耦合振动模型，旨在通过数值方法来研究传热管和流体间的相互作用，为设计者提供可靠、合理的振动分析，为自主设计蒸汽发生器传热管束的流致振动分析提供技术支持。

1 数值模型

对于流体诱导振动问题，任何数值模拟都需考虑以下4个基本问题：流场模拟、结构振动模拟、流体-结构交互作用模拟、数据处理。对这样一个系统来说，控制方程包括流体的质量和动量守恒及描述结构运动的动力学方程。因此，流体诱发振动问题的研究涉及到流体控制方程与结构动力的双向耦合，近年来已有许多关于流体诱发振动的综述性文献[6]报道。

本文采用有限体积法离散三维、黏性、瞬态、不可压缩Navier-Stokes方程，并联合大涡模拟方法求解流场区域，利用有限元方法离散传热管结构体，采用Newmark积分方法求解瞬态动力平衡方程来获得结构的位移、速度等响应。考虑结构大变形以及由大变形带来的流场网格的变形问题，采用基于扩散光顺的Diffusion方法控制运动边界的网格更新，通过流固耦合交界面进行固体域和流体域间的数据传递，文献[7]对解耦和求解过程进行了详细说明，本文不再赘叙。

2 物理模型及计算参数

传热管的物理参数：管长0.5 m，外径D=0.01 m，内径Di=0.009 5 m，弹性模量E=10 GPa，泊松比υ=0.3，密度ρs=6 500 kg/m3，阻尼比ζ=0.047。

流体参数：流体为水，密度ρ=998.2 kg/m3，动力黏度μ=1.003 mPa·s，无量纲入口流速Ur=U/fnD=0.5～17，其中U为来流速度(m/s)，fn为管的固有频率(Hz)。

流场区域及边界条件：流场计算区域及网格如图1所示。选用ICEM CFD作为网格划分工具，采用六面体结构化网格。图1中左边入口采用速度入口边界条件，右端出口采用压力出口边界条件，其他外边界按对称边界和固定壁面处理；管壁面为流固交界面，设为动网格边界。

图1 流场区域及局部网格示意图

时间参数：结构计算与流体计算的时间步长取为0.000 25 s。

3 结果分析

3.1 响应分支

图2 振幅随Ur的变化

图3 Griffin图

图4 升力、阻力、振幅随Ur的变化

在下分支阶段Ⅳ(5

在非相干阶段Ⅴ(Ur≥9)，流动速度的进一步增加导致Ay/D与ClRMS的下落，两者均趋于零，由此也可看出漩涡强度是自限定的；阻力系数维持在一个比较恒定的值。

从图4还可看到，最大的阻力系数和升力系数并不是同时出现的，升力峰值先于阻力峰值。ClRMS曲线的显著特征是存在1个尖峰，其恰好出现在初始分支与拟上端分支的转换处，在转换前，升力急剧增加，而当响应模式发生转换后，升力迅速下落，因此就导致了这样1个尖峰。

3.2 频率锁定

图5示出传热管的旋涡脱落频率fvs和振动响应频率fex与固有频率fn之比fvs/fn、fex/fn随Ur的变化情况。图5b中对角线代表静止管的旋涡脱落频率，即St=0.2，水平线代表固有频率，即fex=fn，实点指功率谱中的主频峰值。

图5 旋涡脱落频率和响应频率随Ur的变化

响应频率fex随Ur的变化规律为：在初始分支，升力与横向位移的频率与静止管的旋涡脱落频率fst相同，遵循St=0.2的线性关系；当进入拟上端分支和下端分支阶段时，旋涡脱落被升力方向的运动所控制，从图5可看出，旋涡脱落频率与振动响应频率偏离了St=0.2这条线，基本上等于管的固有频率，即fex/fn≈1.0，即发生了“频率锁定”，同步区域为3

3.3 相位差

传热管的流致振动特性还可通过升力与位移间的相位角来表征有关响应信息，相位角φ定义为升力相对于横向位移的相位差，φ与时间的函数可通过Hilbert变换得到。Lissajou图表示了传热管运动与流体间的能量传递，其倾角则给出了位移与升力间相位角的估计。

图6 相位角φ随Ur的变化

图6为相位角φ随Ur的变化情况。可明显看到，当Ur=9～13时，相位差经历了从同相到反相的改变，相位差的这种跳跃称为“相位开关”，是一种典型的非线性现象。图7示出了各Ur下的相位角时程，通过综合分析各响应阶段的相位角时程，发现当Ur<9时，相位角时程的变化趋势类似，相位角φ接近于0°，其时程也非常恒定(如Ur=1.75)；在相位差的转变阶段(如Ur=10)，相位角φ的时程开始变得紊乱，也不再接近于0°；而当流速超出“频率锁定”的临界速度，处于高Ur时，相位角φ在360°范围内变化，时程变得非常紊乱(如Ur=16)。

图7 不同Ur下的相位差时程

图8为不同Ur下管的运动轨迹，图9为不同Ur下传热管的Lissajou图。从图8、9可看出，当Ur≤3时(图9a)，相位差接近于0，升力与位移同相；随着流速的增大，平衡位置向下游移动，横向变形逐渐增大，管的运动轨迹是阻力方向和升力方向频率比为2的“8”字形图，当Ur=6时达到最大，相应的流体力也明显增大，Lissajou图为类似椭圆的平行四边形(图9b)；在非相干阶段(Ur≥9)，运动轨迹变得紊乱，不再是“8”字形，升力系数与位移间的相位也由同相变为反相(图9c)，升力、振幅随流速的增大而逐渐减小。

图8 不同Ur下的运动轨迹

3.4 相图

图9 不同Ur下的Lissajou图

图10 不同Ur下位移的相图

图11 不同Ur下升力的相图

图12 不同Ur下位移的Poincare截面

图13 不同Ur下升力的Poincare截面

当2

Lissajou图、运动轨迹、相图、Poincare截面映射都表现了三维弹性管流致振动的周期性，对所有流速范围内各种响应的分析发现，根据不同阶段的响应分支特征，Lissajou图、运动轨迹、相图、Poincare截面映射均呈现不同的形式，联合这4种图与响应分支图对三维弹性管的流致振动进行分析，便可较容易地确定其运动行为及振幅响应特性。

4 结论

本文基于双向流固耦合方法，建立了三维横向流体诱发传热管振动的数值模型，对不同流速下单管的流致振动特性进行了数值计算，得出以下主要结论：

1) 高m*ζ低m*流体-弹性管耦合系统的响应与高m*ζ高m*系统和低m*ζ低m*系统相比，会出现拟上端分支，振幅响应具有更宽的同步区域。

2) 当Ur≤2时，管的振幅很小，流致振动机理为湍流激励；升力系数均方根随Ur的增加而增加，阻力系数均方根随Ur的增加先增大后减小，在Ur约为2时达到最小值；相位角φ接近于0°。

3) 当3

4) 当Ur≥9时，升力系数均方根趋于0，横向振幅随流动速度的增加而下落；当Ur=9～13时，相位差经历了从同相到反相的改变，相位角φ在360°范围内变化，其时程变得非常紊乱。

5) 在均匀湍流流动作用下，三维弹性管的升力与横向位移在Ur=0.5～18的范围内未出现周期解的分叉。

参考文献：

[1] 高李霞，蒋自龙，马建中，等. 螺旋式传热管束流致振动试验研究[J]. 原子能科学技术，2008，42(增刊)：468-471.

GAO Lixia, JIANG Zilong, MA Jianzhong, et al. Experimental study on flow induced vibration of heat transfer helical tube bundle[J]. Atomic Energy Science and Technology, 2008, 42(Suppl.): 468-471(in Chinese).

[2] ANAGNOSTOPOULOS P. Numerical study of the flow past a cylinder excited transversely to the incident stream, Part 1: Lock-in zone, hydrodynamic forces and wake geometry[J]. Journal of Fluids and Structures, 2000, 14(6): 819-851.

[3] PLACZEK A, SIGRIST J F. Numerical simulation of an oscillating cylinder in a cross-flow at low Reynolds number: Forced and free oscillations[J]. Computers & Fluids, 2009, 38: 80-100.

[4] GABBAI R D, BENAROYA H. An overview of modeling and experiments of vortex-induced vibration of circular cylinders[J]. Journal of Sound and Vibration, 2005, 282: 575-616.

[5] 冯志鹏，张毅雄，臧峰刚. 直管束流固耦合振动的数值模拟[J]. 应用数学和力学，2013，34(11):1 165-1 172.

FENG Zhipeng, ZHANG Yixiong, ZANG Fenggang. Numerical simulation of fluid-structure interaction for tube bundles[J]. Applied Mathematics and Mechanics, 2013, 34(11): 1 165-1 172(in Chinese).

[6] SARPKAYA T. A critical review of the intrinsic nature of vortex-induced vibrations[J]. Journal of Fluids and Structures, 2004, 19: 389-447.

[7] 冯志鹏，臧峰刚，张毅雄. 双弹性管流固耦合振动的数值模拟[J]. 原子能科学技术，2014，48(8)：1 428-1 434.

FENG Zhipeng, ZANG Fenggang, ZHANG Yi-xiong. Numerical simulation of fluid structure interaction in two flexible tubes[J]. Atomic Energy Science and Technology, 2014, 48(8): 1 428-1 434(in Chinese).

[8] KHALAK A, WILLIAMSON C H K. Motions, forces and mode transitions in vortex-induced vibrations at low mass-damping[J]. Journal of Fluids and Structures, 1999, 13(7-8): 813-851.

[9] SKOP R A, BALASUBRAMANIAN S. A new twist on an old model for vortex-excited vibrations[J]. Journal of Fluids and Structures, 1997, 11: 395-412.

[10] LI Tian, ZHANG Jiye, ZHANG Weihua. Nonlinear characteristics of vortex-induced vibration at low Reynolds number[J]. Commun Nonlinear Sci Numer Simulat, 2011, 16: 2 753-2 771.