考虑双向流固耦合的换热器弹性管束流致振动特性三维数值模拟研究

2022-09-23 00:58李岳霖艾诗钦刘月婵

振动与冲击 2022年17期

李岳霖，孙超,3，艾诗钦，刘月婵,2

(1.哈尔滨理工大学测控技术与通信工程学院，哈尔滨 150080；2.哈尔滨工程大学水声工程学院,哈尔滨 150001；3.哈尔滨工程大学动力工程及工程热物理博士后科研流动站，哈尔滨 150001)

随着节能技术的发展，换热器在化工领域的应用愈加广泛，带来了显著的经济收益。然而为了满足生产需要，不断追求高热效率，管束挠性越来越大，振动越来越敏感，易发生疲劳损坏[1-3]。因此对换热器设备结构进行设计、评估和检修时，必须对换热管流致振动的相关特性问题进行研究分析，以保证设备的可靠性。

壳程流体流过管束时，管束后部会形成交替脱落的漩涡，对圆柱产生周期性的交变应力，打破管束原有安全状态。长时期受到流体冲刷的影响，能量不断积累，超过安全临界值，易产生剧烈爆炸。轴向流绕流结构体存在流致振动这种现象于1950年末就被发现，可惜并未得到足够的重视[4]。

赖永星[5]研究指出只有在流速远高于正常流速场合，纵向流致振动才需要考虑。正常情况下，横向流速可能引起较大振动，对结构产生危害最大，人们更为关注横向流致振动产生的危害；晋文娟[6]通过水洞实验指出单柔性管进行研究存在不足之处；吴倩倩等[7]对内外充液管束流致振动进行双向流固耦合数值模拟研究，发现单向流固耦合数值模拟结果并不十分可靠；Zhao等[8]对相邻两根管束的相互影响进行分析，得出两相邻管的节径比与位移之间的关系；Katinas等[9]通过实验分析，刘景伟等[10]通过数值计算与实验对比分析，研究了线性排列与交错排列管束两种不同排列方式的振动特性；冯志鹏等[11]针对正方形顺排列弹性管束流固耦合系统三维数值模型，实现了不同弹性管束模型的流体力及振动响应特性分析。

目前，通过大量的理论与实验研究，学术界比较认同的管束流致振动机理有以下四种：漩涡脱落，流体弹性不稳定，紊流抖振与声共振[12]。现有研究基本是以这些理论为支撑进行开展，大多数都是针对外流域不同形状管束、单管及排列的不同方式对流致振动特性进行单向耦合分析，且针对不同节径比开展管束集群研究较少。因此，笔者就换热器管束流致振动问题展开研究，针对包含不同节径比管束的振动特性进行双向流固耦合分析，为换热器安全运行提供重要的设计依据。

1 换热管的固有频率

1.1 单跨管的固有频率

首先进行理论值计算，计算模型如图1，考虑到换热管支撑条件的影响，对其结构进行合理简化，即换热管视为直管，两端管板视为固定支撑，管内部流体介质为水，求解换热管振动问题可以等效为连续梁振动问题。换热管束固有频率计算式[13]为

图1 单跨管模型Fig.1 Single-span-tube model

(1)

式中：fi为i阶固有频率，Hz；l为跨距，m；EI为管的抗弯刚度；λi为由下面式(2)决定

(2)

式中：N为跨数；n=0,1,2,……,(频带为基数时)；n=1,2,3,……,(频带为偶数时)。

m为管的单位长度质量，kg/m。计算方法如(3)所示

m=mt+mi+mo

(3)

式中：mt为单位长度空管质量，kg/m；mi为换热管内单位长度流体(本文为水)的质量，kg/m；mo为单位长度管所排开管外流体(本文为空气)的质量，kg/m，相较于前两项质量可忽略不计。换热管材料为0Cr19Ni9钢，密度：7 900 kg/m3，泊松比：0.247，弹性模量：195 GPa。换热管长：l=1 200 mm，内径：d=15 mm，外径：D=16 mm。将数值代入公式，计算结果列于表1中。

表1 单跨管固有频率理论值与数值计算值Tab.1 Theoretical and numerical values of natural frequency of single span pipe

张杰等[14]对U型充液管道进行了流固耦合模态分析，结果表明随着模态阶数的增加，管道的固有频率逐渐增大，无耦合计算的管道固有频率大于流固耦合状态时的结果。本文运用数值模拟仿真方法对单跨管固有频率进行分析，计算出前八阶固有频率。为保证计算精度提高计算准确性，流体域与单跨管结构划分为六面体结构网格，结构化网格可以很容易实现区域边界拟合，最大单元：24 mm，最小单元：0.24 mm，最大单元增长率：1.3，网格划分总数为44 980个。网格划分结果如图2所示。

通过表1可知，单跨管固有频率的理论值和数值计算值随着阶次的增加误差逐渐增大，前八阶频率误差基本低于3%。结果表明，利用数值仿真对单跨管固有频率计算不会引起较大误差。

(a) 单跨管网格

1.2 四跨管的固有频率

实际生产设计中考虑到管子热膨胀和方便安装等问题，换热管与支撑板间存在间隙，管束因流体绕流后会有支撑没有接触上，被称为非有效支撑。首先假设全部支撑都是有效支撑，换热管材料不变，四跨管折流板简化为简支支撑，两端管板简化为固定支撑如图3所示。四跨管总长：L=1 200 mm，每跨长：l=300 mm，内径：d=15 mm，外径：D=16 mm。固有频率理论计算结果和数值仿真计算结果如表2所示。

图3 四跨管模型Fig.3 Four-span tube model

表2 四跨管固有频率理论值与数值计算值Tab.2 Theoretical and numerical values of natural frequency of four span pipe

通过表2可以看出两种计算方法所得固有频率误差仍基本低于3%，可以采用该方法对多跨管固有频率进行计算。相同阶数下，多跨管固有频率普遍高于单跨管固有频率，说明管束跨数会对固有频率计算有影响。

1.3 非等跨管的固有频率

换热器设备正常运行时，设备各个构件达到最佳状态，具有高度可靠性和稳定性，随着运行时间提升，设备出现不同程度损坏，可靠性和稳定性得不到保证，其换热管束中多跨管某一支撑换热管的折流板失效时，就会出现非等跨的有效支撑现象，出现较大安全生产隐患；以四跨管为例，将失效支撑折流板用数字0代替，有效支撑折流板用数字1代替，对非等跨直管固有频率计算分析，有三种不同位置处失效，即“011”型；“101”型；“100”型如图4所示。

(a) “011”型

采用数值仿真模拟对三种失效型前四阶固有频率进行计算，得出相应失效状态下振型图，图5为“011”振型，图6为“101”振型，图7为“100”振型，固有频率计算结果如表3所示。

(a) 一阶振型

表3 多跨管不同失效型固有频率计算结果Tab.3 Calculation results of natural frequencies of different failure modes of multi span pipes

据表3计算结果显示，相同失效支撑数量下失效位置不同的多跨管固有频率相接近，而随着支撑失效数的增加，同阶固有频率会降低[15]。

2 换热管的流致振动

2.1 网格独立性及数值方法验证

在进行数值模拟计算之前，通常要进行网格的独立性检验。网格的疏密及质量直接影响着计算结果的精度。通常认为当网格密度达到一定程度后，继续增加网格单元数量对于计算结果的影响非常小，此时可认为网格疏密对于计算结果的影响可以忽略。

流场与单跨管结构都采用六面体结构网格进行划分。由于在近壁面速度变化梯度较大，所以本文对近壁面处网格进行细化，为达到计算精度要求，本文设置了5层边界层，第一层网格高度应满足y+≈1，因此在Re=64 000时，第一层网格高度δ≤0.04 mm，近壁面网格划分如图8所示。

图8 圆柱xy面局部加密网格Fig.8 Local encryption grid of cylindrical xy plane

由表4可知，179 496网格单元与232 956网格单元的各项数据十分接近。因此，为了保证计算及后处理的效率和求解精度之间的平衡，采用179 496网格单元进行数值模拟。

表4 不同网格数量下的各项仿真数据对比Tab.4 Comparison of various simulation data under different grid numbers

为了验证数据的准确性，本文与文献[16]中的试验数据进行对比，结果如图9所示。计算了不同外流速度冲击下的St数。由表5可知，仿真数据与试验数据误差较小。因此，可以采用该数值模拟方法进行换热管的流致振动分析。

表5 不同外流速度下St数的对比分析Tab.5 Comparative analysis of St number under different outflow velocities

图9 圆柱绕流的St数与雷诺数Re的关系曲线Fig.9 Relationship curve between St number and Reynolds number Re of circular cylinder flow

2.2 水作用下单跨管的双向耦合

外流速度作用下换热器管束将产生变形，漩涡脱落对管束产生力的作用，从而管束位移呈现一定的变化规律。Khushnood等[17]通过实验研究指出在一定壳程流速范围内，湍流是引起管内振动的主要激励机制。

为研究三维管束的流致振动规律，本文设置外流速度为4 m/s，研究单跨管在该流速冲击下的振动特性。

管束两端保持固定支撑，换热管材料不变。换热管长：L=500 mm，内径：d=15 mm，外径：D=16 mm。为保证流场流动得到充分的发展，流场进口与圆柱中心距离为2D，流场宽度为10D，流场出口距离圆柱中心8D。网格划分如图10所示。

(a) 管结构与管外流场端面网格

根据雷诺数计算公式得出4 m/s速度下雷诺数为Re=64 000，流域特征处于亚临界区，选用k-ε湍流模型可以达到计算要求，考虑到水作用下的双向流固耦合，选取时间步长为5×10-4s，进行瞬态的数值模拟。

图11 速度xy截面云图Fig.11 xy section velocity nephogram

图12 压力xy截面云图Fig.12 xy section pressure nephogram

由图13可以看出，换热管相对于zx平面进行往复摆动，管束中部位移量最大，运动轨迹与后文图18相吻合，呈现“8”字型。

(a) t=0.441 s

据图14与图15可知，升、阻力曲线与x、y向位移曲线在时域上的变化趋势具有相似性，达到一定稳定状态后呈现简谐输出。对其稳定部分分别作FFT变换，即可得两者各自频域曲线，如图16、图17所示。

图14 升、阻力时域曲线图Fig.14 Time domain curves of lift and drag

图15 x、y向位移时域曲线图Fig.15 Time domain curve of displacement in x and y direction

图16 升、阻力频率曲线图Fig.16 Lift and drag frequency diagram

从图16可知，升力频率为49.75 Hz，阻力频率为104 Hz，符合阻力与升力频率成两倍关系理论[18]。据x、y向位移频域曲线显示，y和x向位移响应频率分别与之升、阻力频率一一对应，x向位移主要是由阻力引起，y向位移主要是由升力引起。

管束振动达到稳定状态后，将横坐标对应x向位移，纵坐标对应y向位移，绘制三维管束振动位移曲线，得到振动位移曲线图，如图18所示。可以发现位移曲线呈现“8”字型，这是因为x向位移响应频率是y向位移响应频率的两倍，且振幅较大，出现典型的Strouhal振型，位移曲线图呈现上下对称的两个部分。

图18 振动位移曲线Fig.18 Vibration displacement curve

改变外流域流速大小，得到不同速度下的振动位移曲线如图19所示，可以看出三种速度下同样出现Strouhal振型，振动位移曲线幅值与外流速成正比关系，随流速增大而增大。

图19 不同速度下振动位移曲线Fig.19 Vibration displacement curves at different speeds

2.3 水作用下四跨管的双向耦合

跨度较长的换热器设备在理想条件下折流板视为简支有效支撑，换热器管束两端保持着固定有效支撑。

将2.2中单跨管简化等分为四跨管，中间折流板视为简支支撑，两端管板视为固支支撑模型，从上到下每跨分别命名为Ⅰ跨、Ⅱ跨、Ⅲ跨、Ⅳ跨如图20所示。四跨管同样取来流速度为4 m/s的水对其进行冲击，不同跨振动位移曲线如图21所示。

图20 四跨管模型图Fig.20 Four span pipe model diagram

(a)Ⅰ跨

据图21所示，互为对称跨的振动位移曲线“8”字型相似且位移幅度相等，位处中间两跨的位移幅度要大于两侧跨。

表6与图16、图17所得数据对比可知，相同流速冲击下，多跨管振动频率与单跨管基本一致。四跨管升、阻力频率满足两倍关系，位移响应频率与之升、阻力频率存在一一对应关系。

表6 各跨升、阻力频率及轴向位移响应频率计算结果Tab.6 Calculation results of rise,resistance frequency and axial displacement response frequency of each span

3 管束的流致振动

3.1 网格独立性验证

本节针对3节径比的管束进行网格验证。选取不同数量的四种网格进行评估。

由表7可知，303 256网格单元与344 215网格单元的各项数据十分接近。因此，为了保证计算及后处理的效率和求解精度之间的平衡，数值模拟计算时采用303 256网格进行。

表7 不同网格数量下的各项仿真数据对比Tab.7 Comparison of various simulation data under different grid numbers

3.2 节径比1.5时管束的流致振动

三维单管能从一定程度上反映振动特性，但并未考虑其他管的影响，存在一定的局限性。受管与管之间的距离影响，每根管的流场不能得到充分发展，管间流体会相互影响，相互作用，导致换热器管束流场和受力与三维单管存在差异。

选取五根管组成的管束群，探讨不同节径比下管束的流致振动特性，分别对节径比为1.5和3.0的两种管束群进行建模和网格划分。每根换热管模型和材料与2.2中相同，流域为水，流场宽度为8D，流场长度为20D，每两根相邻换热管间距为1.5D。计算模型如图22所示。

图22 管束计算模型Fig.22 Tube bundle calculation model

整个流场中管束周围存在较大的速度梯度变化，为使计算结果更加精确，对管束周围进行加密细化网格处理，流域与换热管结构选择六面体结构化网格划分，换热管网格最大单元：3.29×10-3m，最小单元：5.06×10-5m，最大增长率：1.05；流域网格最大单元：9.36×10-3m，最小单元：1.01×10-3m，最大增长率：1.1，网格总数为303 256。网格划分如图23所示。