不同二维声学边界单元在敏感度分析上的数值表现

(信阳师范学院 建筑与土木工程学院，河南 信阳 464000)

优化方法可分为直接方法和近似方法[1]，直接方法的代表有遗传算法和退火算法[2-3]，这类方法是通过合适的参数变化寻找全局最优值，但该类方法往往需要大量的函数计算。作为对比，近似方法只需少量的函数计算就可以达到收敛，但该方法往往只收敛到局部最小值。大多数局部近似方法都是基于目标函数在某一设计点附件的泰勒展开进行寻优计算的，因此在每个迭代步都需要进行敏感度计算，这一步往往是十分耗时的。有限差分法，因其很容易的被实施，已被广泛应用于敏感度计算，然而该方法计算精度不高，并对每一个设计变量都需要进行差分计算，计算量大。解析与半解析方法可以有效用于敏感度计算，例如直接微分法[4-6]。该方法将目标函数对设计变量直接进行微分计算，计算精度高，然而对于多设计变量的优化分析，在每一个迭代步对每一个设计变量都要进行微分计算，这会大大降低计算效率。伴随变量法[7]，通过引入与设计变量无关的伴随方程，在进行敏感度分析时，该伴随方程只需求解一次，然后间接计算目标函数的敏感度值，因此有效提高了计算效率。本文推导出基于伴随变量法的声学边界元敏感度表达式，然后采用不同类型边界单元进行声学敏感度计算，并对比其计算效率。

实际上在早期的积分方程应用中，包括最近的许多研究[8-10]，常量单元得到了广泛的应用，因其实施起来简单。但其计算精度差，在实际工程计算中，为了达到一定的精度要求，往往要加大离散，增加了计算量。很多学者使用连续线性和二次单元进行数值计算，以提高计算精度。然而连续单元对网格的适应性不高，在角点处需要特殊处理，难以进行复杂边界的积分计算，然而高阶非连续单元可以很好的克服这些问题。文献[8]和文献[9]分别给出了二维和三维非连续高阶单元的计算过程。非连续单元的插值节点位于单元内部，插值型函数的表达依赖节点位置。文献[10]给出了不同类型连续与非连续边界单元的计算精度对比，并研究了单元误差与频率、网格尺寸和插值节点位置的关系，并得到在相似自由度下非连续元比连续元表现更有效的结论。上述文献都是基于传统边界元法进行不同类型单元的计算精度比较的，使用传统边界元法进行外声场分析时，在某些虚假本征频率处可能得不到精确结果。另外，这些文献都没有进行敏感度分析。本文使用Burtion-Miller[11]法克服虚假本征频率问题，首先推导出基于不同类型边界单元离散下的无奇异Burtion-Miller边界积分表达式，然后使用带有解析解的简单算例对比不同类型边界单元在进行声学敏感度分析时的计算表现，并研究在进行敏感度分析时的数值误差随非连续单元内部插值节点位置变化关系，以考察得到优化节点位置参数值。

最后本文使用MMA方法进行典型声屏障结构优化分析，以验证本文发展算法的有效性。

1 二维声学边界元

1.1 无奇异边界积分方程

流体中声场分布满足Helmholtz控制方程，通过格林公式变换后，可得到如下的边界积分方程：

(1)

(2)

式中：

(3)

(4)

(5)

(6)

如果在点处边界光滑，系数为1/2。单独使用式(5)或(6)求解外声场问题，在虚拟频率处会产生非唯一解问题，通过组合这两个方程构成Burton-Miller表达式可以有效克服这一问题[11-12]。然而式(5)和(6)含有奇异积分项，为了得到精确结果必须对此做特殊处理，本文使用Cauchy主值积分和Hadamard有限积分法直接精确的计算奇异积分。接下来给出不同类型单元离散下的，无奇异积分表达式：

(7)

(8)

(9)

式中m表示每个单元的插值节点数目，Φ表示插值函数。函数f(x)可以为声压 φ(x)或声压通量q(x)等。

接下来首先给出不同类型边界单元的插值形函数表达式：

(1)常量单元

对于常量单元来说，满足m=1,Φ1=1。

(2)线性单元

对于线性单元来说，插值形函数Φ表达为：

(10)

式中ζ表示单元内部局部坐标，β表示非连续单元内部插值节点的局部坐标。当β=1时，式(10)为连续线性单元的插值函数表达式。

(3)二次单元

对于二次单元来说m=3，插值函数Φ如下表达：

(11)

当β=1时，上式为连续二次单元的插值形函数。

使用不同单元离散边界时，可得到方程(7)与(8)中的不同b(xj)和D(xj)的表达式，分别如下：

(1) 常单元

(12)

式中，系数C1和C2分别表达如下：

(13)

(14)

(2)线性单元：

(15)

式中， xl表示边界单元中不同于节点xj的另外一节点。α表示节点xj的局部坐标，当β≠1 时，β=|α|，上式中的系数可表达如下：

(16)

(17)

(18)

(19)

式中：

(20)

(21)

(22)

当β=1时，非连续线性单元退化为连续线性单元，式(15)中的系数表达式变化为：

(23)

(24)

(25)

(26)

(3)3二次单元：

(27)

式中，xl和xm表示边界单元中除了源点xj外另两个插值节点。当β≠1 时，式(27)中系数表达式为：

(28)

(29)

式中，

(31)

在此，我们推出了基于不同类型边界单元离散下的无奇异边界积分表达式，离散并组合(7)式与(8)式后，可得到如下的线性系统方程：

(32)

将未知量转移到方程左边，已知量转移到方程右边，式(32)可重新表达如下：

(33)

1.2 基于伴随变量法的声学敏感度分析

声学设计敏感度分析能够提供结构形状变化和结构声学表现之间的关系，因此敏感度分析时结构声学优化设计中极其重要的一个环节。任意目标函数W可以定义为：

(34)

(35)

同时边界和域的微分表达式可写为：

(36)

(37)

(38)

(39)

(40)

伴随方程可以定义为：

(42)

根据下面边界条件：

(43)

(44)

将式(42)、(43)和(44)代入式(41)中，可以得到下面表达式：

同样的，

(46)

将式(46)代入式(45)，可以得到：

可以看出式(47)不含有边界上未知量的敏感度。边界声压的梯度和伴随函数 被引入到目标函数的敏感度方程中。

3 算例分析

3.1 单元类型和误差定义

非连续单元的插值节点位于单元内部，插值节点位置参数值决定了插值形函数的表达，并决定了数值计算精度，因此研究节点位置对计算精度的影响以得到优化节点位置参数值是十分重要的。在图1中'DBEmn'和'CBEmn'分别表示带有'm'个几何插值点和'n'个物理插值点的非连续和连续单元类型。例如，图中'DBE21'表示线性几何近似的常量边界单元，'DBE22'表示非连续线性边界单元，'CBE22'表示连续线性边界单元，而'CBE33'表示连续二次边界单元，'DBE33'表示非连续二次边界单元。对于三维不同连续和非连续单元类型的比较可以参照文献[10]。常单元的物理插值节点位于单元中心处，而非连续单元的物理插值节点位置由参数 (0< <1)决定，见图1。另外，本文使用基于复数形式的声压面误差来表征数值计算精度，对于面误差的详细定义，请参考文献[10]。

图1 不同类型二维边界单元的几何与物理插值点分布Fig.1 Distribution of geometrical nodes and interpolation nodes

3.2 无限长圆柱壳散射算例

本小节进行无限长圆柱壳在平面波入射时的散射声场分析，平面波入射方向垂直于圆柱壳轴线，该模型可以简化成二维问题，平面波幅值为1Pa，圆柱壳半径 m。文献给出了该模型的声场解析表达，如下：

(48)

式中 代表纽曼符号， 代表对 进行求导。选取圆柱半径 为设计变量，将上式对设计变量进行微分，可得到如下声压敏感度表达式：

(49)

图2给出相似自由度下不同类型边界单元计算精度的对比，用于计算面误差的曲线选取为距圆心 处的圆环上，计算频率为5000Hz。观察该图可以发现，面误差随着计算节点数增加而降低，连续线性单元CBE22精度最差，非连续二次单元DBE33精度最高，非连续单元比连续单元表现更好。

图2 不同类型单元计算精度对比Fig.2 Comparison of numerical performance for different types of boundary elements

图3给出了分别使用单Helmholtz边界积分方程(5)式(亦称作CBIE方程)和Burton-Miller方程进行数值计算得到的结果与解析解的对比，DBE33单元用于数值离散，离散单元数为10000.观察该图可以发现使用CBIE方法在某些虚假频率处计算结果偏离解析解，然后使用Burton-Miller方法得到的结果与解析解符合一致。

对于非连续单元，选取合适的节点位置参数值以提高单元在计算精度方面的表现是十分关键的。图4呈现非连续线性边界单元DBE22插值节点位置参数值对计算精度的影响，依次选取4个计算频率进行数值误差分析，离散单元数为10000。观察该图可以发现，随着计算频率的增加面误差相应增加；对于每条频率误差曲线面误差都是先随着节点参数值增加缓慢降低，然后迅速降低，在达到最低值时迅速增加，最后又开始缓慢增加。优化节点位置参数值一致逼近0.57，该值非常接近勒让德多项式的零值。图5给出DBE33单元插值节点位置参数值对计算精度的影响，优化节点位置参数值一致逼近0.77，该值非常接近勒让德多项式的零值。

图3 计算点 处声压敏感度实部与虚部分别随频率变化关系Fig.3 Acoustic pressure sensitivity at point with different frequencies

图4 DBE22单元面误差随节点位置参数变化关系Fig.4 Surface error in terms of nodal position for DBE22 element

图5 DBE33单元面误差随节点位置参数变化关系Fig.5 Error in terms of nodal position for DBE33 element

图6 不同单元声压敏感度误差对比Fig.6 Surface error of acoustic pressure sensitivity for different types of elements

图6给出了相似自由度下不同类型单元数值精度对比。观察该图可以发现，同等阶次的非连续单元比连续单元表现的更有效，比如DBE33比CBE33精度高，DBE22比CBE22精度高。另外，通过对比DBE33与DBE22可以发现，二次单元比线性单元表现更好。因此，同等自由度下二次非连续单元表现最有效。

3.3 多个无限长圆柱壳散射算例

在这一小节里，我们进行多级散射问题分析。首先考虑四个无限长圆柱体在平面波作用下的声散射问题，圆柱半径为0.2m，圆心位置和坐标原点位置如图7所示。

图7 DBE22单元面误差随节点位置参数变化关系Fig.7 Error in terms of nodal position for DBE22 element

360个计算点均匀分布在以原点为圆心，半径为2m的圆环上。声压幅值敏感度值如下表达：

(50)

图8给出在这些计算点处的声压幅值敏感度值分布，设计变量为圆柱壳的半径，DBE33单元被用于数值计算，波数k=4。观察该图可以发现，使用快速多极边界元法得到的计算结果与传统边界元法得到的结果符合一致，表明快速算法的使用保持了传统边界元法的高计算精度特点。

图8 半径r=2m圆环上计算点声压幅值敏感度分布Fig.8 Sensitivity for the amplitude of sound pressure at the sample points distributed on a circle with radius r=2m

图9和图10分别给出4个和400个圆柱壳在平面波入射时的散射声压敏感度分布 ，DBE33单元用于数值计算，波数为k=4。图10中的圆柱壳分布在1m×1m的方格子里，每个圆环被离散成1000个单元，总共形成4×105个单元，快速多极算法被用于加速边界元系数矩阵与向量相乘的计算和系统方程的求解。这两幅图显示出在进行大规模多域问题的分析时，本文发展的二维快速敏感度算法具有较好的应用潜力。

图9 4圆柱散射时的声压幅值敏感度分布Fig.9 Sensitivity field for the amplitude of sound pressure from 4 rigid cylinders scattering

图10 400圆柱散射时的声压幅值敏感度分布Fig.10 Sensitivity field for the amplitude of sound pressure from 400 cylinders scattering

3.4 Y型声屏障结构优化分析

随着我国汽车保有量的迅速增加，引发的交通噪声污染问题也日益严重。因此，降低降低噪声是一项越来越引起重视的研究课题。在噪声源和居民区之间放置声屏障可以有效的降低保护区内的声压级[13]。虽然声屏障的高度对其降噪效果影响很大，然而考虑到车内人员的视觉感受和建筑成本，声屏障的高度是受到限制的。在高度一定情况下，声屏障顶端结构对其降噪表现也有很大的影响，常用的顶端结构是Y型。本文对典型的Y型声屏障进行二维优化分析，以得到降噪效果更好的顶端结构外形。如图11所示，声源距离地面高1.0m，离声屏障10.4m，Y型声屏障顶端左右两边角度分别为θ1和θ2，长度分别为l2和l1，低端高为3.0m，分布在声隐区的6个计算点坐标分别为(15,2)、(30,2)、(45,2)、(15,5)、(30,5)和(45,5)。

图11 Y型声屏障截面图Fig.11 Cross section of the Y-shaped noise barrier

本文采用MMA方法[14]对Y型声屏障进行优化分析，目标函数为6个计算点处的平均声压级，设计变量为顶端角度和长度。约束条件如下：

(51)

设定设计变量初始值为

图12依次给出设计变量和目标函数随迭代次数变化关系，可以看出在第500次迭代后收敛的比较好。优化后的角度值为θ1=27°和θ1=￥2°，目标函数优化值为60.42dB。在优化分析过程中，为了保证顶端结构可被用于边界元离散，强迫设置边长大于0.1m，观察图12中的第二幅图可以发现，顶端边长的优化值趋向于满足l1=0和l2=3，该结果表明在满足一定材料成本限制时，Y型声屏障的优化后顶端结构外形是半Y型，实际上半Y型声屏已被广泛应用于工程应用中。图13给出半Y型声屏障周围的声压级分布，频率为200Hz。使用基于MMA的优化分析时，在每一个迭代步中都需要进行目标函数的敏感度计算，实际上这一步是十分耗时的，本文采用伴随变量法进行敏感度计算并使用FMM算法加速优化分析，使得本文发展算法对大规模实际问题有较高的应用潜力。

图12 设计变量及目标函数随迭代次数变化Fig.12 Values of design variables and objective function with number of iteration

图13 半Y型声屏障周围的声压级分布Fig.13 SPL contour plot of half-Y noise barrier

4 总结

本文首先给出基于Burton-Miller方法的不同类型单元无奇异计算表达式，然后给出基于伴随变量法的声学敏感度方程。通过带有解析解的数值算例对比不同类型边界单元的计算精度，发现非连续单元比连续单元表现的更有效，同时非连续单元的优化节点位置参数值非常接近勒让德多项式的零值。最后通过多极算射例子和声屏障算例验证本文发展算法的有效性，尤其在结构优化分析上的高应用潜力。

致谢：国家自然科学基金项目(11172291)；高等学校博士学科点专项科研基金项目(20133402110036)；中国科学技术大学校青年创新基金项目(WK2090000007)。