耦合边界元法和等效源法的稳健CHIEF法

包英超, 向 宇, 陈 洁, 石梓玉

(1. 广西科技大学 广西汽车零部件与整车技术重点实验室,广西 柳州 545006;2. 广西科技大学 机械与汽车工程学院,广西 柳州 545006;3. 合肥工业大学 机械工程学院,合肥 230009)

在工程声学问题计算中,边界元法(boundary element method,BEM)[1-2]由于其计算精度高、系数矩阵主对角占优、求解稳定好、降维求解、适合无限域/半无限域问题的优点,被广泛应用于飞机、列车、船舶、建筑以及环境的噪声预测和控制[3-5]。

然而,常规边界元法(conventional boundary element method, CBEM)存在两个缺陷:(1) 解的非唯一性问题[6-7];(2) 奇异积分问题[8]。目前,有两种经典方法解决Helmholtz边界积分方程的非唯一性问题。第一种是组合亥姆霍兹积分方程公式(combined helmholtz integral equation formulation,CHIEF)法,先通过在内域中选取若干个CHIEF点形成补充方程,再与常规边界元方程联立形成超定方程组,利用最小二乘法求解该超定方程即可获得唯一解。该方法原理简单,易于实现。但是,若CHIEF点位于内部本征问题的节点线上时,该方法失效。另一种是Burton-Miller方法,即将常规边界积分方程与其法向导数方程线性耦合,当耦合系数取为复数时就可得到全频段唯一解。然而,法向求导增加了积分核函数的奇异性,在计算中需处理强奇异积分和超强奇异积分。虽然可以采用区间分割法[9]、解析与半解析法[10-11]和非线性变量转换法(包括双曲正弦变换[12]、距离变换[13]和指数变换[14-15])等方法消除或计算奇异性,但这些方法大多处理过程较繁琐或者得到的边界积分方程形式复杂(尤其对曲面或曲线单元),增加了应用实施的难度。

对于Helmholtz外声场问题的求解,Ochmann提出一种源强模拟技术[16-17],其基本思想源于叠加原理和Helmholtz外问题解的唯一性定理,即虚拟源强叠加产生的外声场与振动体自身产生的辐射或散射声场等效。因源强的作用位置、数量及类型的不同,又可分为单点-多极子法、多点-单极子法、多点-偶极子法、多点-多极子法等,在应用中这类方法通称等效源法[18-19](equivalent source method, ESM)。由于等效源点(虚拟源强)与真实边界上的配点不重合,因此ESM完全避免了BEM中的奇异积分计算,且可以很容易采用单极子源和偶极子源线性组合的方法来获得唯一解。此外,ESM的插值形函数在实边界的分布精确满足控制方程[20],其计算精度高于BEM。但ESM矩阵通常具有较大的条件数,呈病态性,其数值求解稳定性不如BEM。尤其是当等效源位置设置不当时,严重病态的系数矩阵会导致求解过程对输入误差极度敏感,很难得到稳定的数值解[21]。

综上所述,CBEM存在解的非唯一性问题;CHIEF法有内点选择失效的风险;Burton-Miller法存在超强奇异积分处理和计算精度下降的问题;ESM存在系数矩阵病态,求解稳定性和鲁棒性较差的问题。为此,本文提出一种既能克服解的非唯一性问题,又无需奇异积分处理,且具有高计算精度和稳定性的耦合CHIEF法。文中首先简要阐述了CBEM和ESM的基本理论,然后详细地推导了耦合CHIEF法的理论公式并分析了其作用机理,最后利用经典的声辐射和声散射算例验证了所提方法的有效性和优越性。

1 基本理论

1.1 边界元法

考虑一个三维空间声辐射外问题,如图1所示。其中,S为结构体的实际表面边界,D-为结构体的内部域,D+为结构体的外部域(分析域)。在区域D+内,声场控制方程为非齐次Helmholtz方程

图1 边界元法示意图Fig.1 Schematic diagram of the boundary element method

(∇2+k2)p(r)=-f(r) ,r∈D+

(1)

式中:p(r)为r处的声压;k=ω/c为波数;ω为角频率;c为声速;f(r)为空间分布源。应用格林第二恒等式可以将方程(1)表示为常规边界积分方程(conventional boundary integral equation,CBIE)为

(2)

式中:r和rS分别为场点和源点位置矢量;pin(r)为场点处的入射声压;nS为结构表面S上背离分析域的法向;C(r)为奇异性系数,当场点r∈S且S为光滑边界时,C(r)=0.5;当场点r∈D+时,C(r)=1;当场点r∈D-时,C(r)=0。G(r,rS)为三维自由空间格林函数

将场点在实边界面S上进行配点[C(r)=0.5],离散式(2)可得到CBEM方程

ASP=BSQ+Pin

(3)

以常数单元为例,式(3)中,AS、BS的系数矩阵元素分别为

(4)

其中,

(5)

当场点和源点位于同一个单元内时,式(4)中的Green函数及其法向导数分别具有O(1/R)的弱奇异性和O(1/R2)的强奇异性,因此,CBEM在形成系数矩阵的过程中需要处理对角元素的奇异积分。

通常利用边界积分方程(2)求解外部声学问题时,对应的内部声学问题的特征频率处不能保证解的唯一性。Burton等提出一种边界积分方程及其法向导数方程线性组合的方法,消除解的非唯一性问题。式(2)的法向导数边界积分方程(normal derivative boundary integral equation,NDBIE)为

(6)

式中,n为边界S上场点r处的内法向,如图1所示。将式(6)与式(2)进行线性组合并离散,可得Burton-Miller法的求解方程

(7)

式中,ABM,BBM的系数矩阵元素分别为

式中,α是组合系数[22],当α为复数时,式(7)可在全频域内克服解的非唯一性问题,但由于边界积分方程的法向求导引入了超强奇异积分,增加了计算的难度。

1.2 等效源法

ESM是在声源内部布置一系列满足Helmholtz方程和Sommerfeld辐射条件的虚拟源强来模拟声源向外辐射的声场。为了便于计算,可采用简单点源,并将这些点源布置在与实边界共形并向内回缩距离为d的内部虚拟面上,点源的数量和位置通常与实边界配点一一对应,如图2所示。其中,rE为等效源位置,nE为等效源的单位法向。根据等效源的思想,其边界场点r处的声压可表示为

图2 等效源法示意图Fig.2 Schematic diagram of the equivalent source method

(8)

P=AEW+Pin

(9)

Q=BEW+Qin

(10)

式中,W为等效源源强列向量。AE和BE系数矩阵元素分别为

由于实边界上的场点和等效源点不重合,因此,AE和BE中不存在奇异积分,其计算难度远低于BEM。此外,在ESM中,实边界的声压是由内部各个源点辐射的声压叠加形成,其分布精确满足Helmholtz方程[23];而在BEM中,无论采用常数单元还是高阶单元,其实边界声压分布均为近似假设,不满足Helmholtz方程。因此,在相同的单元数和单元类型情况下,ESM的计算精度高于BEM。

2 耦合CHIEF法

2.1 公式推导

本文借鉴CHIEF法的思想,将ESM方程作为补充方程,并与CBEM矩阵方程(3)组合,建立一种耦合算法。

联立式(9)和式(10),消去虚拟源强列向量W,经整理后可得

(11)

在ESM中,矩阵AE和BE通常是病态的。为了获得较高的计算稳定性,一般需要避免直接计算单个矩阵的逆。分别利用AE和BE左乘式(11),使两个矩阵成对出现,得到

(12)

(13)

将式(12)和式(13)分别与式(3)联立,可得以下两种组合方式:

组合1

(14)

组合2

(15)

方程(14)和方程(15)即为本文提出的耦合CHIEF法的方程。本文的组合方程虽然在形式和原理上与CHIEF法相似,但ESM方程和CBEM方程是不相关的两个独立方程,其中,CBEM方程对应实边界面,ESM方程对应虚拟等效源面,两组方程有着不同的本征频率,所以,它们各自的方程不会在同一波数下同时失效。因此,相对于传统CHIEF法,本文方法具有更高的可靠性和鲁棒性。

2.2 CBEM奇异矩阵的间接计算

Xiang等[24]利用系数矩阵的等价关系,提出了一种奇异矩阵的间接计算方法,用于间接求解CBEM的奇异矩阵BS。将CBEM的方程(3)改写为

(16)

根据解的唯一性原理,通过对比式(12)和式(16)右侧第一项的系数矩阵,可进一步得到CBEM和ESM系数矩阵之间的等价耦合关系

ASAE=BSBE

(17)

矩阵BS可由式(17)间接计算

(18)

使用式(18)计算BS不需要直接计算奇异积分,且间接矩阵BS包含了ESM系数矩阵AE和BE的信息。因此,在一定程度上还可进一步继承ESM的高精度优点。

3 数值算例

3.1 球面模型的声辐射和声散射

通过计算几个典型球面模型的声辐射和声散射问题,验证本文方法的有效性和优越性。在算例中均采用三角形平面常数单元离散真实边界,虚拟边界缩进距离d取为真实边界平均单元大小的1.0倍,可获得与真实边界离散单元数量相同且共形的虚拟球面。球面模型的平面三角形离散如图3所示(990个单元)。实边界离散为青色网格,等效源布置在虚拟面上(取三角形单元质心)。

(1) 声辐射

情形1:半径为a的脉动球距离球心r处的声压解析值为

(19)

情形2:半径为a,径向振速为vncosφ的振动球对距离球心r处的声压解析值为

(20)

式中:vn为脉动球表面质点法向振速幅值;ρ为介质的平均密度;c为介质中的声速;φ为场点矢径与z轴正方向的夹角。

(2) 声散射

设有沿着z轴正向入射的平面波pin=eikz被一个中心位于坐标原点、半径为a的球体阻挡,考虑如下三种球面边界:

情形1:当球体表面为刚性边界时(即表面总振速∂(psc+pin)/∂r|r=a=0),球面无量纲散射声场的解析值为

(21)

情形2:当球体表面为柔性边界时(即表面总声压(psc+pin)|r=a=0),球面无量纲散射声场的法向导数的解析值为

(22)

情形3:当球体表面为阻抗边界时(即表面总声压和表面总振速满足如下关系: [psc+pin+ξ(∂psc+∂pin)/∂r]|r=a=0),球面无量纲散射声场的解析值为

(23)

3.2 本文方法的有效性

在下面的算例中,球面声辐射与声散射问题均仅计算真实边界上的表面声压或振速,取r=a=1 m,声阻抗取ξ=10,计算波数范围设置为k=0～10。

对于平面常数单元,在形成CBEM和Burton-Miller的系数矩阵时,非对角元素采用25点Gauss-Legendre积分计算,AS的主对角元素为0.5,BS、ABM和BBM的主对角元素采用无奇异的线积分计算[25-26]

(24)

(25)

(26)

式中:R(θ)为单元边界到配点ri的距离;θ为以配点ri为原点的局部极坐标系中的角度变量。但需要注意的是,式(24)、式(25)和式(26)的奇异积分解析算法仅适用于平面单元,曲面单元一般采用奇异相消法,计算更为繁琐。

3.2.1 解的非唯一性及计算误差比较

首先对比CBEM、Burton-Miller法和本文方法的计算误差。其中,本文方法分别考虑了替换BS矩阵前后的组合1和组合2。结果如图4所示,声压相对误差ε定义为

图4 球面模型在不同边界条件下的表面声压和声压法向导数的误差对比Fig.4 Error comparison of surface sound pressure and normal derivative of sound pressure in spherical model under different boundary conditions

(27)

由图4可知,CBEM在本征频率处失效,Burton-Miller法和本文所提出的方法在整个计算波数内均能有效克服解的非唯一性。从计算误差来看,Burton-Miller法的误差最大,其在非本征频率处的误差高于CBEM,说明超强奇异积分核函数影响了计算精度。而本文方法的误差明显低于Burton-Miller法,尤其是替换BS之后,其在整个波数范围的误差均低于CBEM。此外,虽然在辐射问题中,组合1和组合2的计算误差相近,但在声散射问题中,组合2的计算精度要显著高于组合1。因此,在下文中将采用组合2+替换BS矩阵的方式进行计算比较。

3.2.2 条件数比较

CBEM、ESM、Burton-Miller法和本文方法的系数矩阵条件数,如图5和图6所示。由图5和图6可知,CBEM在本征频率处出现奇异性的跳跃,其它三种方法均为平滑曲线。进一步表明了ESM、Burton-Miller和本文所提方法均能克服解的非唯一性。

图5 声压系数矩阵的条件数对比Fig.5 Comparison of the coefficient matrix condition numbers of sound pressure

图6 声压法向导数系数矩阵的条件数对比Fig.6 Comparison of the coefficient matrix condition numbers of Normal derivative of sound pressure

在图5的声压系数矩阵条件数中,本文方法的条件数虽然略大于CBEM和Burton-Miller中的条件数,但数值在30以内,但远小于ESM。同样,在图6的声压法向导数系数矩阵条件数中,本文方法的条件数最小,其次是Burton-Miller和CBEM,ESM的条件数则大得多。由此可见,本文方法的求解稳定性与BEM相当,显著优于ESM。

3.3 计算精度和计算效率对比

为了进一步的对比本文方法和Burton-Miller法的计算精度,在不同单元数下分别采用这两种方法计算振动球辐射和柔性球散射,其相对误差如图7和图8所示。由图7和图8可知,本文方法的计算误差要小于Burton-Miller方法。且随单元数量的增加,本文方法的收敛速度略快于Burton-Miller方法。

图7 不同单元数振动球表面声压误差对比Fig.7 Comparison of dimensionless sound pressure errors on the surface of the oscillating sphere with different number of elements

图8 不同单元数柔性球表面声压误差对比Fig.8 Comparison of scattered sound pressure normal derivative errors on the flexible spherical surface with different number of elements

再来对比计算效率。以阻抗球散射为例,表1中记录了Burton-Miller法、BEM、ESM和本文方法的计算时间,包括不同单元数下生成系数矩阵的时间,以及求解计算整个问题的总时间。其中,求解计算整个问题的总时间包括了输入计算参数、生成系数矩阵以及求解矩阵方程的时间。

表1 不同方法对应的计算时间Tab.1 The calculation time corresponding to the different methods

由表1可以看出,Burton-Miller法的矩阵计算时间略大于BEM的计算时间,而ESM的矩阵计算时间要远小于Burton-Miller法的计算时间。其次,对比BS矩阵替换前后的计算时间,可以看出,间接替换计算形成BS所需的时间远小于直接计算BS的时间。最后可以发现,本文方法的系数矩阵的计算时间及问题求解的整体时间要显著小于Burton-Miller法的计算时间。其主要原因是ESM中的AE和BE矩阵在形成的过程中不需要积分计算和无需处理对角元素的奇异积分,因此,其计算效率远高于Burton-Miller法或BEM中需要数值积分的矩阵。且替换计算形成BS矩阵是由AS,AE,BE的简单运算得到,较于直接计算BS,无需通过数值积分逐个计算矩阵所有元素,故所需时间也非常短。综上所述,本文方法在计算效率方面相比于Burton-Miller方法具有一定优势。

3.4 不同结构的声辐射和声散射问题的计算结果

为了进一步考察本文方法适应任意形状结构体的情况,本节采用CBEM、Burton-Miller法以及组合2+替换BS计算光滑连续边界的两端带球帽的圆柱模型、环面模型和“音叉”模型,模型及参数如图9～图11所示。在计算三种模型时,将等效源布置在由真实边界向内共形回缩的虚拟面上,缩进距离d为1.0倍的真实边界平均单元大小。考虑如下五种边界条件:脉动辐射、振动辐射、刚性散射、柔性散射和阻抗散射,然后计算模型上P点的声压及其法向导数。由于两端带球帽的圆柱模型、环面模型和“音叉”模型难以获得解析解,因此,将CBEM的计算结果作为参考和比较。结果如图12～图14所示。

图9 两端带球帽的圆柱面的三维离散模型及三维视图参数Fig.9 3D discrete model and 3D view parameters of cylinder with ball cap at both ends surface

图10 环面的三维离散模型及三维视图参数Fig.10 3D discrete model and 3D view parameters of torus surface

图11 “音叉”面的三维离散及三维视图参数Fig.11 3D discrete model and 3D view parameters of tuning fork surface

图12 在不同边界条件下两端带球帽的圆柱面模型表面P点的无量纲声压幅值及无量纲声压法向导数幅值Fig.12 Dimensionless sound pressure amplitude and normal derivative of sound pressure amplitude of P point on the surface of cylindrical model with ball cap at both ends under different boundary conditions

图13 在不同边界条件下环面模型表面P点的无量纲声压幅值及无量纲声压法向导数幅值Fig.13 Dimensionless sound pressure amplitude and normal derivative of sound pressure amplitude of P point on the surface of torus under different boundary conditions

图14 在不同边界条件下“音叉”面模型表面P点的无量纲声压幅值及无量纲声压法向导数幅值Fig.14 Dimensionless sound pressure amplitude and normal derivative of sound pressure amplitude of P point on the surface of tuning fork under different boundary conditions

由图12～图14可知,在特征频率下,本文方法有效克服解的非唯一性。在非特征频率下,本文所提方法的计算结果与CBEM方法的计算结果吻合较好,说明本文方法可以用于不同结构模型的声辐射及声散射问题计算。

4 结 论

针对边界元法的非唯一性问题和奇异积分问题,本文借鉴CHIEF法思想,将CBEM和ESM方程联立形成超定方程组,并利用两组方程系数矩阵的耦合等价关系进行替换计算,提出了一种在全波数均具有唯一解的高精度、高稳定性的声学耦合CHIEF法。文中对所提方法进行了详细的推导和阐述,通过典型三维声辐射和声散射球源算例进行了验证,得到结论如下:

(1) 由于采用单极子和偶极子线性组合的ESM方程作为补充方程,所提方法在计算波数域内均可获得唯一解,有效避免了传统CHIEF法中的内点方程失效问题。

(2) 利用奇异矩阵的间接计算,不仅避免了直接对奇异积分的计算,而且还一定程度上继承ESM高精度的特性。且所提方法仅需对一个系数矩阵进行积分计算,与Burton-Miller法相比,减少了三个积分核函数的数值积分计算,可以大幅减少计算时间,提高计算效率。

(3) 本文方法的计算精度远高于Burton-Miller法,且具有更好的收敛性。特别是对于球面模型的声辐射问题,在相同计算精度下,本文方法极大地减小了计算规模。

(4) 本文方法的系数矩阵条件数远低于ESM,其条件数值较好,具有良好的求解稳定性。

此外,文中虽然只考虑了平面常数单元,但所提出的耦合算法及奇异矩阵间接替换的思想完全可以方便地应用到曲面或曲线单元中。