0 引　言

管壳式换热器作为舰船的关重设备，广泛应用于动力系统、冷却系统、空调系统等多个重要环节，对舰船的正常运行和性能发挥起着不可或缺的作用。换热管束处于换热器复杂流场环境，通常受到湍流、周期性旋涡脱落以及流体弹性不稳定性等多种流体因素诱导引起的异常振动，其中流体弹性不稳定性，是造成换热管束振动并快速磨损破坏的关键因素之一^[1]。

流体弹性激振最早由Roberts^[2]于1960年发现，后由Connors^[3]和Lever等^[4]分别提出了流体刚度引起的位移模型和流体阻尼引起的惯性模型，目前仍作为换热管束振动失稳工程计算的重要设计准则。近年来，流体弹性不稳定性是引发的换热器管束振动失效频发，同时受2012年美国圣奥诺弗雷核电站（SONGS）发生的蒸汽发生器管束故障^[5]事件影响，换热管束的流体弹性激振研究受到广泛关注。

针对换热管束的流体弹性激振的研究方法主要分为理论分析、数值仿真以及试验研究，其中试验研究偏重于对Connors和Lever方程进行理论修正^{[6 - 7]}；数值仿真主要采用非定常、大涡模拟等CFD与FEM联合仿真对换热管束进行双向流固耦合计算，研究管束的流致振动特性^[1]；理论分析主要针对射流转换模型、准静态模型、非稳态模型^[8]、流管模型和准稳态模型等模型进行理论推导及数值求解^[7]。当然，目前也有很多学者将三者进行结合，开展基于理论、仿真以及试验的综合研究，对换热管束的流体弹性激振机理特性具有更进一步的认识^{[8 - 10]}。

国内外学者近几年对换热管束的流体弹性激振稳定性开展深入研究，其中李先达等^[11]利用高速相机视觉测量技术对不均匀横流作用下的传热管束流致振动响应进行了试验研究，得出了不均匀横流作用下传热管束流弹失稳临界流速，结果表明：Connors方程计算结果偏保守40%；杨世豪等^[12]基于非定常流体力建立了传热管流致振动的理论模型，分析得出了非定常流体力对传热管流弹失稳机制的影响，研究结果表明支承方式会影响失稳临界流速，而且非定常流体力作为一种强迫力在流弹失稳之前，会引起管束的“拍振”现象；Chen等^[10]通过将Hassan时域求解模型与CFD参数获取方法集成，开发了一种换热管束流体弹性不稳定性预测模型，Khalvatti等^[13]测试结果对比，验证了模型的有效性和准确性；Petr^[14]基于准稳态模型，利用流动对流和粘性效应考虑流体弹性激励的延时效应，评估了流体弹性频率对阻尼控制的流体弹性不稳定性预测的具有一定影响；刘建等^[15]利用CFD计算管束流体动力学参数建立换热管束流固耦合模型，该方法可以一定程度反映管束和流体相互作用的主要特征，并对流弹失稳临界流速有一定的预测能力；Zhao等^[1]提出了一种兼顾具有衰减函数和显式流体弹性力半解析时域流体弹性激振模型，支持频域稳定性分析和真实时域响应分析，适用于刚性阵列中单个弹性管或阵列中多个柔性管的配置，并与试验数据对比，验证了其准确性；Tan等^{[16 - 17]}开展了换热管束流致振动试验，并利用Ansys CFD建立了采用大涡模拟和动网格的换热管束数值模型，讨论了流体弹性不稳定的振动特性，仿真与实验结果吻合较好。

以上各类研究中，试验研究需要花费大量人力财力，而且存在试验规律分析存在偶然因素、部分物理量不可测试等问题，另外，使用有限元分析软件（CFD-FEM）来解决多管束流体弹性系统的相关动力学问题并不合适，例如在计算管束几何形状变化的装置时，需要大量的时间和配备强大处理器的计算集群。

因此，本文首先采用粘性涡域方法求解单管束的流体激励载荷；然后基于对换热管束流体载荷近程特性和管束对称性的分析，组建多管束的流体动力相互作用矩阵；其次，结合管束振动微分方程，建立换热管束横向绕流弹性激振数学模型，并利用稳定性条件求解多管束稳定性临界流速；最后，利用相似性准则，研究了由典型管束单元求解大型对称管束稳定性临界流速计算方法。

1 基于配置法的多管束流体激励载荷计算方法

流体绕流静止圆柱体的复势（$ {U_\infty } $为二维空间中无穷远处均匀流体的流速）为：

$ {F(z) = {U_\infty }\left( {z + \dfrac{{{R^2}}}{z}} \right) + \dfrac{1}{{2 {\text{π}} i}}\left( {\displaystyle\sum\limits_{n = 1}^N \Delta {\varGamma _n}\ln \left( {z - {z_{en}}} \right) - \displaystyle\sum\limits_{n = 1}^N \Delta {\varGamma _n}\ln \left( {z - {z_{in}}} \right)} \right)}。$

(1)

式中：$ z\left( {x,y} \right) $为二维空间中的坐标；$ {z_{en}} $、$ {z_{in}} $分别为外部和内部涡的坐标；$ \Delta {\varGamma _n} $为流场中第n个外部涡在逆时针旋转时的强度。公式第一项是稳定流绕圆柱体的复势，第二项考虑了流体中自由涡的存在，第三项则考虑了圆柱体内虚构涡的存在。

从无穷远处以速度$ {U_\infty } $流向封闭任意轮廓的均匀流体，作用在该轮廓上的复力可以由外部和内部涡速度表示：

$ {F_X} + i{F_Y} = - i\rho \sum\limits_{k = 1}^N \Delta {\Gamma _k}\widehat {{V_{ek}}} + i\rho \sum\limits_{m = 1}^N {\frac{{\mathrm{d}\Delta {\Gamma _{im}}}}{{\mathrm{d}t}}} {Z_{im}}。$

(2)

式中：$ \Delta {\varGamma _k} $为流场中第k个外部涡的强度；$ \Delta {\varGamma _{im}} $为位于物体轮廓上第m个涡的强度；$ {{\text{Z}}_{im}} $为位于物体轮廓上第m个涡的坐标；$ {\hat V_{ek}} $为笛卡尔坐标系下的复速度。

考虑到流体绕半径为R的圆柱体流动的复势表达式为：

$ \begin{aligned} F(z) = &{U_\infty }\left( {z + \frac{{{R^2}}}{z}} \right) - i{V_{yy}}\left( {z - \frac{{{R^2}}}{z}} \right) + \\ &\frac{1}{{2{\text π} i}}\left( {\sum\limits_{n = 1}^N \Delta {\varGamma _n}\ln \left( {z - {z_{en}}} \right) - \sum\limits_{n = 1}^N \Delta {\varGamma _n}\ln \left( {z - {z_{in}}} \right)} \right)。\end{aligned} $

(3)

式中：$ {V_{yy}} $为圆柱体旋转速度。

可以得到作用在半径为R的封闭运动圆形轮廓上复力的表达式：

$ \begin{aligned} {F_X} + i{F_Y} = &\rho \sum\limits_{n = 1}^N \Delta {\varGamma _n}\left( {{V_{en}} - {V_{in}}} \right) - \\ &i\rho \left( {\sum\limits_{n = 1}^N \Delta {\varGamma _n}\left( {{U_{en}} - {U_{in}}} \right) - 2 {\text{π}} i{R^2}\frac{{\mathrm{d}{V_H}}}{{\mathrm{d}t}}} \right)。\\ \end{aligned} $

(4)

式中：$ {\hat V_{en}} = {U_{en}} - i{V_{en}} $；$ {\hat V_{in}} = {U_{in}} - i{V_{in}} $。

实际的热交换器管束是由相同的直管道或空间弯曲的多跨弹性管道组成的系统，这些管道与内部和外部的流体相互作用，典型管束截面横向绕流计算模型如图1所示。

其中单根管束k的正面阻力$ {F_X} = {C_X}\rho R{U^2} $，升力$ {F_Y} = {C_Y}\rho R{U^2} $，$ \rho $为热载体密度，R为管道半径，L为管道跨度，C_X为沿X轴方向的力系数，C_Y为沿Y轴方向的力系数，U为远离管束的流体流速，其中F_X和F_Y可以由式(4)快速求得。

2 换热管束横向绕流弹性激振数学模型

考虑流体附加质量和阻尼的情况下，换热器管束在流体中的弯曲振动可以用Bernoulli-Euler方程组描述：

$ m\frac{{{\partial ^2}\bar w(z,t)}}{{\partial {t^2}}} + b\frac{{\partial \bar w(z,t)}}{{\partial t}} + EJ\frac{{{\partial ^4}\bar w(z,t)}}{{\partial {z^4}}} = {\bar F_z}(z,t)。$

(5)

式中：m为考虑管内液体的管道单位长度质量；b为考虑管道材料内部摩擦和结构阻尼的参数；EJ为管道横截面的抗弯刚度；$ (\overline{w}(z,t)=[w_{1x},w_{1y},...,w_{Ny}]^{\mathrm{T}} $为管束中管道轴线相对于未受干扰位置的位移向量；$ {w_{jx}}({w_{jy}}) $为第j根管道轴线在$ Ox(Oy) $方向上的位移；$ {\bar F_z}(z,t) $为横向流绕流管道时作用在管道上的分布流体动力载荷向量，z为沿管道长度方向的坐标，t为时间。

在不具体设定弹性管道端部固定条件，并且仅考虑前几阶弯曲振动中的某一阶振动的情况下，式(5)的解可以表示为：

$ \bar w(z,t) = a\bar X(tu/2R)\psi (z/l)。$

(6)

式中：$ \bar X = col({x_1},{y_1},...,{x_n},{y_n}) $为一个2N维向量，表示管道相对于未受干扰位置的位移；$ \psi (z/l) $是与所考虑的弹性管道（作为两端固定梁）振动形式相对应的固有函数。

假设管束按照其某一阶固有振型进行振动，仅考虑管道前几阶弯曲振动中的某一阶，忽略随机强迫振动，可得到描述管束固有振动的微分方程组：

$ \ddot{x}+2\xi\omega_0\dot{x}+\omega_0^2x=\mu_1\boldsymbol{C}(\tau)。$

(7)

式中：$ {\omega _0} = {\omega _{p1}}R/U $为单根管道的无量纲固有频率；$ {\omega _{p1}} $为所考虑振动形式下的弯曲振动固有频率；$ \xi $为相对阻尼；$ {\mu _1} = \rho \cdot {R^2}/m $为无量纲质量参数；$ \boldsymbol{C}(\tau) $为流体动力相互作用矩阵；$ \tau = tU/(2R) $为无量纲时间；$ x $、$ \dot x $、$ \ddot x $分别为管道截面的位移、速度和加速度向量。

在初始条件下，对式(3)进行拉普拉斯变换，并利用不稳定力呈线性的假设，方程组具有以下形式：

$ ({p^2} + 2\xi {\omega _0}p + \omega _0^2){x^*} = {\mu _1}S(p){x^*}。$

(8)

式中：$ {x^*} $为管道位移向量的拉普拉斯变换像；p为变换参数；$ S\left( p \right) $为流体动力相互作用矩阵$ \boldsymbol C(\tau ) $的变换像。

管束不稳定振动方程(8)的解，与流体动力相互作用矩阵$ \boldsymbol{S}\left(p\right) $的特征向量一致。在将上述矩阵化为对角形式，并消去矩阵$ \boldsymbol{S}\left(p\right) $的特征向量后，同时考虑到$ {\lambda _j}(p) $是矩阵$ \boldsymbol{S}\left(p\right) $的第j个特征值，方程(8)的特征方程可以表示为：

$ \prod\limits_{j = 1}^{2N} ( {p^2} + 2 \cdot \xi \cdot {\omega _0} \cdot p + \omega _0^2 - {\mu _1} \cdot {\lambda _j}(p)) = 0。$

(9)

在假设$ \xi $和$ {\mu _1} $为小参数的情况下研究方程(9)。根据Lyapunov稳定性理论，对于系统初始平衡状态而言，特征方程(9)的所有根都具有负实部是必要且充分条件。在稳定区域的边界上，$ p = i\omega $，此时只需研究方程(9)中的一个因子即可，该因子包含具有最大虚部的特征值$ \lambda (p) $。仅考虑相对于小参数$ \xi $和$ {\mu _1} $不高于一阶的小量时，由方程(9)可得：

$ \left\{ {\begin{aligned} &{{\omega ^2} = \omega _0^2 - {\mu _1} \cdot Re[\lambda (i{\omega _0})]}，\\ &{2\xi \omega _0^2/{\mu _1} = Im[\lambda (i{\omega _0})]} 。\end{aligned}} \right. $

(10)

方程(8)确定了管道按照最不稳定共振的振动频率，以及在给定无量纲频率下参数$ {{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}} $的临界值。如果知道$ \lambda (i{\omega _0}) $的函数关系，那么方程(10)中的第二个等式就可以用来求出所研究管束对应于$ {{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}} $值的临界速度。在相对于小参数$ \xi $和$ {\mu _1} $的一阶近似下，稳定性的必要且充分条件为：

$ \left\{ {\begin{aligned} & {\frac{{2\xi }}{{{\mu _1}}} \geqslant {{\left(\frac{{2\xi }}{{{\mu _1}}}\right)}_{cr}}}，\\ & {{{\left(\frac{{2\xi }}{{{\mu _1}}}\right)}_{cr}} = \frac{{Im[\lambda (i{\omega _0})]}}{{\omega _0^2}}} 。\end{aligned}} \right. $

(11)

方程(11)中的第2个式子确定了参数$ {\left( {{{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}}} \right)_{cr}} $的临界值，即它表征了在稳定区域边界上，阻尼力与流体动力相互作用的不稳定力之间的关系。

3 多管束稳定性边界数值计算 3.1 多管束稳定性边界计算流程

为求解多管束的稳定性边界，首先要计算某一绕流速度下单个管束的流体激励载荷，并组建多管束的流体动力相互作用矩阵$ \boldsymbol C(\tau ) $；然后通过拉普拉斯变换将矩阵$ \boldsymbol C(\tau ) $转换为拉普拉斯变换像$ \boldsymbol{S}(p) $；其次求解矩阵$ \boldsymbol{S}(p) $的特征值$ {\lambda _j}(p) $，并选取虚部最大的特征值，从而确定最不稳定共振的振动频率，以及在给定无量纲频率下参数$ {{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}} $的临界值$ {\left( {{{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}}} \right)_{cr}} $；最后，通过系列化的绕流速度确定不稳定边界。多管束横向绕流稳定性边界条件计算流程如图2所示。具体计算过程如下：

1）求解流体动力相互作用矩阵$ \boldsymbol C(\tau ) $。以单个管束为对象，并保持某一管束振动，计算流体弹性激励，分别在X方向和Y方向各重复2N次（N为管道数量），能得到以流体动力相互作用矩阵$ \boldsymbol C(\tau ) $形式呈现的流体激励数据；

2）求解线性流体动力相互作用矩阵$ \boldsymbol{S}(p) $。将流体动力相互作用矩阵$ \boldsymbol C(\tau ) $经拉普普拉斯变换得到矩阵$ \boldsymbol{S}(p) $。

3）计算流体动力相互作用矩阵$ \boldsymbol{S}(p) $的特征值。求解矩阵$ \boldsymbol{S}(p) $的特征值$ {\lambda _j}(p) $，并从所有特征值中找出虚部最大的特征值。在稳定区域的边界上，$ p = i\omega $是虚数，考虑到相关参数中不高于一阶的小量$ \xi $和$ {\mu _1} $，进而得出临界参数值$ {\left( {{{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}}} \right)_{cr}} $。

4）不同绕流速度下临界参数值$ {\left( {{{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}}} \right)_{cr}} $求解，重复开展步骤1～步骤3计算求解，以此确定管束的稳定区域边界。

3.2 三管束稳定性计算

为确定流体动力相互作用矩阵，根据第2节给出的公式和第3.1节给出的计算流程，对每一根管束进行非定常流体动力的逐次计算。在计算过程中，每次仅让其中一根管道在X轴或Y轴方向上，以频率$ {\omega _j} $、小振幅按简谐规律振动，如图3所示。

图4为流体动力相互作用矩阵的一个示例关系。例如，$ \boldsymbol{C}_{2X1X} $表示第一根管道在X轴方向振动时，作用在第二根管道上沿X轴方向的流体动力；$ \boldsymbol{C}_{2Y1X} $表示第一根管道在X轴方向振动时，作用在第二根管道上沿Y轴方向的流体动力。这里$ \tau = {{tU} \mathord{\left/ {\vphantom {{tU} {2R}}} \right. } {2R}} $为无量纲时间，t为时间。

稳定性条件(11)中包含2个无量纲参数，无量纲频率$ {\omega _0} = {{{\omega _T}R} \mathord{\left/ {\vphantom {{{\omega _T}R} U}} \right. } U} $和无量纲阻尼参数$ {{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}} $，$ {{2\xi } \mathord{\left/ {\vphantom {{2\xi } {{\mu _1}}}} \right. } {{\mu _1}}} $是表征征阻尼力与流体动力之间关系的参数。为了呈现研究结果，并与其他管束稳定性研究文献类似，分别定义无量纲速度$ {V_r} $和无量纲阻尼参数$ \Delta $：

$ {V_r} = \frac{U}{{2{f_T}R}}\frac{q}{{(q - 1)}} = \frac{ {\text{π}} }{{{\omega _0}}}\frac{q}{{(q - 1)}}。$

(12)

式中：$ {V_r} $为管束道间最小间隙处的无量纲平均流速；U为远离管束的横向来流速度；$ {f_T} = {{{\omega _T}} \mathord{\left/ {\vphantom {{{\omega _T}} {2 {\text{π}} }}} \right. } {2 {\text{π}} }} $为管束的固有振动频率；无量纲参数$ q = {s \mathord{\left/ {\vphantom {s {2R}}} \right. } {2R}} $，表征管束排的密集程度。

$ \Delta = {{ {\text{π}} \xi } \mathord{\left/ {\vphantom {{ {\text{π}} \xi } {2{\mu _1}}}} \right. } {2{\mu _1}}}。$

(13)

式中：$ \xi $为管道振动的衰减系数。

图5为15种不同的管道无量纲绕流速度下三管排（管束密集程度为$ q = s/d = 1.41 $）的计算结果。通过求解方程(9)并分析其根来确定三管系统的稳定性（用☆星号表示稳定）或不稳定性（表示不稳定）（见图5）。将不同的ξ值和6个六维线性流体动力相互作用矩阵的特征值代入，并依据Lyapunov稳定性准则，得到了图5中的管束稳定性曲线（实线）。方程(10)得到的点相当精确地落在图4的稳定曲线（实线）上。图5中还展示了根据文献[18]的实验数据，通过直接分析密集程度为$ q = s/d = 1.41 $的单排管束无量纲速度特性得到的临界速度测量结果，可以看出，计算结果与文献[18]中的实验数据具有良好的一致性，相对误差小于15%。因此，所呈现的数值数据证实了关于线性流体动力相似性的主要假设，以及基于这些假设研究大型管束稳定性的方法的正确性。

3.3 五管束稳定性计算

管束中主要是相邻管道之间存在显著的相互作用，也就是说矩阵中的大部分元素值较小，可近似看作零，即为流体动力相互作用的近程作用特性。在实际使用的大型管束中，可以选取一个由3~9根管道组成的系统（管束单元），这些管束是与第l根管道距离最近的相邻管道。实际管束的横截面存在平移对称性，这使得流体动力相互作用也具有对称性，并且这种对称性可扩展到管束的所有典型单元。在简化线性流体动力相互作用矩阵的计算过程中，会用到这一重要特性。线性流体动力相互作用矩阵可根据所阐述的方案通过物理实验或数值实验来确定。对于大型管束，只需研究一个典型单元的流体动力相互作用，然后通过计算扩展到整个管束，进而找出特定的比例因子（对与对称性偏差以及对典型单元所采用的绕流条件的偏差进行多次求和）。

还需注意的是，决定流体动力相互作用的流体物理特性是相邻管道间最小间隙处的平均流速V。这一假设与线性流体动力的近程作用特性相符。而且还需要假设，管束截面具有规则横截面布局，水动力连接的近程作用特性和平移对称性。

因此，在N = 3的系统中，流体动力的值与$ N \gg 1 $（无限多根管道的情况）系统中的流体动力值存在一个比例因子$ {({{{V^{(N)}}} \mathord{\left/ {\vphantom {{{V^{(N)}}} {{V^{(3)}}}}} \right. } {{V^{(3)}}}})^2} $，其中$ {V^{(3)}} $是三根管道组成的一排中最小间隙处的平均流速，$ {V^{(N)}} $是N根管道组成的一排中最小间隙处的平均流速。

对矩阵$ \boldsymbol{S}(p) $随N变化特性的数值研究显示，当$ N \gg 1 $时，在整个频率范围内，具有最大虚部的特征值实际上并不依赖于N。这一特性体现为，大型管束的临界速度主要取决于横截面的类型和布局（管束的密集程度），与管道数量的关联较小。

根据上述假设，可以构建由三管排的流体动力相互作用矩阵组建大型管排线性流体动力相互作用矩阵的算法。线性流体动力相似性的准则是无量纲频率$ {\omega _v} $，其定义为：

$ {\omega _v} = \frac{{{\omega _T}R}}{V} = \frac{{{\omega _0}}}{V}。$

(14)

在通过使用N = 3系统的线性流体动力来确定大型管束的线性气液动力矩阵$ \boldsymbol{S}^{(N)}(i\omega_0^{(N)}) $时，必须确保

$ \omega _v^{(3)} = \omega _v^{(N)} = {\omega _v}。$

(15)

根据所采用的假设，对于$ N \gg 1 $的系统，2N阶矩阵$ \boldsymbol{S}^{(N)}(i\omega_0^{(N)}) $的元素与$ \boldsymbol{S}^{(3)}(i\omega_0^{(3)}) $的元素关系如下：

$ {\boldsymbol S^{(N)}}\left( {i\omega _0^{(N)}} \right) = {\left( {\frac{{{V^{(N)}}}}{{{V^{(3)}}}}} \right)^2}{\tilde {\boldsymbol S}^{(N)}}\left( {i\omega _0^{(3)}} \right),\frac{{\omega _0^{(N)}}}{{{V^{(N)}}}} = \frac{{\omega _0^{(3)}}}{{{V^{(3)}}}}。$

(16)

式中：$ {\tilde {\boldsymbol S}^{(N)}}(i\omega _0^{(3)}) $为根据三管排矩阵组建的N管矩阵。

根据稳定性条件(11)，大型管束的稳定性条件可写为：

$ \left( {\frac{{2\xi }}{{{\mu _1}}}} \right)_*^{(N)} = \frac{{{\text{Im}}\left[ {{\lambda ^{(N)}}(i\omega _0^{(N)})} \right]}}{{{{(\omega _o^{(N)})}^2}}}。$

(17)

式中：$ {\tilde \lambda ^{(N)}} $是根据三管矩阵恢复的N管矩阵的特征值。在这种情况下，N管排的稳定性条件为：

$ \left( {\frac{{2\xi }}{{{\mu _1}}}} \right)_*^{(N)} = \frac{{{\text{Im}}\left[ {{{\tilde \lambda }^{(N)}}(i\omega _0^{(3)})} \right]}}{{{{(\omega _o^{(3)})}^2}}}。$

(18)

因此，通过使用在不同$ \omega _0^{(3)} $值下计算得到的$ \boldsymbol{S}^{(3)}(i\omega_0^{(3)}) $矩阵，可以得到大型管束在相应无量纲频率下的参数值。

因此，五管排稳定性计算可以基于三管排的计算结果进行组合计算。计算流程如下：

步骤1　获取三管排的线性气液动力连接矩阵；

步骤2　找出每个矩阵所有特征值的最大虚部；

步骤3　根据相应的固有频率$ {\omega _0} $，由三管排矩阵组件五管排矩阵（见图5）；

步骤4　确定无量纲系数，方程式(18)；

步骤5　计算模拟参数$ \Delta $；

步骤6　确定绕流速度V_r：$ {V_r} = {{{V^{(5)}} {\text{π}} } \mathord{\left/ {\vphantom {{{V^{(5)}} {\text{π}} } {{\omega _0}}}} \right. } {{\omega _0}}} $ 其中$ {V^{(5)}} $是五管排中管道间最小间隙处的无量纲平均流速。

基于第2节中描述的方法，对由5根管道组成的管束片段的绕流过程进行了数值模拟，图6展示了使用粘性涡域方法对由5根管道组成的管束绕流进行数值试验的结果，在数值实验中，管束的密集程度为1.41，仅对中心管道设定简谐振动，并评估作用在每根管道上的流体动力载荷的时间历程。

图7针对五管管束的流体动力相互作用矩阵的一个元素片段。展示了在第一根管道沿X轴方向振荡时，作用在第五根管道上的流体动力载荷系数，C_5x1x表示第一根管道在X轴方向振动时，作用在第五根管道上沿X轴方向的流体动力载荷系数；C_5y1x表示第一根管道在X轴方向振动时，作用在第五根管道上沿Y轴方向的流体动力载荷系数。

并按照前述由三管束求解五管束稳定性问题的求解方法，得到了五管束线性流体动力相互作用矩阵及其特征值，利用稳定性条件(11)确定了五管束的临界绕流速度，如图8所示，与试验结果^[18]对比，相对误差小于15%。

4 结　语

1）建立了一种换热管束的横向绕流弹性激振数学模型，通过本模型可以通过利用离散粘性涡域方法和稳定性准则快速计算管束的稳定性临界流速，计算结果与试验结果相对误差约15%～20%。

2）基于小参数假设和管束近程作用特性，确定了管束稳定性条件。

3）结合离散粘性涡域方法和配置法实现了管束流体载荷计算方法，并基于相似性准则和稳定性条件，建立了由典型管束单元组建大型对称布置管束稳定性临界流速的过程。

但本文中稳定性条件、线性流体动力相互作用矩阵涉及小参数假设、管束对称、管束近程作用特性等前提条件，因此，后续可针对一般型式的管束布局方案稳定性进行计算研究。