1. 引言

充液管路系统的振动噪声预报是一类典型的流/固/声耦合问题。复杂充液管路系统的声振响应常采用基于“十四方程”理论发展而来的方法求解[1]-[11]。其中，Wiggert等[1] [2]最早提出三维管道振动十四方程模型，实现了管道的轴向、弯曲和扭转运动的同时求解。Tijsseling等[3] [4]对充液管路十四方程的特征线法进行了阐述，首次提出将其与有限元法结合分析系统的响应。李艳华等[5] [6]推导出多管段的频域传递矩阵法，分析了弯管弯曲半径及角度对流固耦合振动传递特性的影响。李帅军等[7]考虑流体引起的离心力，用传递矩阵法分析了输流管道压力波传递特性。而传递矩阵法在求解长距离输流管道时会存在数值溢出和响应计算不稳定。尹志勇等[8]和吴江海等[9] [11]提出了基于ISM的船舶管路系统声振耦合频域响应预报方法，研究了充液管路系统弯管、分支管的振动带隙特性，对结构波的传播进行了研究，包括参数实验获取和理论计算，并在船体结构上进行相关实验验证。基于“十四方程”的ISM在计算复杂空间管路系统时，设计自由度少，计算效率和精度高，在空间充液管路振动响应预报方面已经广泛应用[8]-[11]。

然而在船舶等实际应用场景中，管路常与船体连接，管路系统的振动会通过连接结构传递到船体，并引起噪声辐射。因此，管路与连接的船体结构需要进行耦合求解。吴江海等[12]-[14]采用管路阻抗综合法、船体有限元法与管路支撑测试数据相结合的方法，建立了管路–支撑–船体耦合模型，研究了管路–船体耦合的振动特性，但是仅考虑了管路、平板、圆柱壳结构的支撑连接，且管路支撑边界条件需要实验获取。考虑管路与船体连接结构耦合时，难以单独应用ISM进行求解，采用FEM求解则不仅需要占用大量计算资源，且不容易揭示声振沿管路系统的耦合传递机制。刘碧龙等[15]采用解析计算方法对具有弹性障板的有限长充液直管的管口声辐射进行求解，其中管道部分采用行波法，弹性障板采用模态法，结果表明管口辐射流体声和弹性障板辐射的结构声对耦合系统的声辐射均有较大贡献。虽然解析方法物理过程清晰，便于揭示充液管路的声振耦合传播机制，但实际的船体结构复杂，解析方法不能处理复杂的管路和壳体耦合结构。总体来说，考虑管路与壳体耦合求解的文献较少。

本文给出一种ISM和FEM相结合的计算方法，即管路部分采用ISM、弹性板采用FEM的混合计算方法。以两种管路–弹性板耦合系统为例分析，研究管路与壳体耦合的声振传递机制。通过管路与弹性板耦合处振动响应结果与有限元直接计算结果对比，验证了本文方法的正确性，应用ISM-FEM混合计算方法分析了弯管和管内流体对管路弹性板耦合系统的固有频率和振动响应的影响。本文提出的ISM-FEM混合计算方法为研究与壳体耦合的复杂充液管路的声振耦合传递机制提供了一种高效的计算和分析工具。

2. 理论推导

本文以直管–弹性板和弯管–弹性板耦合系统为例，研究管路与连接结构的声振耦合的求解方法和传递机制。为简化推导，只考虑了管内流体和管壁的耦合，管口外部的流体区域假设为轻流体，其对系统的耦合作用予以忽略。其中充液管路部分采用ISM，弹性板采用FEM求解。阻抗综合法是将部件以节点形式连接，首先用有限元方法得到管路与弹性板耦合处的阻抗边界条件，利用管路与弹性板耦合处力和位移边界条件的连续性，与管路阻抗矩阵组合后得到管路-弹性板耦合系统整体的阻抗矩阵。

2.1. 振动方程及阻抗综合法

充液管路系统按其振动形式可以分为轴向振动、横向振动和扭转振动三类，3个方向的振动互不耦合，可以分别独立求解。以充液直管–弹性板耦合系统模型为例，充液直管内力、声压与位移示意图如图1所示。在直角坐标系(x, y, z)下，u_x、u_y、u_z分别为xoz平面横向位移、yoz平面横向位移、轴向位移，f_x、f_y、f_z分别为横向力和轴向力，m_x、m_y、m_z为各平面力矩，ϕ_x、ϕ_y、ϕ_z为各平面转角，p和u_f为声压和管内流体位移，本文为简化推导，研究yoz平面响应共有8个未知量，管道轴向为Z方向。

Figure 1. Schematic diagram of straight tube-elastic plate coupling system

图1. 直管–弹性板耦合系统示意图

根据文献[2]，充液直管轴向振动方程为：

$\{\begin{array}{c}\frac{\partial p}{\partial z}+{{\displaystyle \rho }}_f\frac{{{\displaystyle \partial }}^2{{\displaystyle u}}_f}{\partial {{\displaystyle t}}^2}=0\\ \frac{\partial p}{\partial t}+{{\displaystyle K}}^{*}\frac{{{\displaystyle \partial }}^2{{\displaystyle u}}_f}{\partial t\partial z}-2\mu {{\displaystyle K}}^{*}\frac{{{\displaystyle \partial }}^2{{\displaystyle u}}_z}{\partial t\partial z}=0\\ \frac{\partial {{\displaystyle f}}_z}{\partial z}-{{\displaystyle \rho }}_p{{\displaystyle A}}_p\frac{{{\displaystyle \partial }}^2{{\displaystyle u}}_z}{\partial {{\displaystyle t}}^2}=0\\ {{\displaystyle f}}_z-\mu {{\displaystyle A}}_p\left(\frac{R}{h}+\frac{1}{2}\right)p-E{{\displaystyle A}}_p\frac{\partial {{\displaystyle u}}_z}{\partial z}=0\end{array}$ (1)

充液直管yoz平面内横向振动方程为：

$\{\begin{array}{c}\frac{\partial {{\displaystyle F}}_y}{\partial z}-\left({{\displaystyle \rho }}_p{{\displaystyle A}}_p+{{\displaystyle \rho }}_f{{\displaystyle A}}_f\right)\frac{{{\displaystyle \partial }}^2{{\displaystyle u}}_y}{\partial {{\displaystyle t}}^2}=0\\ {{\displaystyle f}}_y-\kappa G{{\displaystyle A}}_p\left(\frac{\partial {{\displaystyle u}}_y}{\partial z}-{{\displaystyle \varphi }}_x\right)=0\\ \frac{\partial {{\displaystyle m}}_x}{\partial z}+{{\displaystyle f}}_y-\left({{\displaystyle \rho }}_p{{\displaystyle I}}_p+{{\displaystyle \rho }}_f{{\displaystyle I}}_f\right)\frac{{{\displaystyle \partial }}^2{{\displaystyle \varphi }}_x}{\partial {{\displaystyle t}}^2}=0\\ {{\displaystyle m}}_x-{{\displaystyle EI}}_p\frac{\partial {{\displaystyle \varphi }}_x}{\partial z}=0\end{array}$ (2)

式(1)~式(2)中，R、h、A_p和A_f分别为管路内径、壁厚、横截面积和流体横截面积；ρ_f和ρ_p分别为管内流体密度和管壁密度；K*为管内流体体积模量；I_f和I_p分别为流体和管壁的转动惯量；κ为Timoshenko梁剪切系数；G为管路剪切模量。

对轴向振动方程和yoz平面横向振动方程均采用分离变量法，具体参数见文献[2]，以轴向振动式(1)为例，管路两端(z = 0, z = L)的轴向振动传递矩阵可以写为

${{\displaystyle T}}^{FP}=B\left(L\right)B^{-1}\left(0\right)$ (3)

根据ISM传递矩阵与阻抗矩阵之间的关系[11]，左端为力状态向量，右端为位移状态向量，得到管路轴向阻抗矩阵：

$Z^{FP}=\frac{1}{{{\displaystyle T_3}}_{}}\left[\begin{array}{cc}{T}_1& -\left(T_2{T}_3-T_1{T}_4\right)\\ 1& T_4\end{array}\right]=\left[\begin{array}{cc}{Z}_1& Z_2\\ Z_3& Z_4\end{array}\right]$ (4)

式中，T_i (i = 1, 2, 3, 4)为2 × 2矩阵。

同理可以根据管路横向振动方程式(2)，可以得到横向阻抗矩阵Z^yoz。对于直管，轴向和横向振动互不耦合，将轴向和yoz面横向阻抗矩阵组合后的直管整体阻抗矩阵为：

$\left[\begin{array}{c}P\left(0\right)\\ F_z\left(0\right)\\ M_x\left(0\right)\\ F_y\left(0\right)\\ P\left(L\right)\\ F_z\left(L\right)\\ M_x\left(L\right)\\ F_y\left(L\right)\end{array}\right]=\left[\begin{array}{cccccccc}{{\displaystyle Z_{11}}}_{}& Z_{12}& 0& 0& Z_{15}& Z_{16}& 0& 0\\ Z_{21}& Z_{22}& 0& 0& Z_{25}& Z_{26}& 0& 0\\ 0& 0& Z_{33}& Z_{34}& 0& 0& Z_{37}& Z_{38}\\ 0& 0& Z_{43}& Z_{44}& 0& 0& Z_{47}& Z_{48}\\ Z_{51}& Z_{52}& 0& 0& Z_{55}& Z_{56}& 0& 0\\ Z_{61}& Z_{62}& 0& 0& Z_{65}& Z_{66}& 0& 0\\ 0& 0& Z_{73}& Z_{74}& 0& 0& Z_{77}& Z_{78}\\ 0& 0& Z_{83}& Z_{84}& 0& 0& Z_{87}& Z_{88}\end{array}\right]\left[\begin{array}{c}{U}_z\left(0\right)\\ U_f\left(0\right)\\ U_y\left(0\right)\\ {\Phi }_x\left(0\right)\\ U_z\left(L\right)\\ U_f\left(L\right)\\ U_y\left(L\right)\\ {\Phi }_x\left(L\right)\end{array}\right]$ (5)

将式(5)简写为：

$\left[\begin{array}{c}F(0)\\ F(L)\end{array}\right]=\left[\begin{array}{cc}{Z}_1& Z_2\\ Z_3& Z_4\end{array}\right]\left[\begin{array}{c}W(0)\\ W(L)\end{array}\right]$ (6)

式中，F和W为直管两端状态向量，Z_i (i = 1, 2, 3, 4)为4 × 4矩阵。

弯管采用离散方法，根据弯头处的连续性，将90˚直角弯管离散3段直管，弯管离散模型如图2所示。

Figure 2. Schematic of the discrete model of the bend

图2. 弯管离散模型示意图

传递矩阵和状态向量均在局部坐标系下，将充液管路传递矩阵T从局部坐标系(x, y, z)转换到全局坐标系 $(\bar{x},\bar{y},\bar{z})$ 下得到 $\tilde{T}$ 。管路系统的坐标转换矩阵为：

$\tilde{R}=\left[\begin{array}{cccccccc}\cos \alpha & 0& \sin \alpha & 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ -\sin \alpha & 0& \cos \alpha & 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& \cos \alpha & 0& \sin \alpha \\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& -\sin \alpha & 0& \cos \alpha \end{array}\right]$ (7)

式中， $\alpha$ 为管路直角坐标系和全局坐标系的夹角。

故充液管路在全局坐标系下传递矩阵为：

$\tilde{T}={\tilde{R}}^{-1}T\tilde{R}=\left[\begin{array}{cc}{{\displaystyle T}}_1& {{\displaystyle T}}_2\\ {{\displaystyle T}}_3& {{\displaystyle T}}_4\end{array}\right]$ (8)

式中，T为管路传递矩阵，T_i(i = 1, 2, 3, 4)为4 × 4矩阵。

再由式即可得到弯管中离散直管全局坐标系下的阻抗矩阵，管路两端满足力和位移的平衡，对弯管系统进行组装，得到弯管系统力与位移的关系为：

$\left[\begin{array}{c}{{\displaystyle F}}_A\\ {{\displaystyle F}}_B\\ {{\displaystyle F}}_E\\ {{\displaystyle F}}_F\\ {{\displaystyle F}}_C\\ {{\displaystyle F}}_D\end{array}\right]=\left[\begin{array}{cccccc}{{\displaystyle Z}}_1^{AB}& {{\displaystyle Z}}_2^{AB}& 0& 0& 0& 0\\ {{\displaystyle Z}}_3^{AB}& {{\displaystyle Z}}_4^{AB}+{{\displaystyle Z}}_1^{BE}& {{\displaystyle Z}}_2^{BE}& 0& 0& 0\\ 0& {{\displaystyle Z}}_3^{BE}& {{\displaystyle Z}}_4^{BE}+{{\displaystyle Z}}_1^{EF}& {{\displaystyle Z}}_2^{EF}& 0& 0\\ 0& 0& {{\displaystyle Z}}_3^{EF}& {{\displaystyle Z}}_4^{EF}+{{\displaystyle Z}}_1^{FC}& {{\displaystyle Z}}_2^{FC}& 0\\ 0& 0& 0& {{\displaystyle Z}}_3^{FC}& {{\displaystyle Z}}_4^{FC}+{{\displaystyle Z}}_1^{CD}& {{\displaystyle Z}}_2^{CD}\\ 0& 0& 0& 0& {{\displaystyle Z}}_3^{CD}& {{\displaystyle Z}}_4^{CD}\end{array}\right]\left[\begin{array}{c}{{\displaystyle W}}_A\\ {{\displaystyle W}}_B\\ {{\displaystyle W}}_E\\ {{\displaystyle W}}_F\\ {{\displaystyle W}}_C\\ {{\displaystyle W}}_D\end{array}\right]$ (9)

式中，F包括力和力矩，W包括位移和转角， $Z_j^i$ (i = 1, 2, 3, 4; j = AB, BE, EF, FC, CD)为4 × 4矩阵，弯管系统的阻抗矩阵为24 × 24矩阵。

2.2. ISM-FEM混合方法计算管路–弹性板耦合系统

应用FEM分别对弹性板的内环施加F_z、M_x、F_y单位激励力，得到管路弹性板耦合处位移阻抗边界条件，结合管路ISM即可对整个管路–弹性板耦合系统的振动位移响应进行求解。弹性板为薄板，忽略管路与弹性板耦合处(z = L)管壁位移与流体声压的耦合，根据阻抗矩阵的定义得到弹性板内环的位移阻抗边界条件为：

$\left[\begin{array}{c}{P}^p\left(L\right)\\ F_z^p\left(L\right)\\ M_x^p\left(L\right)\\ F_y^p\left(L\right)\end{array}\right]=\left[\begin{array}{cccc}0& Z_{12}^f& 0& 0\\ Z_{21}^p& 0& 0& 0\\ 0& 0& Z_{33}^p& Z_{34}^p\\ 0& 0& Z_{43}^p& Z_{44}^p\end{array}\right]\left[\begin{array}{c}{U}_z^p\left(L\right)\\ U_f^p\left(L\right)\\ U_y^p\left(L\right)\\ {\Phi }_x^p\left(L\right)\end{array}\right]$ (10)

将式简写为：

$F^p\left(L\right)=Z^p{W}^p\left(L\right)$ (11)

式中 $F^p$ 和 $W^p$ 为弹性板内环状态向量，元素 ${{\displaystyle Z}}_{ij}^p$ 为管路与弹性板耦合处的结构输入位移阻抗， ${{\displaystyle Z}}_{12}^f$ 是管口处流体的声辐射位移阻抗，由下式给出：

$Z_{12}^f=\frac{{\rho }_0{c}_0}{j\omega }\left[1-\frac{{{\displaystyle 2J_1}}_{}\left({{\displaystyle 2k{R}_1}}_{}\right)}{2k{R}_1}+j\frac{{{\displaystyle 2K_1}}_{}\left({{\displaystyle 2k{R}_1}}_{}\right)}{{{\displaystyle \left({{\displaystyle 2k{R}_1}}_{}\right)}}^2}\right]e^{j\omega t}$ (12)

式中，R₁为弹性板内径，k为波数，J₁为一阶贝塞尔函数，K₁为一阶修正贝塞尔函数，本文算例中管口外部的流体区域为轻流体，ρ₀和c₀分别为空气密度和空气中声速。

在管路和弹性板耦合处满足力和位移的连续条件：

$F\left(L\right)=-F^p\left(L\right)$ ， $W\left(L\right)=W^p\left(L\right)$ (13)

将力的连续性边界条件式代入式得：

$\left[\begin{array}{c}F\left(0\right)\\ -F^p\left(L\right)\end{array}\right]=\left[\begin{array}{cc}{Z}_1& Z_2\\ Z_3& Z_4\end{array}\right]\left[\begin{array}{c}W\left(0\right)\\ W\left(L\right)\end{array}\right]$ (14)

由将式代入式，由式位移边界条件的连续性得：

$W\left(0\right)=Z_3^{-1}\left(-Z^p-Z_4\right)W\left(L\right)$ (15)

由式已知：

$F\left(0\right)=Z_{1}W\left(0\right)+Z_{2}W\left(L\right)$ (16)

将式代入式得到管板连接处响应为：

$W\left(L\right)={\left(Z_1{Z}_3^{-1}\left(-Z^p-Z_4\right)+Z_2\right)}^{-1}F\left(0\right)$ (17)

式表明，一旦管道一端激励已知后，可以求得管路与弹性板耦合处的振动位移响应及管路出口流体的振速。

3. 计算结果及分析

为验证本文ISM-FEM混合计算方法的正确性，与直接FEM计算结果对比分析，两种管路–弹性板耦合系统有限元模型如图3所示。直管、弯管总长度相同，管路和弹性板材料均为20号钢，具体尺寸和材料参数如表1所示，有限元软件中管壁和弹性板均采用壳单元，管内流体采用声学单元，网格尺寸均为0.02 m，流体和结构绑定接触，管路弹性板刚性连接，管口外部流体区域为轻流体，水的流体声学体积模量和密度分别为K = 2.14 × 109 Pa和ρ_f = 1000 kg/m³，水中声速为c₀ = 1500 m/s，计算频率为1~1000 Hz，分析步长为1.0 Hz。

Table 1. Material and dimensional parameters of two piping-elastic plate coupling systems

表1. 两种管路–弹性板耦合系统的材料和尺寸参数

Figure 3. Two models of coupled pipe-elastic plate systems (FEM direct calculation)

图3. 两种管路–弹性板耦合系统模型(FEM直接计算)

3.1. 阻抗综合法计算管路振动响应

采用ISM计算两种管路–弹性板系统中的直管和弯管在不同边界条件下的振动位移响应。直管两端标为P₁和P₂如图1所示，P₁端为固支边界条件，管内充满水介质，分别对P₂端管壁和流体端面施加轴向力激励(F_z = 1 N，相位为0)和平面波激励(P = 1 Pa，相位为0)，P₂端振动响应如图4(a)~(b)所示；弯管两端标为P₁和P₂如图2所示，两端为自由边界条件，管内空气，对管路P₁端施加轴向力激励(F_z = 1 N，相位为0)，P₂端振动响应如图4(c)~(d)所示。ISM与FEM在本文算例中计算结果基本一致，证明了ISM-FEM混合计算方法研究管路-弹性板耦合系统的声振传递机制中，ISM预报管路振动响应的正确性。

Figure 4. Comparison of P₂ vibration displacement response results for straight and bent pipes: (a) axial force excitation of straight pipe P₁ wall P₂ axial displacement; (b) plane wave excitation of straight pipe P₁ wall P₂ axial displacement; (c) axial force excitation of elbow pipe P₁ wall P₂ axial displacement; (d) axial force excitation of elbow pipe P₁ wall P₂ lateral displacement

图4. 直管和弯管P₂振动位移响应结果对比：(a) 轴向力激励直管P₁管壁P₂轴向位移；(b) 平面波激励直管P₁管壁P₂轴向位移；(c) 轴向力激励弯管P₁管壁P₂轴向位移；(d) 轴向力激励弯管P₁管壁P₂横向位移

3.2. ISM-FEM混合计算方法验证及结果分析

对直管–弹性板耦合系统的管路与弹性板耦合处的振动响应求解分析。弹性板内环阻抗矩阵由FEM求得，即得到了管路与弹性板耦合处的阻抗边界条件；直管部分则采用ISM求得管路阻抗矩阵，由2.2节中理论推导即可求得直管–弹性板系统整体的阻抗矩阵，表明在管路入口处给定激励，即可求得管路与弹性板耦合处的振动响应。

弹性板外边界钳定，管路P₁端为自由边界，在管路P₁端分别施加轴向力激励(F_z = 1 N，相位为0)和平面波激励(P = 1 Pa，相位为0)。管路与弹性板耦合处的振动响应如图5所示，从图中可以看出本文提出的ISM-FEM混合计算方法均与FEM计算结果非常吻合，验证了本方法的正确性。

Figure 5. Comparison of P₂ vibration displacement response results for straight and bent pipes: (a) axial force excitation of straight pipe P₁ wall P₂ axial displacement; (b) plane wave excitation of straight pipe P₁ wall P₂ axial displacement; (c) axial force excitation of elbow pipe P₁ wall P₂ axial displacement; (d) axial force excitation of elbow pipe P₁ wall P₂ lateral displacement

图5. 直管–弹性板系统振动响应对比：(a) 轴向力激励P₁管板耦合处P₂轴向位移；(b) 轴向力激励P₁管路出口处P₂流体振速；(c) 平面波激励P₁管板耦合处P₂轴向位移；(d) 平面波激励P₁管路出口处P₂流体振速

对弯管–弹性板耦合系统的管路与弹性板耦合处的振动响应求解分析。管路与弹性板耦合处的位移阻抗边界条件在计算直管–弹性板耦合系统时已知，同理弯管部分采用ISM求得管路阻抗矩阵，即可求得弯管–弹性板系统整体的阻抗矩阵。

Figure 6. P₂ axial displacement at the pipe-elastic plate system coupling

图6. 管路–弹性板系统耦合处P₂轴向位移

同样弹性板外边界钳定，管路P₁端为自由边界，在管路P₁端管壁施加轴向力激励(F_z = 1 N，相位为0)，两种管路–弹性板系统管内分别充水和空气，管路与弹性板耦合处的轴向位移响应如图6所示。带有弹性板的充液弯管系统与空气管道相比，在1000 Hz范围内，由于管路结构与管内流体耦合，管内流体质量增加，充液管路固有频率和振动频响曲线向低频偏移。弯管–弹性板耦合系统与直管–弹性板耦合系统相比相同阶数的固有频率显著降低，共振峰的数量增加，整体振动水平有所降低。

Table 2. Comparison of the efficiency of ISM-FEM hybrid computational method with direct FEM computation

表2. ISM-FEM混合计算方法与直接FEM计算效率对比

两种计算方法计算效率如表2所示，直管–弹性板耦合系统计算时长为管内充液，分别施加轴向力、平面波激励求解振动位移响应总时长，弯管–弹性板耦合系统计算时长为管内分别充液和空气，施加轴向力激励求解振动位移响应总时长。结果表明，在本算例中，ISM-FEM混合计算方法约比直接FEM快11倍，在研究末端连接同种弹性板或管路不同参数影响时，ISM-FEM混合计算方法仅需建立弹性板有限元模型即可求解不同参数充液管路与弹性板耦合系统的振动响应，无需重复建立管路–弹性板耦合系统有限元模型，占用计算资源少，计算效率优势会更加显著。

4. 结论

本文采用ISM-FEM混合计算方法，研究管路与壳体的声振耦合的求解方法和传递机制。建立两种管路–弹性板耦合系统模型，通过算例从充液管路与弹性板耦合处管壁的轴向位移、管内流体的振速进行分析。计算结果与FEM对比非常吻合，证明了本文方法的正确性，且计算效率更高，尤其是在计算不同参数影响或连接同种弹性板的充液管路系统时有显著优势。并分析了弯管和管内流体对系统声振耦合的影响，结果表明，管路–弹性板耦合系统管内充液和空气相比，管内流体质量会使充液管路–弹性板系统固有频率和振动频响曲线向低频偏移。弯管使得管路–弹性板耦合系统比直管相同阶数的固有频率显著降低，共振峰的数量增加，整体振动水平有所降低。管口外部的流体区域为重流体时，与壳体耦合的充液管路系统声振耦合传递机制更为复杂，更进一步工作将集中在管口外部重流体负载时，空间管路弹性板耦合系统的声振传递机制和辐射噪声预报。

基金项目

国家自然科学基金(12374447)。

NOTES

^*通讯作者。

参考文献