湍流抖振存在于换热器管束的整个运行过程中，以微幅振动的形式出现，导致结构某些部位与支撑板的碰摩和疲劳损伤等累积失效。管束的微幅振动对强化换热有积极的影响，可以激发管程工质的脉动以提高换热效率^[1]。管束流致振动的早期研究是以漩涡脱落和流体弹性激振为主，湍流激励为辅。针对湍流抖振的工作多是研究其对周期性激励的影响^[2-4]。计算和实验手段的限制导致经常将管束绕流场视为2D，垂直于管轴线方向的平面内进行研究^[5-6]，包括漩涡脱落和流弹激励。然而由于湍流激励复杂的三维特性，在沿管轴线方向上的分布随机，因此换热管表面的每个空间位置都存在一个时间域上的随机湍流力，数学上可以描述为多个随机过程。为了研究湍流所引起的随机振动，文献[7-8]对激励力的频谱标准化方法做过不同的探究尝试，并提出了受单位长度上已知功率谱的随机载荷作用下换热管以一阶振型和固有频率进行湍流抖振的挠度振幅预报经验公式。其中载荷的功率谱需要应用Axisa^[9]提出的等功率谱密度法来设定，该方法假设分布湍流随机载荷是服从正态分布的各态历经随机过程，进而对其进行频谱的标准化，王定标等^[10]基于该方法对大型管壳式换热器进行了随机响应分析。随着计算机技术的发展，通过流场仿真技术得到换热管上的流体力分布，进而利用随机振动方法对其响应做出预报。林家浩等^[11]提出并发展的“虚拟激励法”对于多点激励平稳随机振动方程给出了高效精确的求解方法。对于复杂结构的非平稳随机振动，其计算效率也得到大幅提高。其应用领域目前主要是工程抗震应用^[12-14]、风工程领域^[15-19]、车辆工程领域^[20]以及其他发展的领域^[21-28]。本文以管束壳程流场的CFD仿真结果中的表面压力作为载荷，首先验证载荷的随机性，并利用基于梁模型的传统随机振动计算方法和虚拟激励法分别分析了换热管不同位置的位移响应功率谱。

1 换热管随机振动理论 1.1 换热管的受迫振动

换热管结构振动分析基于将换热管抽象为细长比较大的带有中间支撑的两端固支多跨Bernoulli-Euler梁模型计算，分布载荷作用下的受迫振动微分方程为：

$ \mathrm{EI} \frac{\partial^{4} y(x, t)}{\partial x^{4}}+m \frac{\partial^{2} y(x, t)}{\partial t^{2}}+c \mathrm{EI} \frac{\partial y(x, t)}{\partial t}=f(x, t) $

(1)

式中：EI为梁的抗弯刚度；m为梁单位长度的质量；c为介质阻尼系数；f(x, t)表示时空变化的分布载荷。方程(1)在零初始条件下基于模态叠加法的时域解为：

$ y(x, t)=\sum\limits_{r=1}^{\infty} \varphi_{r}(x) \int_{0}^{l} \int_{0}^{t} h_{r}(\tau) f(x, t-\tau) \cdot \varphi_{r}(x) \mathrm{d} \tau \mathrm{d} x $

(2)

式中：φ_r(x)是横向振动的第r阶模态振型；ω_r为其对应的固有圆频率；h_r(t)是第r阶脉冲响应函数。

1.2 平稳随机响应的计算

线性结构有M个自由度，其脉冲响应函数和频率响应函数均为已知，二者为Fourier变换对。根据随机振动理论，一个平稳随机过程的自功率谱S(f)与其自相关函数R(t)为Fourier变换对。线性系统单输入F(t)、单输出y(t)的情况，响应输出的自谱为：

$ {\mathit{\boldsymbol{S}}_{\mathit{\boldsymbol{yy}}}}(f) = {\mathit{\boldsymbol{H}}^*}(f)\mathit{\boldsymbol{H}}(f){\mathit{\boldsymbol{S}}_{\mathit{\boldsymbol{FF}}}}(f) = |\mathit{\boldsymbol{H}}(f){|^2}{\mathit{\boldsymbol{S}}_{\mathit{\boldsymbol{FF}}}}(f) $

(3)

式中H(f)表示输入点与输出点之间的导纳函数，可以由频响函数求解得到。略去表达式中的自变量频率f，输入与输出之间的互谱可以表示为：

$ {{S}_{\mathit{\boldsymbol{Fy}}}} = {H}{{S}_{\mathit{\boldsymbol{FF}}}} $

(4)

$ {S_{\mathit{\boldsymbol{yF}}}} = {H^*}{S_{\mathit{\boldsymbol{FF}}}} $

(5)

可以推导多输入F(t)=[F₁(t), F₂(t), …, F_M(t)]、多输出y(t)=[y₁(t), y₂(t), …, y_M(t)]的响应功率谱矩阵为：

$ \boldsymbol{S}_{\boldsymbol{y y}}=\boldsymbol{H}^{*} \boldsymbol{S}_{\boldsymbol{F F}} \boldsymbol{H}^{\mathrm{T}} $

(6)

$ \boldsymbol{S}_{\boldsymbol{F y}}=\boldsymbol{S}_{\boldsymbol{F F}} \boldsymbol{H}^{\mathrm{T}} $

(7)

$ {\mathit{\boldsymbol{S}}_{\mathit{\boldsymbol{yF}}}} = {\mathit{\boldsymbol{H}}^*}{\mathit{\boldsymbol{S}}_{\mathit{\boldsymbol{FF}}}} $

(8)

所有涉及响应的功率谱矩阵，包括S_yy、S_Fy以及S_yF都可由激励的功率谱矩阵S_FF与导纳矩阵H的通过矩阵乘法得到。S_FF为一个M阶方阵，对角线元素为每一点输入的自谱，其他元素为2点输入之间的互谱，因为2个随机过程的互谱存在共轭性质，因此S_FF^*=S_FF^T，S_FF为Hermite矩阵。结构受已知功率谱的多点平稳随机激励的作用，借鉴结构动力学的模态叠加法，可以得到响应输出[y(t)]的相关函数矩阵为：

$ \begin{gathered} \boldsymbol{R_{y y}}(\tau)=\mathrm{E}\left[\left[y_{i}(t)\right]\left[y_{j}(t+\tau)\right]^{\mathrm{T}}\right]= \\ \sum\limits_{r=1}^{N} \sum\limits_{s=1}^{N} \boldsymbol{\varphi}_{r} \boldsymbol{\varphi}_{s}^{\mathrm{T}} \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} \boldsymbol{R_{F F}}\left(\tau+\tau_{1}-\tau_{2}\right) \cdot \\ h_{r}\left(\tau_{1}\right) h_{s}\left(\tau_{2}\right) \mathrm{d} \tau_{1} \mathrm{~d} \tau_{2} \end{gathered} $

(9)

对上式应用维纳-辛钦关系进行Fourier变换可得系统响应的功率谱矩阵为：

$ \boldsymbol{S_{y y}}=\sum\limits_{r=1}^{N} \sum\limits_{s=1}^{N} \boldsymbol{\varphi}_{r} \boldsymbol{\varphi}_{s}^{\mathrm{T}} \boldsymbol{H}_{r}^{*} \boldsymbol{H}_{s} \boldsymbol{S}_{\boldsymbol{F} \boldsymbol{F}} $

(10)

式(10)称为CQC方法，计入了所有的参振模态与耦合项。对于M自由度系统，使用N阶模态叠加法对计算进行截断，CQC的计算量为N²M²次实数乘法，在工程中推荐使用一个简化的近似方法SRSS，忽略式(10)中所有的模态耦合项：

$ \boldsymbol{S_{y y}}=\sum\limits_{r=1}^{N} \boldsymbol{\varphi}_{r} \boldsymbol{\varphi}_{r}^{\mathrm{T}}\left|\boldsymbol{H}_{r}\right|{ }^{2} \boldsymbol{S_{F F}} $

(11)

式(11)的计算量为NM²次实数乘法。SRSS仅对于参振频率全部为稀疏分布，且各阶阻尼比很小的均质材料结构才可用。大多数三维结构的参振频率很少为稀疏分布，因此SRSS的近似误差也引起广泛关注和研究。

1.3 虚拟激励法

当结构上的单点随机激励由单位简谐激励e^iωt代替时，其他点的简谐响应为He^iωt，H为两点的跨点导纳函数。令虚拟激励为$\tilde{F}=\sqrt{S_{F F}} \mathrm{e}^{\mathrm{i} \omega t}$，由线性关系可知对应的虚拟响应为$\tilde{y}=\sqrt{S_{{F F}}} H \mathrm{e}^{\mathrm{i} \omega t}$，则$\tilde{y}^{*} \tilde{y}=$$|\tilde{y}|^{2}=|H|^{2} S_{\boldsymbol{F F}}=S_{y y}$，同理$\tilde{F}^{*} \tilde{y}=S_{F y}, \tilde{y}^{*} \tilde{F}=S_{y F}$。已知多点平稳随机激励功率谱矩阵S_FF为一个M阶的Hermite矩阵，则可以按其特征系统分解为：

$ \boldsymbol{S}_{\boldsymbol{F F}}=\sum\limits_{j=1}^{M} \lambda_{j} \boldsymbol{\psi}_{j} \boldsymbol{\psi}_{j}^{* \mathrm{~T}} $

(12)

其中λ_j和ψ_j分别为S_FF的特征值和特征向量，用每一阶特征对构造虚拟激励$\tilde{\boldsymbol{F}}_{j}=\boldsymbol{\psi}_{j}^{*} \sqrt{\lambda_{j}} \mathrm{e}^{i \omega t}$就可以得到$\boldsymbol{S}_{\boldsymbol{F F}}=\sum \tilde{\boldsymbol{F}}_{j}^{*} \tilde{\boldsymbol{F}}_{j}^{\mathrm{T}}$，则虚拟响应为$\tilde{y}_{j}=\boldsymbol{H} \tilde{\boldsymbol{F}}_{j}$，可以推导出：

$ \sum\limits_{j=1}^{M} \tilde{\boldsymbol{y}}_{j}^{*} \tilde{\boldsymbol{y}}_{j}^{\mathrm{T}}=\boldsymbol{H}^{*}\left(\sum\limits_{j=1}^{M} \tilde{\boldsymbol{F}}_{j}^{*} \tilde{\boldsymbol{F}}_{j}^{\mathrm{T}}\right) \boldsymbol{H}^{\mathrm{T}}=\boldsymbol{H}^{*} \boldsymbol{S}_{\boldsymbol{F F}} \boldsymbol{H}^{\mathrm{T}}=\boldsymbol{S}_{\boldsymbol{y y}} $

(13)

式(13)与式(6)的形式相同，因此与CQC具有相同的精度，且计算量仅为M²次实数乘法。此方法被称为虚拟激励法(pseudo excitation method，PEM)。

2 管束绕流场仿真 2.1 结构模型、网格划分与仿真设置

建立19根转置三角排列管束绕流的三维模型，管子外径16 mm，有效长度0.62 m，分为2跨，节径比为1.5。基于大涡模拟(large eddy simulation，LES)湍流模型对管束绕流的计算^[24-26]已经成为目前管束绕流仿真的应用趋势，本文采用LES方法对如图 1的管束绕流模型进行仿真。针对此计算模型共进行了2套网格划分方案：方案1，流场区域采用六面体的结构化网格离散，管束表面边界层厚度0.04 mm，并以每层4%的增长率递增至3 mm，管束区和进、排气道网格细节如图 2所示；方案2，为了提高进、排气道区域网格质量，在方案1的基础上，管束区仍采用全结构化划分，进、出口过渡区域采用非结构化划分。

方案1经过计算，管束表面最大壁面y⁺值为3.83，为保证LES的壁面y⁺≤1的要求，边界层厚度降为0.01 mm，网格厚度增长率与方案1相同，方案2的网格质量更高，数量略有增加。表 1为2种网格方案的对比。

流场介质为常温空气，采用有限体积法求解控制方程，并使用3-D Smagorinsky-Lilly亚格子尺度LES模型，PISO压力速度耦合格式，PRESTO!压力离散格式，二阶迎风动量离散格式，有界二阶隐式离散化的方案应用于此模型。在入口速度为3 m/s的亚临界雷诺数条件下，Re为1.6×10⁴，出口为0 Pa，其余边界包括换热管束外壁、壳体内壁和挡板壁均为无滑移壁面条件，瞬态计算时间步长为10^-4s，总计算时长为2s。

2.2 流场仿真结果讨论

图 3所示为2个网格方案的管束壁面y⁺值分布，方案1由于管束表面流体边界层网格厚度较大，导致壁面y⁺较高，最大值为3.83；改进后的方案2边界层网格厚度为方案1的1/4，因此壁面y⁺可以达到0.73以下，空间分布也更加明显。证明方案2的网格细化使得计算精度有了很大的提高，为后续能够精确计算管束表面流体激励力提供了先决条件。2个网格方案计算得到进出口压降分别为168.5 Pa和171.7 Pa。

考察某一跨管束上压力的分布，图 4展示的是沿流速方向视角，迎流第1排管(1~5号管)、中间排管(11~15号管)与尾管(19号管)上的压力云图和壁面y⁺云图，按图 1(d)中对管子的编号排列。如图 4(a)所示，迎流的第1排管束上所受的压力分布固定，在靠近右侧位置可看出压力没有左侧高，是因为进气道与管束区相连的位置没有覆盖到整个管长，空气绕流在下游管子上会沿轴线方向扩散分布。图 4(b)可以看出在第1排管束上y⁺的分布沿管轴线方向上相当均匀，说明第1排管束的绕流流速分布基本相同，而且图 4(a)和(b)中对称的2对管，2号与3号、4号与5号管上的压力和流速分布对称性很强；中心排管束如图 4(c)、(d)，对称性完全消失，管束表面压力分布混乱，绕流进入充分发展阶段。流场下游管束周期性载荷不再主导，随机激励成为流致振动的主要因素。

为了验证计算模型的准确性，本文对仅有1号管的单管绕流场也进行了模拟，图 5(a)和(b)分别为单管绕流场速度和19管绕流场的速度，在单管尾流中有足够的空间可以形成完整的卡曼旋涡并脱落在单管上产生周期性的脉动升阻力，密排管束则不然，尾流中难以产生固定频率的漩涡脱落。单管仿真得到的升阻力系数时域与频域图如图 6所示，升力的脉动频率为74 Hz，阻力的频率为148 Hz，满足文献[6]中所述的阻力频率是升力频率的2倍关系，此时St数为0.394 67。

在与单管仿真相同的位置和边界条件下，密排管束中的迎流1号管的升阻力系数结果如图 7(a)~(c)所示。密排管束中已不能形成完整的旋涡，流体激励力的能量频带变宽，旋涡脱落受限。图 7(d)~(f)和(g)~(i)表明处于位置对称7、8号管，在阻力系数相近的情况下，升力系数反向，频谱近似。图 7(j)~(l)11号中间管与图 7(m)~(o)19号尾管的频谱中也出现若干主导频率，然而高频中的分量变得越来越多。

考察第19号尾管，每个横截面上升力或阻力分量合成一个合力就是该横截面上的2个垂直的横向作用力。尾管上的阻力和升力在时域内的分布如图 8所示，阻力随时间变化不大，在计算范围内的相关系数不低于0.99，相当于时不变分布载荷；而每个截面的升力在时间上则是随机的。

3个典型位置第1跨中点、全管中点以及第2跨中点升力数据的概率分布如图 9所示，虚线为标准正态分布概率线，基本服从均匀分布的概率曲线。每个截面升力都是一个随机过程，将升力的计算结果截取无关的10次时序样本，每个样本内包含100步连续的瞬态数据，经过计算上述3个位置随机过程的相关函数系数如图 10所示。在时间间隔相等的情况下，相关系数处于同一水平且与计算起始时间点无关，即数据的相关函数保持平稳。各个随机过程的时间平均并不相等，且不等于统计平均，因此本文中研究的管子上M个空间载荷是广义平稳随机过程，属于多点输入平稳随机问题，不同于文献[9]假设的单位长度载荷为各态历经平稳随机过程。

3 换热管随机振动响应计算结果分析 3.1 换热管的固有特性

由1.2和1.3节可知，在计算多跨换热管结构的随机振动响应时，离散分布激励的功率谱矩阵S_FF和结构自身跨点频响矩阵H是必要条件。计算符合国标^[24]的两等跨304不锈钢换热管，其有效管长0.62 m，材料密度为7 930 kg/m³，材料弹性模量为200 GPa，管外径为16 mm，管壁厚8 mm，管内流体密度为998 kg/m³，阻尼比为0.02。利用有限元方法计算其固有频率和模态，计算过程中充液圆管的抗弯刚度即为金属圆管本身的刚度，充液圆管的线质量为单位长度的金属管质量、管内流体质量与排开的管外流体质量之和，管外流体为空气，质量相对于前2项可以忽略。前10阶固有频率结果为580、841、1 484、1 878、3 919、4 545、5 936、6 701、10 226和11 223 Hz。

3.2 随机振动响应算例 3.2.1 载荷完全相干算例

载荷在空间上完全相干，则rank(S_FF)=1。令功率谱为S_FF(x₁, x₂, f)=I，载荷为空间上完全相干的白噪声，对此载荷分别应用CQC、SRSS和PEM方法截断于10阶模态进行随机响应求解，3种方法分别用时441.1、42.4和4.7 s，基本满足计算速度与模态叠加阶数的倍率关系。

由上文中的仿真结果可知，换热管受到湍流随机力作用时的能量主要集中于500 Hz以下的低频段，高频模态无法受激起振，本节中应用空间完全相干的白噪声作为载荷，理论上的响应即为系统的传递函数，且高阶模态的响应幅值远远低于低阶模态，因此计算过程中截断于高达11 223 Hz的十阶模态仅是为了验证计算效率。

图 11(a)、(b)分别为一跨中点位置的位移响应自谱和1/3与2/3位置的位移响应互谱，计算范围为0~5 000 Hz。曲线中各有3个峰值，分别对应于充液管第2、4、6阶固有频率，图 11(b)中CQC和PEM的结果中还存在3 730 Hz附近的反共振峰，SRSS的计算结果中则不存在。

图 11中峰值频率处的幅值相对误差在表 2中所列，相对误差的定义为SRSS和PEM方法与CQC方法在固有频率处幅值的差与CQC结果的比值。对于圆管上的单点响应自谱，PEM可以表现出与CQC十分接近的精度。对于两点响应间的互谱，在低于1 000 Hz的频段内，PEM的精度依旧很好；在高于1 000 Hz的频段内，PEM与CQC的峰值出现分离，并且只出现在反共振峰处，但是误差仍然远远小于SRSS。

3.2.2 载荷不相干算例

对2.2节中仿真得到的尾管升力进行功率谱计算，其rank(S_FF)=M，功率谱矩阵满秩，即载荷在空间上不相干。应用3种方法对仿真载荷作用于充液圆管模型进行响应计算，得到一跨中点位置的自谱如图 12(a)所示，500 Hz以内为载荷控制频段见图 12(b)，用相对范数误差^[29]来表征一个频段内SRSS或PEM的计算精度：

$ {\boldsymbol{\epsilon}}_{n}=\frac{\left\|\boldsymbol{V}-\boldsymbol{V}_{0}\right\|_{F}}{\left\|\boldsymbol{V}_{0}\right\|_{F}} $

(14)

式中：V为SRSS或PEM方法计算得到的某点响应自谱; V₀为CQC方法计算的该点响应自谱。载荷控制频段内，SRSS的相对范数误差为43.5%，PEM的相对范数误差为0.39%，可见PEM的精度远远高于SRSS。

图 12(a)的高频段则回归系统的传递函数，CQC与PEM的响应自谱在2 170 Hz和4 100 Hz附近存在反共振峰，SRSS的响应自谱只有在4 545 Hz处存在与CQC相近的峰值，且精度仍然低于PEM方法。

4 结论

1) 简化换热器中管束绕流的CFD模型瞬态数值模拟获取管束表面湍流激励的瞬态分布，通过本文计算在时间上是多个服从均匀分布的广义平稳随机过程，而不是文献[10]中假设理想的服从高斯分布的各态历经随机过程。作为平稳随机过程，绕流场中密排管束表面，尤其是尾部管表面的流体力主要为能量集中于低频段的宽频作用力。

2) 利用空间完全相干的白噪声作用于充液换热管上以验证3种方法计算横向振动的位移响应谱的效率，PEM的计算速度约等于SRSS的N倍、CQC的N²倍，其中N为CQC和SRSS方法中参与计算的模态截断阶数，且在系统传递函数结果中，PEM在峰值频率处幅值精度与精确算法相当。

3) 空间不相干分布载荷作用于换热管时，造成受迫振动在低频段内载荷控制的响应功率谱或高频段内传函控制的响应功率谱，以CQC结果为参考，应用PEM相对于SRSS方法均高精且高效，因此在已知换热管横向载荷谱的情况下，虚拟激励法在计算其随机响应方面具有重要的工程实际意义，可以为换热器管束无支撑跨距的设计标准提供参考依据。

本文中尚未解决折流板换热器内多壳程流动时管束受力存在的空间相干性问题及其造成的随机振动响应预报方法。