张子祥，王检耀，王鸿东，易 宏

(1.上海交通大学 海洋工程国家重点实验室，上海 200240；2.上海交通大学 海洋智能装备与系统教育部重点实验室，上海 200240)

充液管道在舰船内部分布广泛，其振动所产生的噪声对舰船战斗力有重大的影响：① 管道振动经舰船通海口向外部辐射噪声，是除了螺旋桨噪声、水动力噪声和设备基座噪声以外主要的辐射噪声源[1]；② 管道振动诱发的频率与声纳设备的工作频率相近，降低舰船的探测能力，是舰船自噪声的主要来源[2]。因此对充液管道开展振动特性的研究具有重要意义。

如果管道内没有充满液体而完全被空气包围，空气对结构的作用通常可以忽略。这类不考虑外部流体对结构振动特性影响，求解结构固有频率和模态的问题称为干模态问题。对充液管道而言，壁面上流体(如水、油等)与固体结构的耦合现象通常不可忽略，耦合形式主要包括摩擦耦合、泊松耦合和连接耦合等[3-4]，一般表现为管道固有频率相对干模态的下降。这类考虑流固耦合效应求解结构固有频率和模态的问题称为湿模态问题[5]。对流固耦合的处理方式包括附加质量法、单向流固耦合和双向流固耦合等。附加质量法[6]一般直接将流体在单位长度上的质量等价为管道在单位长度上的附加质量，周名德等[7]考虑流固耦合作用给出了两端简支单管附加质量的解析解。单向流固耦合[8]忽略管壁上的位移协调，首先计算流体对静止管道内壁的作用力，再将充液管道视为有预应力作用的未充液管道进行仿真计算。由于忽略了位移协调条件，当管道刚度较低或激励幅值较大造成管道发生较大位移时，计算与试验结果的误差较大[9]。双向流固耦合同时考虑管壁上的位移协调条件和力平衡条件，通过直接求解流固耦合方程组或流体与固体计算迭代收敛的方法[10]进行求解，是最符合实际情况的流固耦合处理方式。充液管道湿模态的双向流固耦合数值计算方法主要包括传递矩阵法和有限元法。传递矩阵法[11-17]将充液管道离散为若干管段，在联立管道和流体在横向、轴向和周向运动微分方程组的基础上建立相邻管段的传递矩阵，计算效率相对直接有限元法高，但对弯管、变截面管和复杂管系等的处理比较困难。有限元法对分析对象的要求较传递矩阵法低，可以适应各种形状的管道和各种形式的管系，因此无论在理论研究还是工程应用上都更具现实意义。类似于传递矩阵法，一些学者提出了专用于管道分析研究的有限元模型，在计算结果基本相同的条件下，比直接离散化的有限元模型仿真计算成本更低。Hansson等[18]将轴对称的管道壳单元和一维可压无黏流体单元相结合，建立了专用于研究充液管道的有限元模型单元，并与ABAQUS的仿真结果进行了比较。Seo等[19]建立了与理想流体双向耦合的薄圆柱壳单元，用于进行无流速管道和有流速管道的模态分析，并与试验结果和NASTRAN的仿真结果进行了比较。另一些学者将充液管道按固体单元和流体单元分别直接离散，基于有限元法进行复杂管道或管系的双向流固耦合湿模态分析。Li等[20]对飞机液压管道中的充液直管和充液弯管进行了有限元仿真和试验，仿真与试验得到的固有频率吻合良好。黄益民等[21]基于对充液直角弯管的有限元仿真，进行了结构参数对固有频率的灵敏度分析和防共振可靠性计算。Makaryants等[22]对整个冲床液压管系的振动特性进行了有限元仿真，指出了管系结构失效的原因和相应改进措施。姜峰等[23]对浸没在海水中且内部充满石油的海洋立管的振动特性进行有限元仿真，并在此基础上开展了结构设计优化。

为了进一步实现充液管道的减振降噪，常常通过调整管道支承的刚度，实现避免结构共振和减少振动传递到船体的能量[24]。然而这类管道支承的实际装配效果与理论设计可能存在偏差，根据理论设计建立的有限元模型并不能反映实际情况。因此在对充液管道的振动特性进行分析前，需要对这类弹性约束进行修正。目前这类问题的研究对象主要集中于桥梁结构，以充液管道为对象的研究较少，但利用优化法修正有限元模型的思想是基本相同的。以仿真与试验的误差为目标函数，以弹性约束的值作为决策变量，弹性约束的修正就转换为了一般的优化问题[25]。通过仿真与试验分别能得到多阶模态参数，因此仿真与试验之间的各阶误差可以组成一个误差向量。一般或是出于实际工程需要，或是受限于优化手段，通常将该误差向量各分量的平方和作为优化的目标函数，于是原问题简化为一个单目标优化问题。Cao等[26]为了区分对不同阶数的误差要求，在误差向量中引入权重系数，使用优化方法得到了桥梁边界的连接刚度。施洲等[27]在最小二乘优化的基础上，为提升计算效率引入矩阵摄动，对桥梁的附加支撑刚度进行了修正。邹向农等[28]为得到桥梁有限元模型的纵向约束刚度等参数，采用逐步回归法得到桥梁固有频率的响应面代理模型，再通过基于遗传算法的单目标优化进行模型修正。一般而言，单目标优化较多目标优化的操作更简单，但存在难以确定合适的权重，当误差中有多种物理量时目标函数的物理意义不清晰等问题。因此也有部分研究者采取各种手段直接求解多目标优化问题，力图使误差向量的各个分量尽可能最小化，如Ntotsios等[29]通过将求解原多目标优化问题转化为求解大量权重分布不同的单目标规划问题，实现桥台、桥墩等位置约束条件的修正。Jaishi等[30]采用目标规划法，通过为各个优化目标指定目标值和权重，将原多目标规划问题转化为单目标规划问题并实现对各个目标的同时优化，进行桥台刚度、桥墩刚度等参数的修正。总体而言，使用多目标优化方法对充液管道未知约束进行修正的研究较少，有待发掘更多的方法进行探索。

本文开展了弹性约束充液管道的干湿模态试验研究，并基于双向流固耦合的有限元法和多目标遗传算法NSGA-Ⅱ提出了对未知弹性约束充液管道频域特性的预报方法，分析流程如图1所示。首先进行管道干模态试验获得管道的固有频率和振型，并建立相应的有限元模型。然后以未知的弹性约束为决策变量，以仿真与试验固有频率的误差向量为目标函数，建立多目标优化问题，并基于NSGA-Ⅱ进行求解并修正弹性约束。最后将该弹性约束代入充液管道的湿模态有限元模型求解，比较湿模态仿真与试验的固有频率和振型，验证上述预报方法的可行性并分析管道充液前后的振动特性变化。

图1 分析流程

1 有限元建模理论

1.1 干模态理论

结构振动的有限元方程可以表达为[31]：

(1)

对干模态问题，求解激励力向量{FS}为零时结构的运动微分方程，即结构的自由振动方程。一般情况下结构的阻尼较小，对固有频率和固有振型影响不大，所以常常略去总体阻尼矩阵，干模态问题转化为求解：

(2)

1.2 湿模态理论

对湿模态问题，本文假设周围流场中的流体满足以下要求：① 流体为理想流体；② 流场无旋；③ 流场没有宏观流动；④ 流场中无体积力；⑤ 流场中的振动幅度非常小。

此时流体单元内的等效积分方程表达式为

(3)

式中：c为声波在流体介质中的速度；p为声压；ρF为流体密度；ΓFSI为流固耦合面；nFSI为ΓFSI上某点的法向量。

将计算域离散化，各个单元的声压和位移都可以用其节点的声压与位移插值表示：

p=[N]{p}(e)

(4)

u=[N′]{u}(e)

(5)

式中：{p}(e)、{u}(e)分别为各个单元每一个节点的声压和位移组成的向量；[N]、[N′]分别为声压和位移的形函数阵。

将式(4)和式(5)代入式(3)，可以得到流场单元的波动方程：

(6)

其中:

将式(6)中的单元矩阵叠加，形成流体的总体波动方程：

(7)

流固耦合面ΓFSI上的结构单元受到流体声压的作用，单元节点上作用力的大小为：

[R(e)]{p}(e)

(8)

将式(8)中的单元矩阵叠加，形成结构的总体声压作用力的表达式：

{Fpr}=∑[R(e)]{p}(e)=[R]{p}

(9)

则结构的总体双向流固耦合振动方程可以表达为：

(10)

综上所述，湿模态问题即求解式(7)和式(10)联立的双向流固耦合方程组[32]：

(11)

1.3 弹性约束的处理

为了简化问题，本文假设约束的主要成分是径向的、作用在管道与约束接触面上的、均匀的弹性支承约束，忽略其它约束作用。

管道与约束接触面ΓES上的结构单元受到弹性支承约束的作用，单元节点上作用力的大小为：

(12)

式中，k为管道与约束接触面ΓES单位面积上产生单位长度位移所形成的弹性支承约束力。

将式(12)中的单元矩阵叠加，形成结构的总体弹性约束作用力的表达式：

(13)

在这里，将干燥管道与弹性约束的综合称为干模态系统，将充液管道与弹性约束的综合称为湿模态系统。本文假设无论是干模态系统还是湿模态系统，两种系统中对应的弹性约束都没有发生变化。这是因为：① 没有对弹性约束进行替换或调整松紧程度，弹性约束的刚度作为其固有属性没有发生变化；② 没有对弹性约束的位置进行调整，弹性约束对管道的作用位置没有发生变化。所以，干模态系统和湿模态系统中总体弹性约束作用力的表达式是完全相同的。将式(13)分别代入式(2)和式(11)，得到干模态问题和湿模态问题的完整表达：

(14)

(15)

2 试验方案与装置

由振动理论，频响函数矩阵[H(jω)]第p行的各个元素分别为p号测点对各个对应测点处激励的频响函数，表达式为

(16)

式中：ω为输入的激励频率；N为测点数；ζr为模态阻尼比；ωr为第r阶模态的固有频率；{U}r={u1r,u2r,…,uNr}T为质量归一化后的第r阶模态振型。显然，只要获得[H(jω)]某一行的所有元素，就可以将整个频响函数矩阵确定下来。

本文为了获得管道试件在干模态和湿模态下的固有振型和频率，使用锤击法进行模态试验，采用移动敲击力锤，固定加速度传感器的方式进行：敲击力锤依次激励管道外表面上布置的测点，动态信号采集分析系统同时采集敲击力锤的输入信号与加速度传感器的输出信号，可以得到所有锤击测点与传感器测点间的频响曲线。各阶固有频率ωr和模态阻尼比ζr通过对频响曲线进行拟合获得，代入式(16)即可得到各阶振型{U}r。

用于试验的管道试件为316L钢质直管。直管的外径为22 mm，内径为18 mm，长度为1 150 mm。316L钢的密度为8 027 kg/m3，杨氏模量为206 GPa，泊松比为0.3。管道内部的流体为水，密度为998.2 kg/m3，声速为1 482.1 m/s。

试验系统由计算机、LMS 动态信号采集分析系统、三向加速度传感器、敲击力锤、流体循环系统以及试件工装等组成。试验系统简图见图2。试验现场见图3。其中：约束的基座固定在桌面上，使用可调节松紧的箍套实现对管道试件的约束。进行湿模态试验时，关闭水阀并将水槽灌满水，稳定后即可实现流体静止的充液管道。

图2 试验系统简图

图3 试验现场

试件上每间隔50 mm在上、下、前、后各布置一个测点，总计设置84个测点：沿轴向从右向左编号依次加1，沿周向按前、上、后、下的顺时针顺序编号依次加21。测点布置情况示意见图4。三向加速度传感器布置在23号测点的位置，约束A和B分别以20号测点和11号测点所在管道截面为中横剖面布置。

图4 测点的布置

3 干模态仿真计算

3.1 干模态仿真物理模型

建立管道试件的物理模型，见图5。管道试件的模型由7段组成。管道左右端部的2段为连接软管，字母A和B处为约束管道的箍套。忽略管道两端的连接软管对振动的影响，将箍套对管道的约束简化为径向均匀的弹性约束。约束A的左端距管道左端75 mm，约束B的左端距管道左端570 mm。两处弹性约束的长度均为10 mm，刚度未知。使用壳单元进行网格划分，单元数为1 870，节点数为1 880。频率分析范围为0～1 600 Hz。

图5 干模态的物理模型

3.2 多目标优化理论

为了找到最接近真实弹性约束的弹性约束刚度，以管道试件固有频率的仿真结果和试验结果的误差最小化为导向，建立最优化问题：

(17)

式中：k为A和B处的弹性约束刚度组成的向量；K为定义域；ffem(k)和fexp(k)分别为管道试件通过仿真和试验得到的固有频率组成的向量；E(k)为误差向量。

对多目标优化问题，一般难以找到令所有目标同时最优的可行解。因此需要引入非支配解的概念。

对优化问题的可行解k*∈K,若不存在k′∈K且k′≠k*，使目标函数E(k′)满足[33]：

E(k′)≤E(k*)

(18)

且至少有一个严格不等式成立，则称k*为该多目标优化的非支配解。若干非支配解k*形成的曲面称为非支配前沿面，要找的弹性约束刚度就在该前沿面上。

利用多目标遗传算法NSGA-Ⅱ求解上述多目标优化问题。遗传算法的基本原理为通过选择、交叉和变异等手段生成和保留适应度高的个体。对多目标优化问题，适应度的定义方法是在所有个体中寻找所有非支配解，形成第一层非支配前沿面，再寻找剩下个体中的所有非支配解并形成下一层非支配前沿面，直到所有个体被分类完毕。层级越低的个体被分配到越高的适应度。NSGA-Ⅱ的特点包括[34]：① 相比NSGA算法，NSGA-Ⅱ采用的快速非劣排序算法速度更快；② 采用精英策略将父代个体组合入子代，避免损失优秀个体，加速算法收敛；③ 依据拥挤度对同一前沿面上的个体进行排序，无需额外指定参数即可确保个体的分散和多样性。NSGA-Ⅱ的算法流程图见图6。

图6 多目标遗传算法NSGA-Ⅱ流程图

经过多目标遗传算法NSGA-Ⅱ得到非支配前沿面后，根据实际工程经验选择合适的非支配解作为优化结果，修正弹性约束。

3.3 弹性约束刚度优化

干模态振动试验得到的管道试件固有频率见表1的序号A0。当弹性约束刚度无限大时，即直接在约束面上加载理想约束，使约束面不发生位移和变形时，固有频率的有限元仿真结果见序号B0。联合干模态有限元仿真与NSGA-Ⅱ优化得到的弹性约束刚度的部分方案见序号C1到C6。

表1 弹性约束刚度部分优化结果

如果直接使用理想约束进行仿真，干模态固有频率的误差最多可以超过100%，且阶数越低误差越大。而经过NSGA-Ⅱ优化后，非支配前沿面上的刚度方案有效地降低了各阶固有频率的误差水平。其中，C2方案的平均误差δmean为3.06%，是所有方案中较小的。

干模态试验与仿真得到的固有频率对比见图7。

图7 干模态试验与仿真固有频率对比

对比结果显示，C2方案的固有频率误差最高为6.19%，出现在第2阶；最低为0，出现在第4阶。

比较干模态振动试验与C2方案仿真得到的固有振型见表2，可以观察到两者各阶固有振型的形状均相似。综合考虑固有频率和固有振型的吻合程度，采用C2的弹性约束刚度方案，作为修正后的弹性约束开展后续的湿模态分析。

表2 干模态试验与仿真的固有振型

4 湿模态仿真计算

4.1 湿模态仿真物理模型

湿模态计算基于干模态的物理模型与多目标优化算法NSGA-Ⅱ得到的弹性约束刚度，在管道内部建立流场域并生成流固耦合面。湿模态的物理模型见图8。

图8 湿模态的物理模型

使用流体单元对流场域进行网格划分，单元数为39 116，节点数为181 006。频率分析范围为0～1 600 Hz。

4.2 仿真与试验结果分析

湿模态振动试验的试验与仿真结果见表3和表4。

表3 湿模态试验与仿真的固有频率

表4 湿模态试验与仿真的固有振型

湿模态试验与仿真得到的固有频率对比见图9。对比结果显示：采用C2的弹性约束刚度方案，固有频率误差最高为5.30%，出现在第4阶；最低为0.09%，出现在第3阶。

图9 湿模态试验与仿真固有频率对比

干模态与湿模态仿真结果的固有频率与试验结果的平均误差分别为3.06%和2.47%，最大误差分别不超过6.19%和5.30%，可以基本满足实际要求，验证了本文提出方法的可行性。造成误差的主要原因包括：① 忽略了箍套的非线性弹性和阻尼作用；② 弹性支承对管道的作用被简化为是沿管道周向和轴向均匀分布的，忽略了弹性作用的不均匀性。

干模态和湿模态试验得到的固有频率对比见图10。试验结果显示，相同阶的湿模态固有频率均低于干模态固有频率，频率漂移比例的平均值为-6.68%。同时，不同阶数频率偏移的波动非常小，没有明显的增加或减少趋势，基本可以认为在一定范围内频率偏移基本不变。因此，在选择弹性约束方案时，若存在重点关注的阶数或固有频率范围，应该优先选择该阶数下干模态仿真误差小的弹性约束方案，有利于在该阶数上得到更符合实际的湿模态固有频率。

图10 试验固有频率对比

4.3 本文方法的合理性验证

湿模态与干模态的各阶模态存在着对应关系：在相似的固有振型下，湿模态的固有频率相对于干模态发生漂移。所以如果忽略试验误差和仿真误差等因素，同一刚度方案在干模态下对不同阶固有频率误差的相对差别，不同刚度方案在干模态下对同一阶固有频率误差的相对差别都可以通过这种对应关系传递到湿模态上。

观察干模态与湿模态各阶频率平均误差的相关性，从实际计算结果验证本文提出方法的合理性：选取在弹性约束刚度优化过程中计算过的刚度方案，对这些刚度方案下的干湿模态固有频率进行仿真。以干模态固有频率平均误差为横坐标，湿模态固有频率平均误差为纵坐标如图11所示。

图11 干湿模态固有频率平均误差的关系

图11显示干模态固有频率平均误差与湿模态固有频率平均误差近似呈正比例的线性关系：干模态固有频率平均误差越低的刚度方案，其湿模态固有频率的误差水平也越低。从该图可以推论C2方案应用于湿模态预报的可靠性，也验证了本文方法的合理性。

5 结 论

本文开展了弹性约束充液管道的干湿模态试验研究，并提出了对未知弹性约束充液管道频域特性的预报方法。基于文中的假设，该方法总体适用于管道结构的阻尼较小，管内流体无流速或流速较低且约束的主要成分为弹性约束的情况。仿真与试验结果表明：

(1)弹性约束经过多目标遗传算法NSGA-Ⅱ修正后，干模态有限元仿真结果与试验结果的固有频率平均误差为3.06%，最大误差在6.19%以内，振型基本相同。

(2)基于经过修正的弹性约束，湿模态有限元仿真结果与试验结果的固有频率平均误差为2.47%，最大误差在5.30%以内，振型基本相同，验证了本文预报方法的可行性与合理性。

(3)湿模态的固有频率较干模态有明显降低，且各阶频率漂移比例基本相同，约为-6.68%。若存在重点关注的阶数或固有频率范围，应该优先选择该阶数下干模态仿真误差小的弹性约束方案，以获得更加准确的湿模态固有频率和振型。