铅铋流体作用下螺旋管束的流致振动机理研究

王正阳 瞿叶高 余晓菲 陈 果 黄 旋 彭志科

1（上海交通大学机械与动力工程学院 上海 200240）

2（中国核动力研究设计院 成都 610213）

液态铅铋合金具有较低的熔点和较高的沸点，以及具有一定的化学惰性和良好的中子性能、抗辐照性能和传热性能等，正在逐渐替代水和空气，成为一种新型的热交换介质，同时还可以结合螺旋管道在其传热性和安全性方面的优势，成为核反应堆未来应用的重要选择之一。然而，铅铋流体引起螺旋管道产生的流致振动现象可能会导致管道破裂，发生核反应堆泄露事故，严重影响核反应堆安全运行。因此，减小避免螺旋管道在铅铋流体作用下的流致振动对铅铋反应堆的安全运行至关重要。

流体外流经过管道时会导致管道发生流致振动，管道的振动响应又会反作用于流场，形成流固耦合效应。与内流管道振动不同，外流作用在管道上时，会在管道后方产生漩涡，漩涡的脱落和涡的交替出现会导致涡激振动和锁频效应。管束的流致振动与管束的排列方式、流速和流体介质的物理性质密切相关。Paidoussis 等［1］对横向流作用下的弹性结构振动以及其稳定性问题进行了详细的论述。González Mantecón等［2］研究了反应堆核燃料板在轴向流作用下的流固耦合振动。Chen等［3］采用实验方法研究了螺旋管蒸汽发生器扇形管束在外部流体冲击作用下的振动响应特性。Connors［4］研究了直管管束的流致振动问题，并提出了流弹失稳临界流速的半经验公式。Pettigrew 等［5−6］开展了两相流作用下管束的流致振动实验研究，提出了两相流作用下管束的流弹失稳的表达式。冯志鹏等［7］对流体诱发单根和多排直管束产生的动力学响应及流场特性进行了数值模拟，发现并列管有相同响应特性而串列管下游受上游管尾流及涡脱落影响。吴建邦等［8］对单根螺旋管道进行了流固耦合数值计算，计算了不同支撑数量下管道的流弹失稳比和湍流抖振比，发现支撑数大于等于4 时，可以避免流弹失稳和湍流抖振的发生。缪洪康等［9］对换热器U型弯管进行了双向流固耦合数值仿真分析，结果表明，在初始一定时间内，弯管的位移响应和附近的空气流动趋于稳定，同时弯管会出现外扩-回弹等周期性变化特征，管道最终将在一个微小区间内做“8”字往复运动。Yuan 等［10］建立了流体（水）中螺旋盘管的流固耦合数值计算模型，采用有限体积法、大涡模拟和松流固耦合计算方法分析了螺旋盘管的流致振动响应结果，数值计算结果与实验结果符合较好。

针对铅铋流体的流动特性，国内外一些学者开展了相关的研究。Seale［11］采用改进的k-ε湍流模型对内通道铅铋流动进行了数值模拟，得到了轴向速度场、温度场分布，观察到了内通道湍流交混和二次流现象。 Cheng 等［12］利用计算流体动力学（Computational Fluid Dynamics，CFD）软件对三角形栅格和四边形栅格内液态铅铋的湍流性能进行对比分析，分析了不同湍流模型、雷诺数、燃料组件结构条件下，各通道内液态铅铋的水力特性。张勋等［13］研究了液态铅铋双自然循环回路中的瞬态响应和稳定性，并用FLUENT 程序进行了计算验证。吴学民等［14］基于双向流固耦合的有限元方法，开展了液态重金属反应堆的载液堆本体的地震响应分析，研究了冷却剂类型、支撑形式和高径比等因素对堆容器内冷却剂晃动效应的影响。Gu 等［15］开发了铅铋液（Lead Bismuth Fluid，LBE）冷却反应器多物理场耦合程序（Multi-physics coupling-Lead Bismuth Fluid，MPC-LBE），将点动力学模型和燃料销传热模型集成到自主开发的CFD程序中，并利用该程序对池式反应器进行整体模拟，同时从整体性能和系统性能方面进行了验证，结果表明：该程序可对LBE 冷却池式反应堆涉及的过程进行热工水力和安全分析。

虽然，关于管束的流致振动问题已有很多的研究，但针对螺旋管束的研究较少，而且研究的流体介质多为空气和水，对铅铋流体的研究多集中在传热方面。铅铋流体作用下的螺旋管束流致振动研究尚未见报道，其流致振动物理规律和机理仍有待揭示。本文采用数值仿真方法，以铅铋流体外掠单排螺旋管束为研究对象，研究不同铅铋流速工况下管束中各管道的振动响应和流场中的压力脉动响应机制，揭示流场中涡量的演化过程，阐明螺旋管束在铅铋流体作用下的流致振动机理，对实际工程中铅铋流体作用下的螺旋管道设计提供参考依据。

1 流固耦合数值方法

1.1 流体控制方程

铅铋流体在流动过程中可以按照不可压缩牛顿流体来近似处理，其流动遵循质量、动量和能量守恒定律。本文主要研究的是管道的流致振动问题，其热交换量可以忽略不计，故不考虑能量守恒方程，控制方程如下：

连续性方程：

动量方程：

式中：ρ为管外流体的密度，kg·m−3；ui和uj是速度分量，m·s−1；p是压力，Pa；t是时间，s；fi为单位质量体积力，N·kg−1；μ为动力粘度，Pa·s。

本文利用压力基求解器进行求解，基于改进后的压力耦合方程组的半隐式算法（Semi-Implicit Method for Pressure Linked Equations，SIMPLE），采用二阶迎风格式将控制方程在各网格点处进行离散，求解获得各个时刻流场速度与压力，以及流体作用在结构上的流体力。

1.2 管道振动控制方程

在铅铋流体作用下，螺旋管束会产生流致振动。这里不考虑热传导效应，采用有限元法可建立弹性管道的离散动力学方程：

式中：Q为管道的节点位移向量；Ms、Cs、Ks分别为管道的质量矩阵、阻尼矩阵和刚度矩阵；F为作用于管道表面的流体外载荷向量。

通过对流体控制方程的求解可以得到流体作用在管道表面上的流体载荷，将流体载荷施加到管道上，采用Newmark方法在离散的时间点对上式进行求解。

1.3 耦合控制方程

基于有限体积法和有限元法建立铅铋流体和管道的数值离散动力学模型后，在进行流固耦合计算时，有限元求解器通过其耦合交界面向流体域传递位移数据，更新流固耦合边界，流体有限体积求解器通过其耦合交界面向管道表面传递载荷数据。在流固耦合交界面上，流体和固体的热流密度、速度、温度、应力等要满足协调和平衡关系。本文不考虑热流量与温度的交换，耦合界面上的控制方程如下：

式中：n是流体和管道表面耦合交界面的法向矢量；σs和σf是铅铋流体和管道在耦合交界面上的应力。

2 计算模型及参数

2.1 几何模型

铅铋流体作用下的螺旋管束模型（图1）。考虑5根壁厚1 mm、螺距30 mm的管子组成的单排螺旋管束，管子的外径（D）为20 mm，5 根螺旋管道半径分别为0.325 m（1 号管道，Tube 1）、0.355 m（2 号管道，Tube 2）、0.385 m（3 号管道，Tube 3）、0.415 m（4号管道，Tube 4）和0.445 m（5 号管道，Tube 5），管子之间的间隙（P）为30 mm。流体区域为一个以扇形为底面的柱体，扇形的宽度（W）为0.15 m，内侧半径（R1）与外侧半径（R2）分别为0.31 m 和0.46 m。为了能够捕捉单排螺旋管束上游与下游的流场分布特征，流体区域的高度（H）为1 m，管道位于模型的中间部位，即L=0.5 m。螺旋管的旋转角度与扇形柱体的角度θ设置为50º。建立了如图1所示的笛卡尔坐标系，铅铋流体自上而下横向掠过5根螺旋管道，沿着坐标系中的−y方向，与重力方向一致（本文计算中考虑了流体的重力作用）。扇形柱体上顶面为铅铋流体入口面，下底面为出口面，四周的外表面为固定壁面，流体域与管束的交界面为流固耦合面，螺旋管道两端面均为固定壁面，即管道两端固定于壁面上。

图1 单排螺旋管束模型Fig.1 Single row spiral tube bundle model

2.2 物性参数

在实际工作中，铅铋流体和管道的物性参数随着温度发生改变。本文计算中不考虑热传递效应对系统的影响，取温度T=270 ℃时流体和固体的物性参数（表1）。

表1 物性参数Table 1 Physical property parameter

铅铋是低普朗特数Pr流体，在热交换方面与传统的水介质有很大的不同，但本文研究中忽略热传导效应。文献［16］采用k-ε湍流模型，建立了弹性管束壳程流体诱导振动的双向流固耦合模型并对其进行数值仿真计算。计算结果表明：不同流速条件下的数值求解结果与实验数据基本一致，验证了该模型的有效性。因此，本文采用标准k-ε湍流模型来模拟铅铋流体的湍流流动。

2.3 网格划分

采用ICEM软件对铅铋流体进行结构化网格划分（图2）。为了捕捉到螺旋管周围的流动细节，对每根管子附近的网格进行了局部细化，以保证边界层网格的厚度满足仿真计算的要求。

图2 管道和流体网格划分 (a) 流体网格，(b) 流体局部细化网格，(c) 管道网格Fig.2 Meshing of pipe and fluid (a) Fluid mesh, (b) Refining mesh, (c) Pipe mesh

为验证计算结果对于网格单元数目的依赖性，本文对5组不同网格单元数目的计算结果进行了对比。为了实现管道与流体的双向流固耦合计算，将螺旋管附近的流体域设置为动网格区域，采用平滑（Smoothing）和网格重划分（Remeshing）方法，弹性系数设置为0.5。选取的进口流速为1.93 m·s−1，综合考虑流体网格和固体网格单元数量对计算结果的影响，对比5 号管在y方向的加速度均方根值，结果如表2所示。

从表2 可以看出，当弹性管道有限元网格单元数目一定时，流体网格112万和195万对应的加速度均方根值分别为0.082和0.088，相对误差较小，可以认为取流体网格数目为112 万时计算结果已经收敛；当流体网格一定时，固体网格4.6 万和10.5 万的计算结果分别为0.082和0.085，相对误差较小，可以认为固体网格4.6万时计算结果已经收敛。因此，本文数值计算选择流体网格112万、固体网格4.6万的模型进行计算。

表2 网格无关性验证Table 2 Grid independence validation

3 结果讨论与分析

本文研究了不同进口流速工况时，铅铋流体作用下单排螺旋管束的振动加速度响应特性和流弹失稳现象。 计算的流速工况共有7 个，分别为0.83 m·s−1、1.4 m·s−1、1.93 m·s−1、2.4 m·s−1、2.73 m·s−1、3.2 m·s−1和3.7 m·s−1。出口处采用压力出口：P=0。计算的时间步长为0.001 s，总时间为2 s。

3.1 螺旋管束振动响应分析

以进口流速为0.83 m·s−1的工况为例，分析各根管道在铅铋流体作用下的振动响应。提取1 号（Tube 1）、3号（Tube 3）和5号（Tube 5）3根典型管道跨中位置处的监测点沿着y方向的振动加速度时域响应数据。截取时域响应数据中的稳态部分并对该部分数据进行快速傅里叶变换（Fast Fourier Transform，FFT），得到三根管道的振动加速度频域响应曲线（图3）。结果表明：管道的振动响应以单一频率为主，同时伴随有小的倍频出现；1 号和5 号管道的振动响应最大值发生在基频10 Hz附近，而3号管道的振动则是在基频的二倍频20 Hz处取最大幅值。这是因为螺旋管束在铅铋流体的作用下会出现漩涡脱落的现象，同时不同管道之间的漩涡脱落会产生干涉，且不同位置处的管道受到干涉的影响程度也不同，导致各根管道的振动主频出现差异。

图3 加速度响应曲线Fig.3 Frequency response of acceleration

为了分析螺旋管束在铅铋流体作用下产生的漩涡脱落及其演化过程，将扇形柱体切片，研究漩涡脱落在该平面内的演化情况。在该平面的3 号管道（Tube 3）附近设置压力监测点Point 1～6，提取该区域的压力脉动响应数据进行分析，其中监测点4（Point 4）与监测点6（Point 6）落在漩涡脱落干涉区域。如图4所示，选取φ=25º处沿着扇形半径方向且与xoz平面垂直的平面进行切片，记为S平面。

图4 流场的截面示意图Fig.4 Schematics of the cross section of flow field

图5给出了进口流速为1.93 m·s−1的工况下S平面上的5 个不同时刻对应的流场涡量图。图5 中的流场涡量云图由Q准则计算得到。通过对比不同时刻的涡量云图可以发现，流场的涡量分布具有明显的周期性。根据涡量的变化，可以得到漩涡脱落的频率为25 Hz，与该工况下1 号和5 号管道的振动加速度响应的主频一致。这说明在上述来流速度下，管道的振动响应主频率主要是铅铋流体绕流经过管道时出现的周期性漩涡脱落导致的。

图5 不同时刻涡量云图 (1) t=1.905 s，(2) t=1.925 s，(3) t=1.944 s，(4) t=1.964 s，(5) t=1.984 sFig.5 Vorticity contours at different times (1) t=1.905 s, (2) t=1.925 s, (3) t=1.944 s, (4) t=1.964 s, (5) t=1.984 s

下文分析低流速工况下中间管道（Tube 3）在基频的二倍频处的振动响应机理。以进口流速为1.4 m·s−1的工况为例，提取各压力监测点的数据，对其稳定部分的数据做FFT 分析，得到图6 所示的压力响应曲线图。结果表明：4号点（Point 4）与6号点（Point 6）的最大频率峰值位于基频的二倍频36 Hz处，而其他监测点的主频均为基频18 Hz。这表明在漩涡脱落干涉区域、压力脉动的振动频率以基频的二倍频为主。该区域的压力脉动对3号管道产生影响，导致3 号管道也在二倍频处出现最大的振动峰值。为了说明漩涡脱落干涉对3 号管道振动的影响，对一个漩涡脱落周期内流场中涡量变化进行分析。图7 为1.872～1.926 s 内S截面上涡量随时间的变化情况。在1.872 s 时刻，漩涡1（Vortex 1）从2 号管道（Tube 2）右侧脱落，漩涡2（Vortex 2）从3号管道（Tube 3）右侧脱落；在1.884 s时，第一次有漩涡经过漩涡脱落干涉区域；在1.892 s 时，漩涡3（Vortex 3）从3号管道（Tube 3）左侧脱落，漩涡4（Vortex 4）开始从4 号管道（Tube 4）左侧脱落，在1.91 s 时，第二次有漩涡经过漩涡脱落干涉区域。在整个流动周期内，3 号管道附近的流体漩涡脱落干涉区域同时受到2号管道与4号管道附近流体漩涡脱落的影响，会有两个脱落的漩涡经过，导致该区域的压力脉动在二倍频36 Hz 处的幅值最大，流场的压力脉动作用在3 号管道上，致使其也在对应的二倍频处出现最大的振动响应峰值。

图6 脉动压力响应曲线Fig.6 Pressure response

图7 涡量随时间的演化过程Fig.7 Evolution of vorticity contours with time

下文研究管道的振动响应频率随来流速度变化规律。提取各工况下1号、3号和5号管道跨中位置处沿y方向的加速度响应数据，再对其稳定区间内的数据进行FFT 变换，得到每种工况下3 根管道各自振动的主频率，绘制频率-速度图（图8）。结果表明：随着进口流速的升高，1 号管道和5 号管道的振动响应主频逐渐升高，且每根管道振动的主频都对应着该流速工况下的漩涡脱落频率。

图8 不同管道振动主频随进口流速的变化曲线Fig.8 Variation of dominant vibration frequency of different tubes with inlet velocity

根据漩涡脱落频率经验公式［17］可知，漩涡脱落频率与管径、流体流速和斯特罗哈数有关，而斯特罗哈数又只与雷诺数有关，因此在流场中雷诺数一致的条件下，管道的漩涡脱落频率相同。3 号管道则有所不同，在进口流速为2.73 m·s−1时，管道的振动响应主频出现了下降，在此之前，管道的振动响应主频为各个流速工况下的漩涡脱落频率的2 倍，但当流速到达2.73 m·s−1之后，管道的振动主频与对应工况下的漩涡脱落频率一致。这是因为随着来流流速的增大，管道之间的间隙流速也相应地升高，流体速度增大后，其对漩涡脱落干涉区域的冲击干扰作用增强，导致脱落的两个漩涡融合为一个大涡（图9），且大涡的运动频率与漩涡脱落频率一致。因此，当来流速度增大到一定值时，中间管道的振动频率锁定于管道周围漩涡脱落的频率。

图9 不同时刻涡量云图(v=3.2 m·s−1) (1) t=1.944 s，(2) t=1.955 s，(3) t=1.966 s，(4) t=1.977 s，(5) t=1.988 sFig.9 Vorticity contours at different times (v=3.2 m·s−1) (1) t=1.944 s, (2) t=1.955 s, (3) t=1.966 s, (4) t=1.977 s, (5) t=1.988 s

3.2 螺旋管束流弹失稳

螺旋管束在流体作用下，随着来流速度增大，管道的振动也会发生相应的变化。当流速超过某一值时，管束振动的振幅会突然增大，出现流弹失稳现象，该流速对应于管束的流弹失稳临界流速。图10为不同进口流速工况下，1号、3号与5号管道y方向上的加速度振动幅值随流速的变化关系曲线。图10 中的加速度幅值为每个工况下振动稳定区间内加速度时域响应数据的标准差。结果表明：在低流速情况下，各根管道的加速度幅值接近于0 且各不相同，当来流速度超过3.2 m·s−1时，管道的加速度幅值出现了一个激增，当铅铋来流速度为3.7 m·s−1时，三根管道的振动幅值分别为来流速度3.2 m·s−1时的管道振动响应幅值的50倍、130倍和220倍。此时，螺旋管束在铅铋流体的冲刷下出现了流弹失稳的现象。

图10 不同管道加速度幅值随进口流速的变化曲线Fig.10 Variation of acceleration amplitude of different tubes with inlet velocity

进一步分析管束出现流弹失稳时，各根管道的振动加速度响应特性。提取来流速度为3.7 m·s−1工况下1 号、3 号和5 号管道跨中位置处沿y方向的加速度时域响应，对其稳定区间内的数据进行FFT 变换。管道的振动加速度响应曲线（图11）。结果表明：三根管道的振动响应主频均在400 Hz 左右，这表明各管道的振动通过流场形成了强耦合作用。图12 给出了单排螺旋管束在铅铋流体中的湿模态计算结果。结果表明：螺旋管束的第4 阶流固耦合振动固有频率为407 Hz，与上述3 根管道的振动响应主频一致。因此，当管束发生流弹失稳时，管道的振动响应频率与湿模态频率接近。

图11 不同管道加速度频域曲线Fig.11 Frequency domain curve of acceleration

图12 单排螺旋管束在铅铋流体中的湿模态Fig.12 Wet mode of single row spiral tube bundle in cross lead bismuth fluid

4 结语

本文基于有限体积法和有限元法建立了铅铋流体作用下螺旋管束的流固耦合振动数值模型，并基于该模型研究了管束的振动响应与流场压力脉动之间的关联关系，阐明了不同流速下螺旋管束振动响应和流场涡量的演化规律。主要结论如下：

1）在一定的来流速度下，铅铋流体掠过单排螺旋管束时会导致管道发生涡激振动，两侧管道的振动频率与漩涡脱落的频率一致，且该频率随着来流速度的上升逐渐增大。

2）低流速工况下，螺旋管之间的漩涡脱落存在干涉现象，中间螺旋管振动响应受干涉影响最大，位于流场中间的管道的振动主频率为漩涡脱落频率的两倍。中间管道下方两侧存在漩涡脱落干涉区域，该区域会同时受到上游两根管道的漩涡脱落的影响，在一个漩涡脱落周期内先后出现两个漩涡，导致压力脉动频率是漩涡脱落频率的两倍，进而导致中间管道在二倍频处产生最高的响应峰值。随着流速的升高，漩涡脱落干涉对管道的振动响应影响逐渐减弱。

3）当铅铋流体的来流速度超过一定值时，螺旋管束会出现流弹失稳现象，管道的振动幅值呈现大幅增长，管束的振动响应频率与其湿模态固有频率接近。