张志飞，梁 涛，徐中明，夏小均，陶能发

前言

目前，声学的数值分析方法主要是基于单元技术的有限元法(FEM)、边界元法(BEM)、统计能量法(SEA)和它们的混合建模方法。由于有限元法和边界元法是基于单元的方法，因此随着分析频率的增加，为保证求解精度，网格尺寸必须减小，从而导致求解模型随着分析频率的增加而迅速增加，因此，有限元法大多是用于低频分析中[1-2]。而统计能量法由于其高模态密度的要求，只能在高频范围内使用，此外，统计能量法是一种统计方法，因此它是非确定性预测方法[3]。目前在中频段缺乏广泛应用的确定性中频预测方法。

WBM是以Trefftz理论为基础推导而来的间接Trefftz方法[4]。WBM不需要对结构和声域划分网格[5]，而是采用严格满足控制微分方程的波函数的线性叠加来描述系统的动态响应变量[6-7]，该方法在求解区域内不存在近似误差，但是有可能不满足边界条件和连续性条件，为此该方法采用加权余量公式使其在边界和连续界面处的加权余量等于零，这样就能构造一个矩阵方程组，通过计算即可得到基波函数的加权系数，从而求得系统的响应[8]，因此，WBM具有模型规模小的特点。此外，WBM还具有高收敛性，与FEM相比其计算效率更高，故可以应用于更高频率范围的声学分析预测[9-11]。但是，WBM收敛的一个充分非必要条件是所求声域为凸形域，因此其几何适应能力远不如FEM。

针对车内声腔这类边界几何复杂的非凸形域内部声学问题，无法满足WBM求解域为凸形域的收敛条件。若采用FEM，为保证求解精度，随频率增加网格尺寸要求越精细，导致模型规模迅速增大，计算的时间成本急剧增加，从而限制了其在更高频率范围的应用。针对上述两种方法的特点，为高效的实现此类复杂几何域的中低频声学响应，本文中以声压和法向速度为连续性条件，采用了兼具两种数值方法优点的混合FE-WBM方法，并完成其理论推导。以某二维车内声域为例，运用混合方法建立了数值求解模型。通过与传统有限元方法的对比，验证了混合FE-WBM的正确性与高效性。

1　声学问题

对内部有源的凸形封闭声学空间Ω，任意位置r的稳态声压p满足Helmholtz控制方程[12]:

对于内部声学问题，如图1所示，边界Γ由3部分构成(Γ＝Γp∪Γv∪Γz)，在Γp边界施加声压p，在Γv边界施加法向速度vn，在Γz边界施加法向阻抗Z[13]，即

1.1　声学有限元法

图2是有界声腔Ω的网格，每个单元Ωe的边界由4部分组成，即Γe＝Γep∪Γev∪Γez∪Γei。Γep，Γev和Γez是声学边界与单元边界的交集(Γep＝Γe∩Γp，＝Γe∩Γv，Γez＝Γe∩Γz)，Γei表示相邻单元的公共边界。

图2　有限元划分与单元定义

在每个单元Ωe内部，采用简单基函数的线性组合p＾(r)来近似代替其精确解p(r):

式中:贡献因子pea组成的列向量pe是未知的自由度，即节点声压。

在相邻单元的公共界面施加速度连续性条件，得到的单元模型为

式中:Ze＝ - ω2Me＋ jωCe＋ Ke；fe＝jω(qe-ven-是成对出现，在组装整体FE模型中相互抵消。将所有的单元矩阵组合便得到整个模型的FE模型:

式中:Z＝ -ω2M ＋jωC ＋K；f＝ jω(q-vn)；M，C和K分别表示声学质量、阻尼和刚度矩阵；q和vn分别表示与点生源q和法向速度vn相关的载荷向量。

1.2　声学波函数法

WBM是采用全局Trefftz基函数Ψa和声压特解函数的线性组合p＾来近似声压精确解p:

式中:行向量Φ由Trefftz基函数Ψa组成，列向量p是由未知加权系数pa组成，na表示波函数的个数。本文中选取的声学基函数[14]为

式中kxa，kya和kza满足如下关系:

声腔中rq处放置一声源时，在r处的声压特解函数为

理论上，当波数趋于无穷大时，式(6)的解收敛于精确解，但在数值计算中仅能取有限数目的基波函数，因此需要采用截断法则。

基于SPIEGEL针对半幅傅里叶级数提出的理论，即余弦函数的最小周期不应大于其近似函数周期的一半，即截断系数T≥2。最小空间周期λi，min和声波波长λ定义为

式中Li(i＝x，y，z)表示包含求解域最小立方体在x，y，z方向的物理尺寸。

故截断法则Tλi，min≤λ可表示为

近似声压p＾(r)满足Helmholtz方程，因此，通过加权积分使压力近似值p＾(r)满足式(2)中的边界条件，从而得到波模型，求解后便可得到未知波函数加权系数p。加权余量公式为

其中:

2　混合有限元-波函数法

为实现复杂系统在中频的高效、高精度声学预测，采用了充分结合FEM和WBM各自优势的混合FE-WBM。FE-WBM的基本思想是将WBM用于大型均匀的简单几何体，充分利用WBM数值模型小、计算效率高的特点，同时将FEM用于几何复杂的区域，充分利用FEM较强的几何适应能力，最终创建的混合FE-WBM模型的自由度将比纯FEM模型大幅减少，在提高模型计算效率的同时，保证了模型的良好几何适应性。在FE(finite element)域和WB(wave based)区域的耦合界面处Γi，采用压力连续和法向速度连续的耦合条件[15]，即将压力连续性条件以边界条件的形式施加到WBM区域，同时将法向速度连续性条件施加到耦合界面的FEM侧，如图3所示，即

图3　FE-WB直接耦合方式

混合模型的FEM部分为

令 Zfe＝-ω2M ＋jωC ＋K，ffe＝jω(q-vnfe)，则上式可表示为

式中:Qfw为FEM对WBM的耦合矩阵；pw为波函数的加权系数矩阵。

混合模型的WBM部分为

其中:

式中:Qwf为WBM对FEM的耦合矩阵；Cw和cw分别为WBM的反馈矩阵和反馈向量。

由式(18)和式(21)得到FE-WBM耦合模型:

求解式(23)便可得到有限元部分的节点压力pfe和 pw。

式(23)的系数矩阵不再是带状稀疏矩阵，因此若采用直接法求解方程组，求解效率低且未充分利用有限元矩阵Zfe的带状稀疏特性。因此，将式(23)的求解分为以下步骤[16]。

(1)将pw视为已知，求解节点声压矩阵pfe:

将式(24)带入式(21)，整理得

通过求解如下两个线性方程得到矩阵H和h:

因为Zfe是带状稀疏矩阵，因此可以采用高度优化的稀疏带状矩阵求解算法。

(2)求解式(25)，得到 pw。

(3)节点声压矩阵pfe为

3　数值算例

为验证上述方法的可行性，建立了一个简化后的二维车内声腔模型，如图4所示。由图可知，车内声腔的边界几何较复杂，尤其是在声腔与座椅的交界处。如前所述若想采用WBM，则必须将其求解域划分为互不重叠的凸形域，然而该声腔与座椅交界处的几何形状复杂，且座椅头枕等区域无论如何划分都无法保证其WBM子域为凸形域，因此仅采用WBM法是不可行的，但可以采用混合FE-WBM数值方法求解该车内声腔的声压响应。将几何较复杂的座椅区域附近采用FEM，几何较简单的其余区域采用WBM，如图5所示。图中填充区域表示采用有限元法的区域，未填充区域表示采用波函数法的区域，虚线表示各个区域之间的公共边界，将整个车内声腔划分为4个WB域和一个FE域。速度边界施加在WB域，声压响应点在FE域。在左侧边界上施加法向速度v-n＝10m/s，其余边为刚性壁边界，即法向速度为零。介质为空气，其材料参数:密度ρ＝1.225kg/m3，声速c＝340m/s。声压响应点取为驾驶员耳朵附近。

图4　车内声腔模型

图5　车内声腔区域划分

3.1　结果分析

FE-WBM混合模型中，WBM的截断参数T＝4.0，各WB域之间在耦合边界处采用速度-声压连续性条件。为说明FE-WBM的有效性，建立了由4 936个线性四边形单元组成的粗有限元模型和由12 669个线性四边形单元组成的精细有限元模型，并将精细有限元模型作为参考模型。各模型在600和900Hz的声压分布云图如图6和图7所示。图6(b)和图7(b)中的虚线是声腔区域划分线，车内声腔声压分布云图连续，在各WBM子域与FEM域的公共边界处压力分布连续，未出现明显的畸变。从600和900Hz各模型的声压云图对比可知，FE-WBM模型和参考模型的计算结果，无论是声压分布还是声压数值上均符合得较好，而粗有限元模型与参考模型的声压分布存在一些差异，在900Hz时差异较600Hz时的更加明显。

图6 600Hz时声压分布云图

图8(a)和图8(b)分别是二维车内声腔中响应点在1～1 000和500～1 000Hz的声压频率响应曲线。由图8(a)可知，在1～1 000Hz内，三者的整体趋势一致，FE-WBM与参考模型的计算结果，除在个别频率附近存在一些频率和幅值的偏差外，在整个频率范围内均符合得较好。而粗有限元模型在稍高频率范围内相对参考模型出现明显的频率偏移，这一点在图8(b)中更加直观。通过声压分布云图和声压响应点频率响应曲线说明FE-WBM的计算结果精度较高。

图7　900Hz时声压分布云图

图8　二维车内声腔响应点声压频响曲线

3.2　收敛性分析

为比较各方法的计算精度和效率，分析了290Hz时混合FE-WBM与FEM模型的声压响应收敛性。相对声压误差[14]定义为

式中:p＾表示混合FE-WBM或FEM计算的响应点声压；pref表示参考声压，即精细有限元模型计算的响应点声压。

混合FE-WBM模型和FEM模型均是在i5-6500CPU(3.20GHz)，8GB运行内存，64位操作系统，Mat lab R2016a平台运行。图9和图10分别为290Hz时相对声压误差随模型自由度和CPU运行时间增加的变化曲线。虽然WBM的收敛性随着域分解的增加会有所下降[4]，但从图9和图10可知，达到相同的计算精度，FE-WBM的模型自由度和CPU运行时间较FEM依然存在明显的优势。

可见，采用FE-WBM时，将FEM应用于几何复杂的区域，WBM用于简单几何的区域，在保证其良好的几何适应能力的同时也保证了计算精度和效率，因此混合FE-WBM可用于更高频率的分析。

图9　相对误差-自由度收敛曲线

图10　相对误差-CPU运行时间收敛曲线

4　结论

本文中采用结合FEM的强几何适应能力和WBM的高收敛率的FE-WBM，通过声压和法向速度连续性条件实现二者的耦合。以简化后的二维车内声腔模型为例，将几何复杂的区域采用FEM，将几何简单的区域采用WBM，建立了二维声腔的混合FEWBM模型。通过对比粗有限元、精细有限元和混合FE-WBM的声压云图和响应点声压频率响应曲线，证明了FE-WBM的有效性。通过收敛性分析，说明了FE-WBM相对FEM的高收敛率，且FE-WBM中有限元部分的网格越精细，其计算精度越高，与此同时计算时间虽有所增加，但相对FEM依然存在明显的优势。今后的工作是将混合FE-WBM用于三维声腔和结构-声学耦合系统分析中。

[1] 高印寒，张澧桐，梁杰，等.基于声-固耦合模型的车内噪声响度参数的预测与分析研究[J].机械工程学报，2016，52(12):136-143.

[2] 徐中明，何治桥，贺岩松，等.车内低频噪声预测模型的改进[J].汽车工程， 2016，38(7):878-882.

[3] MACE B，DESMET W，PLUYMERSB.Mid-frequency methods in sound and vibration-part A[J].Journal of Sound&Vibration，2013，332(8):1895-1896.

[4] WIM D.A wave-based prediction technique for coupled vibro-acoustic analysis[D].Leuven:University of Leuven，1998.

[5] DECKERSE，ATAK O，COOX L，et al.The wave based method:an overview of 15 years of research[J].Wave Motion， 2014，51(4):550-565.

[6] KITA E，KAMIYA N.Trefftz method:an overview[J].Advances in Engineering Software， 1995，24(1-3):3-12.

[7] BERGEN B，VAN GENECHTEN B，VANDEPITTE D，et al.An efficient trefftz-based method for three-dimensional helmholtz problem in unbounded domains[J].CMES， 2010，61(2):155-175.

[8] PLUYMERS B，VAN HAL B，VANDEPITTE D，et al.Trefftzbased methods for time-harmonic acoustics[J].Archives of Computational Methods in Engineering， 2007，14(4):343-381.

[9] DESMET W，HAL B V，SASP，et al.A computationally efficient prediction technique for the steady-state dynamic analysis of coupled vibro-acoustic systems[J].Advances in Engineering Software，2002，33(7):527-540.

[10] PLUYMERSB，DESMET W，VANDEPITTE D，et al.Application of an efficient wave-based prediction technique for the analysis of vibro-acoustic radiation problems[J].Journal of Computational and Applied Mathematics， 2004，168(1-2):353-364.

[11] VERGOTE K，VANMAELE C，VANDEPITTE D，et al.On the use of an efficient wave based method for steady-state structural dynamic analysis[C].Proceedings of LSAME.08.Leuven，Belgium:KU Leuven Press，2008:433-459.

[12] 何雪松.结构—声学耦合分析及其波函数法研究[D].武汉:华中科技大学，2008.

[13] MARESSA A，JONCKHEERE S，VAN GENECHTEN B，et al.Hybrid finite element-wave based techniques for interior vibro-acoustics.application to a large-sized vehicle model[C].Proceedings of ISMA2012-USD2012.2012，8(12):3927-3940.

[14] VAN HAL B，DESMET W，VANDEPITTE D，et al.Hybrid finite element-wave based method for acoustic problems[J].Computer Assisted Mechanics and Engineering Sciences， 2003，10(4):375-390.

[15] PLUYMERS B.Wave based modelling methods for steady-state vibro-acoustics[M].2006D04.Leuven:KU Leuven，2006.

[16] PLUYMERS B，DESMET W，VANDEPITTE D，et al.Wave based modelling methods for steady-state interior acoustics:an overview[C].Proceedings of ISMA2006.Leuven，Belgium:KU Leuven Press，2006:2303-2358.