引 言
核反应堆蒸汽发生器等换热器中存在大量的管束结构，横掠管束的流体可能导致管束结构发生振动，是蒸汽发生器设计中的重要课题．管束结构的流致振动会导致管束与支撑件以及管束间的碰撞与磨损，磨损严重时会导致一回路冷却剂流出，造成巨大的经济损失．因此，预测横向流作用下的管束结构的稳定性对于反应堆的安全设计具有重要意义．
目前认为，横向流作用下弹性管束结构有四种典型的流致振动机制，分别为周期性漩涡脱落、紊流抖振、声共振以及流体弹性不稳定(流弹失稳)[1]．在上述四种机制中，对于管束结构危害最为严重的是流弹失稳[1-5]．流弹失稳现象的产生，是由于横掠管束的流速超过临界速度，导致流体内的动能不断累积，使得管束振幅迅速增加．这种大振幅的剧烈振动在蒸汽发生器等换热装置设计及服役过程中是必须避免的．
由于管束结构流弹失稳问题的重要性，近几十年来诸多学者开展了大量的理论和实验研究，建立了各种流弹失稳理论模型来预测流弹失稳的临界速度和失稳边界，用于指导蒸汽发生器的工程设计，但这些理论模型都需要依赖实验确定相关模型参数．Connors[4]、Blevins[5-6]先后对一、二维管排建立了流弹失稳分析的拟静力理论模型，在该模型中通过实验获取了不同质量阻尼系数下的临界流速，建立了管束参数与失稳流速的简化模型．Price和Paidoussis[7]、Blevins[5]在拟静力模型的基础上，进一步考虑了弹性管的位移对流体作用力的影响，提出了流弹失稳分析的拟定常模型，通过实验获取管阵的升力与阻力系数以及相应空间导数作为模型的输入参数．Chen[8-9]基于非定常流体理论，给出了较全面的流弹失稳预测非定常模型，但该模型中需要实验确定的参数太多，限制了其在工程中的应用．Lever和Weaver[10-11]采用一维非定常流理论，导出了横向流作用下管束结构流体力的解析表达式，建立了半解析模型．该模型中关于流体弹性力的部分参数需要由实验确定，其中相位滞后(phase lag)是半解析模型中需要通过实验确定的最关键的参数．因为管阵流弹失稳的理论模型中均需要实验来获取输入参数，使得流弹失稳理论分析的周期长、成本高、对实验的依赖性强．近年来，随着CFD和数值风洞技术的飞速发展，采用CFD获得各种流体力参数，从而部分或完全代替风洞或水洞实验成为可能．因此，发展工程中适用的无需实验输入的管束结构流弹失稳预测模型，成为迫切的工程需求，同时也具有重要的学术价值，是流致振动研究的前沿和热点问题．Khalifa等分别通过实验[12]与CFD模拟[13]提出了在半解析模型中采用速度相位滞后的新思路，并验证了用CFD获取速度相位滞后的可行性，给出了一种统一的但非常复杂的速度相位滞后函数表达式．由于在半解析模型的积分运算中会频繁涉及到相位滞后函数，因此这种复杂的相位滞后函数会在编程实现和工程应用中带来困难．
本文以横向流作用下平行三角形管阵与正三角形管阵为例，提出了一种改进的CFD-半解析模型混合的流弹失稳预测方法．采用CFD获得流域内速度场的变化情况与弹性管振动之间的相位差，并根据流管局部坐标与折合速度分段建立简化的速度相位延迟函数，代入半解析模型来预测流弹失稳的临界速度，理论预测的失稳速度阈值与文献[14-16]的实验结果吻合良好，验证了所发展的解析-数值混合流弹失稳预测方法的正确性和可行性．
1 横向流下管束结构流弹失稳的半解析模型
1.1 流体弹性力推导
以平行三角形管阵为例，流管模型示意图如图1所示[10-11]．曲线坐标s用以表示流动路径，流管的平均截面积在整个流管上假设为常数，等于最小间隙处的面积A0(对于平行三角形A0=Pr-D，其中Pr=P/D为管阵节距、D为管子直径)．瞬时流管面积A(s, t)会随着管子的振动而改变，A(s, t)包含两项：平均项A0和波动项a(s,t)，
A(s,t)=A0+a(s,t)．
(1)
同样地，速度和压力定义为
U(s,t)=U0+u(s,t), P(s,t)=P0+p(s,t)．
(2)
图1 流管模型
Fig. 1 The streamtube model
假设在足够远的上游位置sl，扰动可以忽略，速度和压力可以分别视为常数U0和P0．流管假定在sa位置开始贴附于管子并在ss位置分离．
上游扰动函数a(s,t)可以表示为
a(s,t)=a(sm,t)f(s)eiφ(s)，
(3)
其中a(sm,t)是最小间隙位置的面积扰动，是管子几何尺寸的函数．f(s)为人工衰减函数，表示从贴附位置到其他位置的扰动衰减．采用相函数φ(s)来考虑扰动延迟效应(即相位滞后)．相位滞后最早由Chen等[17-18]提出，指的是周期振动的管子与其周边伴随振动的流体存在的相位差．经典半解析模型采用线性相位滞后函数来表示管子振动对流管面积扰动的影响，在实验中通过流场图像后处理获得，实验和CFD模拟表明线性相位函数与实际存在较大误差[12-13]，线性相位滞后函数表示为
(4)
式中Ur=U0/(ωnl0)为折合流速，ωn是振动的固有频率．Khalifa等[12-13]研究表明可以通过CFD获取不同管阵中的速度相位滞后函数φ(s)，代替经典模型中流管面积扰动延迟函数来推导流弹力并预测流弹失稳．
取流管中的一段不可压缩流体的控制体，其连续性方程为
(5)
式中U(s,t)是流体速度向量，n(s)是垂直于控制体表面的单位向量．将相位滞后与面积扰动代入式(1)后再代入式(5)，得到流管中的速度U(s,t)．控制体的线性动量方程如下：
(6)
其中A=Ads表示控制体，∑F=-∮P(s,t)n(s)dA为作用在控制体的外力之和．将相位滞后、流管面积、流管速度代入式(6)，得到流管中的压力P(s,t)．很多研究已经证明流体弹性不稳定通常首先在升力方向产生[2-3]，因此本文只考虑管子在升力方向的运动，在管子与流体贴附的面积上积分压力项可以得到管子受到的流体弹性力Fy：
Fy=
P(s,t)cos(β(s))ds，
(7)
式中β(s)是流体附着表面的法向矢量角度．
1.2 振动方程与失稳判据
假设弹性管在升力方向即y方向以固有频率做简谐运动：
(8)
其中
为流管1与流管2作用在弹性管上的合力(对于平行三角形管束
将y方向的运动方程表示为
y(t)=Yeiωt，
(9)
式中ω为振动的复频率ω=ωR+iωI，ωR实部代表振动频率，而虚部ωI表示振幅的指数衰减或增长．即ωI<0为失稳的必要条件，ωR=0对应于静力失稳，ωR≠0对应于动力失稳．将式(9)代入式(8)并整理，得到
(10)
其中δ为阻尼对数衰减率， 其与结构固有频率
阻尼比ζ及阻尼系数c的关系为δ=2πζ=πc/(mωn)．将流体力
代入式(10)，并对实部与虚部进行分离，可以得到如下流弹失稳临界条件：
(11)
其中
由上式通过迭代求解就可得到发生动态失稳情况下的折合速度Ur和质量阻尼参数mδ/(ρD2)之间的关系．开始假设频率比ω/ωn为1，并给定Ur值，可以计算出一个流体力Fy，然后通过式(11)求出mδ/(ρD2)和ω/ωn，如果ω/ωn与初始假设值的差大于预定的残差，则更新ω/ωn后重复上述计算过程．这样便可计算出各个给定的Ur值对应的质量阻尼参数，从而对相应管阵的流弹失稳阈值进行预测．
2 CFD-半解析流弹失稳分析方法
2.1 CFD数值模拟
以正三角形排列管束为例详细说明CFD计算过程，平行三角形管束的模拟完全一致．正三角形管束节距比为Pr=1.375，流体介质为空气，密度为1.197 kg/m3，动力黏度系数为1.807 5×10-5 kg·m·s．平行三角形与正三角形管阵计算模型图如图2(a)、2(b)所示，其中中心处标记管为弹性管，其余管为刚性管．由HYPERMESH生成正三角形管阵网格模型如图3(a)所示，采用非结构化网格，边界层网格为四边形网格，第一层网格高度为0.09 mm，确保y+为1，共11层，增长率1.2．经过网格无关性验证最终网格数量约为45万．图3(a)中区域放大效果如图3(b)所示，其中浅色区域网格为弹性管周围流体区域网格，为保证该区域数值模拟时的准确性，对网格进行了加密划分．该区域为后文瞬态分析中的动网格区域，所以被单独划分出来．CFD求解器采用FLUENT，湍流模型为RANS模型中的RNG k-ε模型．
(a) 平行三角形管阵,Pr=1.54 (b) 正三角形管阵Pr=1.375
(a) The parallel triangular array, Pr=1.54 (b) The equilateral triangular array, Pr=1.375
图2 平行三角形与正三角形管阵的CFD计算模型示意图
Fig. 2 CFD simulation model of parallel triangular array and equilateral triangular array
图3 正三角形管阵的CFD网格模型
Fig. 3 The CFD meshing model for the equilateral triangular array
为获取速度相位滞后函数，通过FLUENT的UDF给弹性管在升力方向(y方向)施加简谐振动位移，Y(t)=0.01Dsin(ωt),其中D为弹性管直径，平行三角形和正三角形排列弹性管固有频率ωn分别为25 Hz和9.42 Hz．对整个流场进行瞬态分析，获取流域内不同位置处速度随时间的变化情况．
图4展示了正三角形管阵，当Ur=6时一个周期内的弹性管周围区域的速度云图．可以明显看出管的位置发生了改变．弹性管正上方与正下方的高流速区域也在随着弹性管的振动，周期性地发生变化．当弹性管向上运动时(如图4(b)所示)，挤压上方流管，导致上方流管截面变小，上方流速变大；当弹性管向下运动时(如图4(d)所示)，挤压下方流管，导致下方流管截面变小，下方流速变大．
(a) t=0 (b) t=T/4
(c) t=T/2 (d) t=3T/4
图4 正三角形管阵中弹性管振动的一个周期内的速度场变化情况
Fig. 4 The velocity field variation over one cycle of tube vibration in the equilateral triangular array
注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同．
2.2 速度相位延迟函数获取及流弹失稳分析
为了探寻流域内速度场的变化情况与弹性管振动之间的相位差，需要确定获取流体流速的位置．图5为正三角形管阵内流体流速检测点在流管坐标系下的位置．其中s*=s/D，中心处的管子为弹性管．通过管阵几何模型的数学推导，可以获得这些点的位置．
从图5中的任一位置s*处获取流域的速度时程曲线，并与弹性管的位移时程曲线进行对比，可以得到该位置处速度时程曲线与弹性管位置变化时程曲线的相位差，即为该位置处的相位滞后φ(s*)．相位延迟不仅与流域内的位置有关，还与流域内的流速有关．因此，本文对不同流速下的管束流域情况进行了数值模拟．当Ur=1，2，3，4，5，6时，平行三角形管阵中不同位置处的相位滞后如图6(a)所示；当Ur=1，2，3，4，5，6时，正三角形管阵中不同位置处的相位滞后如图6(b)所示．
图5 正三角形管阵中获取相位滞后的位置
Fig. 5 The position to obtain phase lag in the equilateral triangular array
(a) 平行三角形管阵 (b) 正三角形管阵
(a) The parallel triangular array (b) The equilateral triangular array
图6 不同流速下不同位置处的相位滞后
Fig. 6 The phase lag between different positions and different velocities of the parallel triangular array and the equilateral triangular array
根据图6的相位滞后曲线、不同流管局部坐标s*与折合速度Ur，本文给出了两种管束排列方式相位滞后函数的分段函数．对平行三角形管束：
(12)
正三角形管束中，折合速度对于相位滞后函数的影响较小，可得到关于流管局部坐标s*的分段函数为
(13)
式(12)和(13)是非常简单的分段函数，非常易于编程实现和半解析模型的积分运算，将其取代式(4)中经典的线性相位滞后函数代入流弹失稳半解析模型，求解两种管阵在节距比Pr均为1.375情况下的流弹失稳临界速度曲线，如图7所示．其中平行三角形管束实验结果取自文献[14-15]的实验结果，正三角形管束的实验值取自文献[16]的实验结果，两种管束结构流弹失稳理论预测结果和实验结果吻合良好．
图7 节距比Pr为1.375平行三角形与正三角形管束CFD-半解析模型混合流弹失稳临界速度预测
Fig. 7 Fluidelastic in stability velocity prediction based on the phase lag model for a parallel triangular array (Pr=1.375) and a equilateral triangular array (Pr=1.375)
3 结 论
本文提出了一种改进的CFD仿真与半解析方法混合的管束结构流弹失稳预测方法．采用CFD方法计算了不同流速和不同排列方式下管束结构全域流速和压力分布．提取弹性管周围不同位置的流域由于弹性管振动引起的流速变化及速度相位延迟，并根据速度和流管位置坐标将其近似为简单的分段函数，代替经典的半解析流弹失稳模型中的线性相位延迟函数，成功预测了横向流作用下间距比Pr为1.375的平行三角形与正三角形管束结构的流弹失稳阈值．理论结果与文献已有实验结果吻合良好，验证了所发展的CFD-半解析模型混合流弹失稳预测方法的正确性和可行性．用类似的研究思路和方法，可以将CFD数值模拟和流弹失稳分析的拟定常模型、非定常模型结合，发展出一套无需实验输入的管束结构流弹失稳分析方法，可极大地节约工程设计和分析所需的时间和成本．
[1] 姜乃斌, 冯志鹏, 臧峰刚. 核工程中的流致振动理论与应用[M]. 上海: 上海交通大学出版社, 2018. (JIANG Naibing, FENG Zhipeng, ZANG Fenggang. Theory and Application of Flow-Induced Vibration in Nuclear Engineering[M]. Shanghai: Shanghai Jiaotong University Press, 2018. (in Chinese))
[2] WEAVER D S, ZIADA S, AU-YANG M K, et al. Flow-induced vibrations in power and process plant components: progress and prospects[J]. Journal of Pressure Vessel Technology, 2000, 122(3): 339-348.
[3] PAIDOUSSIS M, et al. Real-life experiences with flow-induced vibration[J]. Journal of Fluids and Structures, 2006, 22(6/7): 741-755.
[4] JR CONNORS H J. Fluid elastic vibration of tube arrays excited by cross flow[C]//ASME Symposium on Flow-Induced Vibration in Heat Exchanger, Winter Annual Meeting. 1970.
[5] BLEVINS R D. Flow-Induced Vibration[M]. New York: Van Nostrand Reinhold Co, 1977.
[6] BLEVINS R D. Fluid elastic whirling of a tube row[J]. Journal of Pressure Vessel Technology, Transactions of the ASME, 1974, 96(4): 263-267.
[7] PRICE S J, PAIDOUSSIS M P. An improved mathematical model for the stability of cylinder rows subject to cross-flow[J]. Journal of Sound and Vibration, 1984, 97(4): 615-640.
[8] CHEN S S. Instability mechanisms and stability criteria of a group of circular cylinders subjected to cross-flow, part Ⅰ: theory[J]. Journal of Vibration and Acoustics, 1983, 105(1): 51-58.
[9] CHEN S S. Instability mechanisms and stability criteria of a group of circular cylinders subjected to cross-flow, part Ⅱ: numerical results and discussions[J]. Journal of Vibration and Acoustics, 1983, 105(2): 253-260.
[10] LEVER J, WEAVER D. On the stability of heat exchanger tube bundles, part Ⅰ: modified theoretical model[J]. Journal of Sound and Vibration, 1986, 107(3): 375-392.
[11] LEVER J, WEAVER D. On the stability of heat exchanger tube bundles, part Ⅱ: numerical results and comparison with experiments[J]. Journal of Sound and Vibration, 1986, 107(3): 393-410.
[12] KHALIFA A, WEAVER D, ZIADA S. An experimental study of flow-induced vibration and the associated flow perturbations in a parallel triangular tube array[J]. Journal of Pressure Vessel Technology, 2013, 135(3): 030904.
[13] KHALIFA A, WEAVER D, ZIADA S. Modeling of the phase lag causing fluidelastic instability in a parallel triangular tube array[J]. Journal of Fluids and Structures, 2013, 43: 371-384.
[14] WEAVER D S, EL-KASHLAN M. The effect of damping and mass ratio on the stability of a tube bank[J]. Journal of Sound and Vibration, 1981, 76(2): 283-294.
[15] WEAVER D S, GROVER L K. Cross-flow induced vibrations in a tube bank-turbulent buffeting and fluid elastic instability[J]. Journal of Sound and Vibration, 1978, 59(2): 277-294.
[16] AUSTERMANN R, POPP K. Stability behaviour of a single flexible cylinder in rigid tube arrays of different geometry subjected to cross-flow[J]. Journal of Fluids and Structures, 1995, 9(3): 303-322.
[17] CHEN S S. Flow-induced in-plane instabilities of curved pipes[J]. Nuclear Engineering and Design, 1972, 23(1): 29-38.
[18] CHEN S S, JENDRZEJCZYK J A, LIN W H. Fluid-elastic instability of rectangular tube arrays subjected to liquid cross flow[C]//Fourth International Conference on Pressure Vessel Technology. Vol 2. 1980.