基于分片传递函数法的声振耦合

金鹏飞，史治宇，向育佳; JIN Pengfei; SHI Zhiyu; XIANG Yujia

2025年5月2日 6:22 星期五

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

基于分片传递函数法的声振耦合 PDF

- ORCID：
金鹏飞 ¹
- ORCID：
史治宇 ¹
✉
- ORCID：
向育佳 ²

1. 南京航空航天大学航空航天结构力学及控制全国重点实验室, 南京 210016； 2. 上海索辰信息科技股份有限公司, 上海 201204

中图分类号： TB535

最近更新：2024-08-19

DOI：10.16356/j.1005-2615.2024.04.020

摘要

强耦合条件下，复杂系统的有限元分析受到自由度的影响计算成本高昂。虽然模态叠加技术能减少耦合系统的自由度，但基于忽略不同子系统间高阶模态和低阶模态耦合作用的假设，若使用非耦合模态可能导致收敛性差。基于无阻尼声振耦合方程，本文采用分片传递函数法（Patch transfer function method， PTFM）将耦合面处理成一系列独立的分片，每个分片上所属单元节点值的平均值定义为分片传递函数，通过使用连续性关系计算耦合系统的分片传递函数。将耦合矩阵的求逆简化为源点到响应点的传递函数，可以快速计算耦合系统响应值。分别利用直接耦合法（Direct coupling method， DCM）和分片传递函数法计算板和空气声腔耦合模型，验证分片传递函数法的有效性，并讨论了分片传递函数的算法原理和计算误差。

关键词

分片传递函数法; 声振耦合; 频率响应; 有限元分析

过去几十年的声振耦合研究中，有限元法^［

1］、有限元‑边界元法^{［参考文献 2

百度学术}2］、波函数法^{［参考文献 3

百度学术}3］等直接构造耦合方程的数值方法已经发展得十分成熟，然而数值算法的计算精度和试验值之间仍然存在差异^{［参考文献 4

百度学术}4］，主要有两方面原因，一是耦合模型的计算规模，二是耦合面的处理。对于复杂的声源边界以及声腔，为了获得精确的计算结果，可以通过加密网格来实现，但随之而来的问题是自由度增加，计算规模庞大，并且方程系数矩阵不对称，不能使用稀疏对称性提升计算效率。在工程中分析复杂声振耦合系统的数值方法一般会使用模态叠加法^{［参考文献 5

百度学术}5］降低模型的计算规模，但在分析频段范围内有效的模态叠加计算需要取合适数量的正交模态作为基向量，且产生的截断误差会影响计算精度^{［参考文献 6

百度学术}6］。另一方面的误差源于对耦合面的处理，使用有限元方法通常需对耦合方程组进行解耦。Kim等^{［参考文献 4

百度学术}4］提出了基于模态叠加理论和耦合系数的概念实现了方程的解耦。然而，该方法的实质是单向解耦，必然导致系统性偏差。分片传递函数（Patch transfer function， PTF）是一种子系统划分方法^{［参考文献 7

百度学术}7］，通过线性假设和叠加原理实现耦合方程的求解，但该方法在子系统划分^{［参考文献 8

百度学术}8］、误差来源与收敛性分析^{［参考文献 7

百度学术}7，9］等方面仍有待探索。

本文针对无阻尼声振耦合系统，将耦合面处理成若干独立分片，每个分片的计算域内分别包含不同的声学和结构单元。通过所属分片的节点声压均值和速度均值的连续性关系定义从结构激励源到结构响应点的分片传递函数，进而解耦计算得到耦合系统的响应。最后，以一个四边简支板和空气声腔耦合模型作为分析对象，对比分片函数传递法（Patch transfer function method， PTFM）、文献［

4］方法和直接耦合法（Direct coupling method， DCM）结果以验证PTFM的有效性，并应用算例讨论PTFM在有限元方法中求解频响函数的误差来源。

1 模型建立

如图1所示，弹性结构体Ω_s和声腔Ω_a组成声振耦合系统，结构与腔体耦合面为S_c。

图1 声振耦合系统

Fig.1 Acoustic-vibration coupling system

在Ω_s上施加简谐激励，结构体Ω_s的振动微分方程和声场域Ω_a的亥姆霍兹方程为

([\begin{matrix} K_{s} & - A \\ 0 & K_{a} \end{matrix}] - ω^{2} [\begin{matrix} M_{s} & 0 \\ ρ_{0} A^{T} & M_{a} \end{matrix}]) [\begin{matrix} U_{s} \\ P \end{matrix}] = [\begin{matrix} F \\ 0 \end{matrix}]

（1）

式中：F为施加在结构上的简谐力幅值，U_s表示弹性结构位移幅值，ρ₀为声场介质密度，K_s、K_a和M_s、M_a分别为经有限元法离散结构和声场的刚度矩阵和质量矩阵^［

1］，A表示声场对结构的耦合矩阵^{［参考文献 1

百度学术}1］。若结构和声场形函数分别为N_s和N_a，耦合面外法线向量为n，则

A = \iint_{S_{c}}^{} N_{s}^{T} n N_{a} d S

（2）

\{\begin{array}{l} M_{a} = ∭_{Ω_{a}}^{} N_{a}^{T} N_{a} d V \\ K_{a} = c_{0}^{2} ∭_{Ω_{a}}^{} (\nabla N_{a})^{T} \nabla N_{a} d V \\ M_{s} = ∭_{Ω_{s}}^{} N_{s}^{T} N_{s} d V \\ K_{s} = \frac{1}{ρ_{s}} ∭_{Ω_{s}}^{} (\nabla N_{s})^{T} D \nabla N_{s} d V \end{array}

（3）

式中：微分算子 $\nabla = {[\begin{matrix} \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \end{matrix}]}^{T}$ ， $c_{0} = \sqrt[]{\frac{K}{ρ_{0}}}$ 表示声速，K为绝热体积弹性系数，ρ_s为弹性结构体密度，D为结构体弹性模量矩阵。

2 分片传递函数法

如图2所示，将整个声振耦合系统独立为两个系统各自离散。两个系统间仅通过若干分片利用速度连续性和应力边界建立耦合场传递函数。

图2 分片传递函数子系统示意图

Fig.2 Diagram of PTF subsystem

结构和声场在耦合面的任意分片j上的平均速度和平均声压满足连续性条件，即

\{\begin{matrix} {\bar{p}}_{j}^{s} = {\bar{p}}_{j}^{a} = {\bar{p}}_{_{j}}^{c} \\ {\bar{v}}_{j}^{s} = {\bar{v}}_{j}^{a} = {\bar{v}}_{_{j}}^{c} \end{matrix}

（4）

式中： ${\bar{(\cdot)}}_{j} = \frac{1}{S_{j}} \iint_{S_{j}}^{} (\cdot) d S$ 表示面域平均计算式；S_j为分片面积；上标a代表声场、s代表结构、c代表耦合面，即分片上的交接面。

对于结构，第j个分片受到第k个分片的声压激励产生平均速度响应，根据导纳 $Y_{j k}^{s}$ 线性关系，有

{\bar{v}}_{j}^{s} = Y_{j k}^{s} \cdot {\bar{p}}_{k}^{s}

（5）

同时对于声场，有

{\bar{p}}_{j}^{a} = Z_{j k}^{a} \cdot {\bar{v}}_{k}^{a}

（6）

式中 $Z_{j k}^{a}$ 表示声场阻抗。因此，耦合系统结构上的平均速度响应 ${\bar{v}}^{s}$ 为未耦合系统结构上平均速度响应 ${\tilde{\bar{v}}}^{s}$ 和耦合时因声场作用产生的附加平均速度响应的叠加，即

{\bar{v}}_{j}^{s} = {\tilde{\bar{v}}}_{j}^{s} + Y_{j k}^{s} {\bar{p}}_{k}^{c}

（7）

声场内分片j的平均声压响应 ${\bar{p}}^{a}$ 等于未耦合声场平均声压响应 ${\tilde{\bar{p}}}^{a}$ 与结构对声场耦合产生的平均声压响应之和，即

{\bar{p}}_{j}^{a} = {\tilde{\bar{p}}}_{j}^{a} + Z_{j k}^{a} {\bar{v}}_{k}^{c}

（8）

联立式（4~8），得到耦合面内任意分片上的平均速度响应与未耦合结构的平均速度响应关系，并用矩阵表示为

{\bar{v}}^{c} = (I - Y^{s} Z^{a})^{- 1} {\tilde{\bar{v}}}^{s}

（9）

式中：I代表单位矩阵，Y^s、Z ^a分别表示有限元离散后的结构导纳矩阵和声场阻抗矩阵。通过式（9）获得耦合面任意分片上的平均速度响应后，再由声场阻抗关系式（6）和平均速度连续性条件式（4）计算出声场任一点的平均声压响应。

为得到结构的PTF，由式（5）的定义，求解每个结构分片在单位均布压力作用下的平均速度响应。将第k个分片上的均布压力激励等效为节点压力，即

F_{k} = \{\begin{matrix} \iint_{S_{k}}^{} [N_{s}]^{T} {\bar{p}}_{k} d S & S_{k} \\ 0 & S_{c} / S_{k} \end{matrix}

（10）

式中S_c/S_k表示耦合面除第k个分片外的面域。节点位移响应由式（11）求得，进一步得到单元平均速度响应。

(K_{s} - ω^{2} M_{s}) U_{s} = F_{k}

（11）

同样地，将单个声场分片的单位平均速度激励等效为对应的输入节点速度，可以表示为^［

1］

Q_{k} = \{\begin{matrix} - j ρ_{0} ω \iint_{S_{k}}^{} {[N_{a}]}^{T} {\bar{v}}_{k} d S & S_{k} \\ 0 & S_{c} / S_{k} \end{matrix}

（12）

依据式（13）计算平均压强响应，从而得到声场的PTF

[K_{a} - ω^{2} M_{a}] P = Q_{k}

（13）

分片传递函数中的结构导纳和声阻抗的计算都是一系列分片内的平均值，因而符合声振测试要求^［

10］。在实际工程中，可以用测试值计算出结构导纳Y^s和声场阻抗Z^a，从而大大减少建模计算时间同时也能获得与测试结果一致的仿真结果。

3 仿真算例

为了验证PTFM的有效性，设计一个如图3（a）所示的四边简支板和空气声腔耦合模型，表1为模型参数。如图3（b）所示，在结构和声场内插入节点分别离散成平面四边形单元和六面体单元，共计N_s个结构节点和n_s个结构单元，N_a个声学节点和n_a个声学单元。

图3 声振耦合模型

Fig.3 Acoustic and vibration coupling model

表1 四边简支板‑空气声腔参数

Table 1 Parameters of four‑sided simply supported plate‑air acoustic cavity

参数	数值
声腔长L/m	0.304 8
声腔宽W/m	0.152 4
声腔高H/m	0.152 4
空气声速c₀/（m·s^-1）	340
空气密度ρ₀/（kg·m^-3）	1.21
平板厚度h/m	1.63×10^-3
平板密度ρ_s/（kg·m^-3）	2 690
平板弹性模量E/Pa	7×10¹⁰
平板泊松比ν	0.3

首先通过求解特征方程（1）计算该结构的低阶模态，前九阶固有频率及振型结果如图4所示。特别地，频率分别为217.839 3和540.350 8 Hz时，顶面中心处于共振状态，且该位置处响应在实际中也便于测量。基于此，为比较PTFM、文献［

4］方法与DCM得到的平均声压响应结果，选取声腔上顶面中心为参考点之一。另外选取体心为参考点反映有限元方法对声场的计算情况。

图4 声振耦合模型固有模态

Fig.4 Intrinsic modes of the acoustic and vibration coupling model

以1 Hz为步长，分别使用3种方法在频率为1~1 000 Hz范围内对结构参考点声压进行计算。图5，6（N_s=91，n_s=72，N_a=637，n_a=432）分别为声腔顶面中心参考点与体心参考点处声压随频率变化曲线。图5曲线峰值处的频率分别为218和540 Hz，对应图4耦合场的固有频率，且在该频率下，体心处的声压曲线同样到达峰值，符合第二阶、第四阶固有振型的对称性。综合两个参考点的声压响应曲线，PTFM的计算精度较文献［

4］方法更高。

图5 声腔顶面中心声压-频率曲线

Fig.5 Sound pressure level-frequency curves at the center of the top surface of the cavity

对相同结构设置不同的有限元划分如表2所示。为了说明在整个频带范围内的计算精度，以式（14）计算应用PTFM与DCM在每个频率点处平均声压响应的绝对误差（Absolute error， Ae），结果如图7所示。

表2 有限元节点和单元数目

Table 2 Number of finite element nodes and elements

单元划分	N_s	n_s	N_a	n_a
1	77	60	539	360
2	91	72	637	432
3	97	84	679	504
4	105	96	735	576

图6 声腔体心声压-频率曲线

Fig.6 Sound pressure level-frequency curves at the body center of the cavity

图7 绝对误差‑频率曲线

Fig.7 Ae‑frequency curves

A e = P_{P T F M} - P_{D C M}

（14）

式中P_（._）分别表示两种方法的声压响应值。在整个频带上，除曲线峰值和谷值处对应的频率点外，两种方法得到的频响曲线几乎吻合。

而在固有频率处，使用PTFM计算得到的声压响应绝对误差和按式（15）得到的相对误差大小（Relative error， Re）结果如表3所示。

R e = \frac{P_{P T F M} - P_{D C M}}{P_{D C M}} \times 100 %

（15）

表3 固有频率点计算误差

Table 3 Calculation error at natural frequency

固有频率/Hz	单元划分	绝对误差/dB		相对误差/%
固有频率/Hz	单元划分	顶面中心	体心	顶面中心	体心
218	1	0.89	0.80	0.78	0.69
	2	2.64	2.55	2.11	2.01
	3	-3.93	-4.00	-2.93	-2.93
	4	-1.84	-1.89	-1.41	-1.43
540	1	2.10	2.60	2.23	3.01
	2	-5.30	-4.93	-4.96	-5.00
	3	-1.76	-1.47	-1.81	-1.65
	4	-1.17	-0.88	-1.29	-1.07

依据计算结果，两参考点声压响应相对误差相近。根据式（9），耦合方程的振动响应是由未耦合场结构的响应计算得到，图8展示了部分未耦合场的固有频率和振型，未耦合场的固有频率始终与耦合场固有频率有所偏差，因而存在计算误差。当划分的单元更小时，刚度矩阵、质量矩阵的稀疏性和式（13）求逆运算的奇异性问题显著，综合影响计算误差。

图8 未耦合场的固有频率及固有模态

Fig.8 Natural frequency and natural mode of the uncoupled field

顶面中心与体心处的谷值频率分别为376、657和502 Hz，根据机械振动理论^［

11］，曲线谷值对应的频率值为该参考点的反共振频率。表4列出这些振频率点的声压响应计算误差。由式（11，13），分片上的传递关系定义为单元平均速度、平均声压响应的比，因此对某一点一定存在近似偏差。当提高有限元的数目时，分片单元面积减小，式（11，13）的平均效应被削弱，但仍存在上述稀疏性和奇异性问题影响计算精度。

表4 反共振频率点计算误差

Table 4 Calculation error at anti‑resonant frequency

单元划分	绝对误差/dB			相对误差/%
	顶面中心		体心	顶面中心		体心
	376 Hz	657 Hz	502 Hz	376 Hz	657 Hz	502 Hz
1	37.83	23.59	36.25	-1 254.37	351.41	-378.43
2	27.37	-10.66	19.03	382.54	-36.30	349.31
3	20.50	-2.81	8.66	146.84	-7.79	65.92
4	13.64	-8.44	2.72	58.84	-23.41	7.63

虽然在频响曲线谷处的响应有所偏差但总体的趋势是与直接法的趋势一致。其次，反共振频率是一种局部现象，即只对个别点的测量误差有影响，且由于计算步长为1 Hz，对于实际工程中更宽泛的频率间隔，精度也是可以保证的。

最后，在计算精度满足要求的情况下，使用PTFM计算时分片单元的类型和尺寸与结构单元相同，单元数目远小于声学单元数目，因此由式（11）计算的线性方程组维度较式（1）减少了N_a个，降低了有限元计算的复杂程度。关于PTFM在耦合方程的降维机理有待进一步研究。综上所述，表明了该方法的高效性。

4 结论

本文建立了一般结构声场的无阻尼声振耦合模型，引入了分片传递函数方法实现耦合方程的解耦。设计了一个四边简支板和空气声腔耦合模型，对比在声腔顶面中心和体心处文献［

4］方法、分片传递函数法和直接声振耦合有限元法的平均声压响应结果，分析并讨论了应用PTF原理在共振频率与反共振频率处的声压响应计算误差原因，验证了该方法的有效性。

参考文献

Sandberg G， Wernberg P A， Davidsson P. Fundamentals of fluid-structure interaction［M］. ［S.l.］：Springer， 2008. [百度学术]

Lee C M， Royster L H， Ciskowski R D. Formulation for an FE and BE coupled problem and its application to the earmuff-earcanal system［J］. Engineering Analysis with Boundary Elements， 1995， 16（4）： 305-315. [百度学术]

Desmet W. A wave based prediction technique for coupled vibro-acoustic analysis［D］. ［S.l.］：［s.n.］， 1998. [百度学术]

KIM S H， Lee J M. A practical method for noise reduction in a vehicle passenger compartment［J］. Journal of Vibration on Acoustics Transactions of the Asme， 1998， 120（1）： 199‑205. [百度学术]

Davidsson P. Structure-acoustic analysis， finite element modelling and reduction methods［D］. Lund， Sweden： Univerisity Libraries， 2004. [百度学术]

Xiao W， Li L， Lei S. Accurate modal superposition method for harmonic frequency response sensitivity of non-classically damped systems with lower-higher-modal truncation［J］. Mechanical Systems and Signal Processing， 2017， 85： 204-217. [百度学术]

Ouisse M， Maxit L， Cacciolati C， et al. Patch transfer functions as a tool to couple linear acoustic problems［J］. Journal of Vibration and Acoustics， 2005， 127（5）： 458‑466. [百度学术]

Maxit L， Aucejo M， Guyader J. Improving the patch transfer function approach for fluid-structure modelling in heavy fluid［J］. Journal of Vibration and Acoustics， 2012， 134（5）： 051011. [百度学术]

Aucejo M， Maxit L， Totaro N， et al. Convergence acceleration using the residual shape technique when solving structure-acoustic coupling with patch transfer functions method［J］. Computers & Structures， 2010， 88（11/12）： 728-736. [百度学术]

Mohamady S， Raja Ahmad R K， Montazeri A， et al. Modeling and eigen frequency analysis of sound-structure interaction in a rectangular enclosure with finite element method［J］. Advances in Acoustics & Vibration， 2009， 13： 1-9. [百度学术]

胡海岩. 机械振动基础［M］. 北京：北京航空航天大学出版社，2018. [百度学术]

Hu Haiyan. Mechanical vibrations fundamentals［M］. Beijing： Beihang University Press， 2018. [百度学术]

请使用 Firefox、Chrome、IE10、IE11、360极速模式、搜狗极速模式、QQ极速模式等浏览器，其他浏览器不建议使用!