二维多频声辐射问题的频域边界元分析方法

钦 宇，周枫林,，王炜佳，张玉良，袁小涵

（1.湖南工业大学 机械工程学院，湖南 株洲 412007；2.衢州学院 机械工程学院，浙江 衢州 324000）

1 背景知识

离心风机出口气动噪声、机械设备运转时摩擦噪声等声源的主要频率成分在频谱上具有一定宽度，而多频声源的辐射问题一直是工程师关注的重要问题。适用于声辐射分析的工程数值方法以边界元法（boundary element method，BEM）[1-3]、有限元法（finite element analysis，FEM）[4]和矩量法（method of moments，MoM）[5]为主要代表，其中BEM 和MoM 都以边界积分方程为理论基础，特别适用于无限域的声辐射分析。

边界元法是一种基于边界积分方程为数学基础的数值分析方法，其边界积分方程能够自动满足无穷远场的辐射条件，并且仅需要对边界进行离散。较之其他的数值计算方法，边界元法引入了基本解，使其具有了解析和数值相结合的特点，从而使其计算精度相对较高，是一种比较适用于分析声学问题的计算方法[6]，目前被广泛应用于求解裂纹问题[7]、扩散波问题[8]、热传导问题[9]等。

目前对声学问题的分析主要采用频域法[10]和时域法[11]，时域类方法具有直接且通用性好的优点，然而在选择时间步长时会遇到数值不稳定问题。频域法通过Helmholtz 控制方程和基本解，推导出边界积分方程，再将边界积分方程进行线性离散，建立矩阵，求解线性方程组来获得相应的解[12]。频域法在频率成分预判较准确的基础上能够获得较高的计算精度，通过反变换获得时域解。

目前，数值反变换的方法主要有 Laplace 数值反变换法[13-14]，Stehfest 算法[15]，离散傅里叶反变换法[16-18]，工程中常利用 Laplace 数值变换结合边界元法计算结构动力响应问题，首先需要得到变换域中的一组所对应的位移和应力，然后通过 Laplace 数值反变换获得时域中相应量的解[19-21]。然而对于多频噪声的分析，需要使用多个采样频率并借助数值反变换技术才能获得满足要求的时域解答。

本研究中先应用频域边界元法分析二维多频声场辐射问题，将时域描述的控制方程通过 Fourier变换转换为频域描述的控制方程（Helmholtz 方程），再求解等间隔采样频率下的 Helmholtz 方程， 然后采用离散傅里叶反变换（inverse discrete Fouriertransform，IDFT）方法获得时域解。最后以两个不同区域的内声场问题算例，通过与参考声压进行比较，论证了所提方法的正确性。

2 二维声学频域边界元法

在均匀理想介质E中，二维声学问题在空间内传播的波动方程表示形式为

x为配置点空间坐标；

c为介质中声音传播的速度。

对式（1）进行Fourier 变换，假设时间项选取为e-jωt，可将瞬态问题从时间域转化到频率域进行简单稳态分析，则频域声压可以表示为

j 为虚数单位，j2=-1；ω为声源简谐振动角频率。

将式（2）代入（1），消去时间项，即可得如下Helmholtz 控制方程。

式中k为波数，且k=ω/c。

控制方程的边界条件一般满足以下3 类。

Dirichlet 边界条件，即声压已知：

式中符号¯表示已知值。

Neumann 边界条件，即法向速度已知：

式中：q(x)为声压在x处的法向导数；n为单位外法向，指向背离声域；ρ为声学介质的质量密度；vn为法向速度。

Robin 边界条件，即阻抗边界条件已知：

式中Z为声阻抗。

将方程（3）转化为等效积分形式，并将试函数取为二维声学问题的基本解：

式中：c(x)为自由项；

y为场点空间坐标；

q(y)为声压在y处的法向导数；

Sy为结构边界S的子边界；

Gk(x,y)和为二维和三维声学问题的基本解，且

r为源点与场点之间的距离，且r=|x-y|；

x为配置点（即计算点或源点）；

y为场点。

为了得到式（7）的数值解，需要根据不同的要求与情况，将边界划分为一些小单元。本文主要采用线性离散方法，将边界S划分为Ne个线性单元，则边界积分方程可以离散为

式中：u为任一边界的节点，u=1, 2, …,Ne；

Sv为边界上的单元v；

qv为在单元v上的法向导数；

将源点置于所有场点上可得到：

式中Huv和Guv均为影响系数矩阵，为遍及所有节点单元上的积分总和。

也可将式（1）写为

可以将上述方程组表示为如下矩阵形式：

此时，利用边界条件，将已知量移到方程右边，未知量移至左边可得：

式中：A为未知声压值和位置法向导数组成的Ne×Ne系数矩阵；

b为已知量组成的Ne阶已知向量；

h为由包含未知声压值和未知法向导数组成的Ne阶向量。

3 离散傅里叶逆变换

为实现离散傅里叶反变换运算，需将连续信号进行截断、采样。

在频域边界元法中，首先取N个等间隔不同的采样频率，分别求解线性方程组（14），得到一组解x(ω)，再通过离散傅里叶逆变换（IDFT）得到时域解x(t)。设x(t)是周期为T的离散函数，每个周期内有N个等间隔频率离散点。其中Δω=ω0=2π/T。为实现IDFT 运算，必须通过频率采样使频域函数为有限离散值，此处采用δ函数：

式中：Δt为时间步长，其中Δt=T/N；

ω0为初始频率点。

为避免产生频率混叠现象，信号采样频率

同时为避免频混现象，可通过选择较小的采样间隔T（即较高的采样频率）来减小这种误差。

为保证离散后的信号能唯一确定原连续信号，采样一般等间隔取值，其频率采样信号为

其Fourier 变换：

根据δ函数的筛选性，化简求得其Fourier 系数Ck为

将式（20）代入式（19），可求得：

一个周期内N个采样点的复数值为

式中N为采样频率点的数量。

由T=2π/ω0=NΔt可得：

傅里叶级数可以写成复指数形式，根据欧拉公式

可得到：

4 数值算例

为验证二维频域边界元法结合离散傅里叶反变换求解多频声辐射问题的准确性，设计了两个验证算例，将计算结果与解析解进行对比以验证该方法的正确性。

4.1 旋转L 形边界区域内部声压变化

在一个旋转L 形边界区域条件下进行内点声压变化情况分析，其边界节点之间间距设为0.1 cm，x、y分别为内部点的横坐标与纵坐标（算例中都将如此使用），计算求解过程中一共采用了64 个边界节点，旋转L 形边界区域节点分布情况如图1所示。

图1 旋转L 形边界节点分布图Fig.1 Distribution diagram of rotating L-shaped boundary nodes

图1所示边界上的声压分布按照如下公式给出：

在y=0 的边界上，

在x=0.8 的边界上，

在y=0.8 的边界上，

在x=1.6 的边界上，

在y=1.6 的边界上，

在x=0 的边界上，

域内声压精确解为

计算采用在频率区间1~4 内等间隔取7 个频率点，计算时间区间为0.5~3.0 s，在域内选取两个观察点，分别为Q1(0.4, 1.2)与Q2(0.5, 1.2)，计算两个内部观察点声压值随时间变化的过程，Q1、Q2声压值随时间变化的情况，如图2所示。

图2 Q1、Q2 声压值随时间变化情况Fig.2 Variation of Q1 and Q2 sound pressure value with time change

从图2中可以看到，在该方法下两个取样点上的声压计算值与精确值吻合较好。为更准确直观地展示计算结果和精确解的差异，于表1~2 中列出了随机选取的部分时间点上，对于Q1和Q2两个内部观察点的声压以及通量计算结果与精确解数值。

从图2以及表1与表2中可以得出，在该分析方法下旋转L 形边界区域域内点的计算结果和精确解不仅在变化趋势上基本保持一致，而且其具体计算数值与精确解下的精确值较为吻合，误差较小，充分说明该方法的准确度较高。

表1 内部取样点Q1 和Q2 上声压计算值与精确值Table 1 Calculated and exact values of sound pressure at internal sampling points Q1 and Q2

表2 内部取样点Q1 和Q2 上通量计算值与精确值Table 2 Calculated and exact values of flux at internal sampling points Q1 and Q2

4.2 凸形边界区域内部声压变化

在一个凸形边界区域条件下进行内点声压变化情况的分析，其边界节点之间间距同样设为0.1 cm，区域声速设为1 个标准单位，计算求解过程中一共采用了80 个边界节点，凸形边界区域的节点分布情况如图3所示。

图3 凸形边界节点分布图Fig.3 Distribution diagram of convex boundary nodes

图3所示凸形边界上声压分布按如下公式给出：

在y=0 的边界上，

在x=2.4 的边界上，

在y=0.8 的边界上，

在x=1.6 的边界上，

在y=1.6 的边界上，

在x=0.8 的边界上，

在x=0 的边界上，

此域内声压精确解为

计算时间区域为 0.5~4.0 s，取两个域内点Q3(1.5,0.4)、Q4(1.6, 0.4) 作为该边界条件下的观察点，在频率区间 1~4 内等间隔取 7 个频率点计算。分析Q3及Q4声压值随时间变化的过程，比较结果如图 4所示。表3中列出Q3和Q4的部分计算结果与精确解的对比。

由图4和表3可知，在凸形边界条件下选取的两个域内取样点上声压计算值与精确值吻合情况也是较为理想的。

表3 内部取样点Q3 和Q4 上声压计算值与精确值Table 3 Calculated and exact values of sound pressure at internal sampling points Q3 and Q4

图4 内部观察点Q3、Q4 声压值随时间变化情况Fig.4 Variation of sound pressure value at internal observation point Q3 and Q4 with time change

比较旋转L 形与凸形算例的结果，发现在不同边界条件下其域内点的计算结果都能够和精确解基本保持一致，而且在合理的采样区间内能达到较高的计算精度，这充分说明了所提方法的可行性以及其具有较高的时域解计算准确度。

5 结语

针对含有复杂频率成分声源的辐射问题，采用离散傅里叶反变换和频域边界元法相结合，研究了一种二维多频声辐射问题的频域边界元分析方法。采用傅里叶变换将时域描述的波动方程转化为Helmholtz 方程，并选取多个等间隔频率点作为采样频率，应用边界单元法求解各特征的Helmholtz 方程，获得不同位置在各采样频率下的声压，采用离散傅里叶反变换将频域声压幅值和相位转化成时域声压。两个数值算例分析下的域内声压计算值和精确解吻合较好，误差较小，证明该方法具有可行性，并具有较高的时域解计算准确度。