汽车的噪声与振动性能(NVH)是评价汽车品质的重要指标之一，在产品设计阶段初期用数值方法准确预测车内辐射噪声，从而进行具有量化依据的噪声优化设计，对于提高产品竞争力具有重要意义．当进行汽车噪声与振动性能分析计算时，通常先简化模型和划分网格，不仅耗费大量精力，而且可能丢失几何细节信息，导致计算误差[1]．等几何分析是全新分析理念，通过用几何建模的样条曲线直接进行分析计算，简化了网格划分过程．边界元法(BEM)与等几何分析理念结合形成等几何边界元法(IGABEM)[2]，统一了几何模型和分析模型，为汽车噪声与振动性能分析计算带来了新的发展前景．等几何边界元法的相关研究主要集中在简单几何和连续单元的研究中．对于汽车噪声与振动性能计算问题，车身结构一般由多个曲面拼接构成，在结构的拐角或者曲面拼接处，法向量很可能并不连续．当用连续单元表征主要变量和导数时，若相邻单元的法向量不唯一或者边界条件不同，则共享节点处的边界值便无法确定，这就是边界元计算中的角点问题．数值实验结果表明，若仍然用连续单元来表征主要变量和参数值，则在角点处会得到较大的计算误差[3]；因此，提高等几何边界元法处理复杂几何和混合边界条件下的计算精度，是推动等几何分析在车内噪声计算中广泛应用的必经之路．在等几何边界元法中，关于角点问题的研究较少，Scott等[4]和Marussig等[5]在基于等几何边界元法计算力弹性力学问题时，采用非连续单元表达主要变量的导数-牵引力，但是对于主要变量-位移仍然保留了连续的表达形式．Andrade等[6]针对疲劳裂纹扩展问题中的尖端裂纹，提出了一种强化等几何边界元方法，通过在网格细分时插入重复的节点，增加非均匀有理B样条(non-uniform rational B-splines，NURBS)的不连续性．王杰[7]基于几何的参数空间和物理空间相互独立的思路，用不同的样条插值描述几何场和物理场，但在解决角点问题的同时增加了计算的复杂度．本研究提出一种完全非连续等几何边界元法(dIGABEM)，将传统边界元法中的“双节点”思路引入等几何边界元解决角点问题．本研究是对文献[8]的改进和完善，在前文基础上深入总结并阐述了角点问题的形式，对非连续单元的配置点位置进行了细致研究，改进了非连续等几何边界积分方程，结合数值算例对计算精度和效率深入讨论，并且补充设计了声压响应实验，对算法进行了证明．本研究提出的非连续等几何边界元法不仅实现了几何模型和分析模型的统一，而且解决了角点问题，结合并行算法提高了车内辐射噪声计算的精度和效率．1 声学等几何边界元法1.1　边界积分方程构建理想流体介质中的声学波动方程为∇2U(s,t)-(1/c2)(∂2U/∂t2)=0    (∀s∈Ω), (1)式中：∇2为拉普拉斯算子；s为求解域Ω内的一点，U(s,t)为t时刻s点处的速度势函数；c为声波在流体中的传播速度．对于简谐声场，U(s,t)=ϕ(s)e-iωt，(2)式中：ϕ(s)为仅与位置有关的速度势幅值函数；ω为圆频率．将式(2)代入式(1)，可以得到∇2ϕ(s)+k2ϕ(s)=0    (∀s∈Ω), (3)式中k=ω/c为波数．若将声源置于s点处，则在场点x处的声压响应就是声学问题的基本解．若s和x不重合，则格林函数G(s,x)满足∇2G(s,x)+k2G(s,x)=0    (∀s,x∈Ω)．三维声学问题的格林函数及其导数的形式分别为G(s,x)=eikr/(4πr);∂G(s,x)∂n=eikr4πrik-1r∂r∂n,式中：r=s-x；n为法向量．通过格林第二特征方程可以将亥姆霍兹(Helmholtz)微分方程转化为亥姆霍兹积分方程，从而将三维问题转化为对求解域边界的二维积分．在求解域Ω内的亥姆霍兹边界积分方程可表达为    ∫Ωϕ(x)∇2G(s,x)-G(s,x)∇2ϕ(x)dΩ=∫Γϕ(x)∂G(s,x)∂n-G(s,x)∂ϕ(x)∂ndΓ(x), (4)式中：ϕ(x)为场点x处的声势能；Γ为求解域Ω的边界；∂ϕ(x)/∂n(x)为势能梯度．将式(3)代入式(4)，可以得到    ∫Ω∇2G(s,x)+k2G(s,x)ϕ(x)dΩ=ϕ(x)∂G(s,x)∂n-G(s,x)∂ϕ(x)∂ndΓ(x). (5)为了求解边界上未知的ϕ和∂ϕ/∂n，须要把式(5)改写成只有边界的方程形式．让源点s趋近边界Γ，可以得到如下常规形式的边界积分方程    C(s)ϕ(s)=∫Γ-G(s,x)[∂ϕ(x)∂n]∙ϕ(x)[∂G(s,x)/∂n]dΓ(x), (6)式中系数C(s)取决于s点所处的位置和边界形状，若s位于求解域内，则C(s)=1；若s位于求解域外，则C(s)=0；若s位于光滑的边界表面，则C(s)=0.5．在声学问题中，边界条件一般有以下三种类型：a. 狄利克雷(Drichlet)条件，即部分边界上声学势能已知，ϕ(x)=ϕ¯(x)；b. 诺依曼(Neumann)条件，即部分边界上势能导数已知，∂ϕ(x)/∂n=∂ϕ¯(x)/∂n；c. 罗宾(Robin)条件，即势能导数可由势能的线性方程来表达，∂ϕ(x)/∂n=α(x)ϕ(x)+β(x)．1.2　等几何分析理论1.2.1　非均匀有理B样条构建NURBS是计算机辅助设计(CAD)中的标准建模曲线，也是等几何分析中最常用的曲线．NURBS由B样条曲线演化而来，B样条曲线基函数以节点向量为基础．一维空间节点向量是一系列单调不递减的实数Ξ=ξ1,ξ2,…,ξi,ξi+1,…,ξn+p+1    (ξi∈R),式中i=1,2,…,n+p+1，n为基函数或者控制点的数量，p为曲线阶数．根据Coxde Boor递推公式，B样条曲线的基函数Di,p(ξ)可以表示为               p=0:Di,0(ξ)=1      (ξi≤ξ≤ξi+1),0      (其他);    p∈Z+:Di,p(ξ)=[(ξ-ξi)/(ξi+p-ξi)]Di,p-1(ξ)+[(ξi+p+1-ξ)/(ξi+p+1-ξi+1)]Di+1,p-1(ξ).引入权重因子wi，NURBS基函数Ri,p(ξ)可利用B样条曲线基函数定义为Ri,p(ξ)=Di,p(ξ)wi /∑j=1nwjDj,p(ξ)．1.2.2　核函数的等几何离散形式在等几何边界元法中，分析域Ω的边界被划分为E个不重叠的子域Γe，每个子域Γe上的局部坐标映射定义为Γe={Fe(u,v):u,v∈[0,1]}．等几何分析中的积分计算基于每个节点跨距，例如[ξi,ξi+1)×[ηj,ηj+1)，但是高斯积分的参数空间Y=(ξ¯,η¯)却定义在[-1,1]×[-1,1]上，因此须要定义一个从全局坐标到参数空间的映射实现等几何边界元中的坐标转换，如图1所示．在三维问题中，雅克比矩阵可以表达为J=∂x∂F∂F∂Y，式中：∂x/∂F为从全局坐标映射到局部坐标的雅可比矩阵；∂F/∂Y为从局部坐标映射到参数空间的雅可比矩阵．10.13245/j.hust.230251.F001图1等几何单元的映射空间利用等参元概念，等几何分析中的势能及其导数可表达为基于NURBS基函数的离散形式，即ϕ(x)=∑i=1n∑j=1mRi,jp,q(u(x),v(x))ϕ¯i,j; (7)∂ϕ(x)/∂n=∑i=1n∑j=1mRi,jp,q(u(x),v(x))q¯i,j, (8)式中：ϕ¯为狄利克雷边界条件；q¯为诺依曼边界条件；n和m，p和q分别为u和v方向上控制点的个数及曲线阶数；Ri,jp,q为二维NURBS基函数，    Ri,jp,q(ξ,η)=Ni,p(ξ)Mj,q(η)wi,j/∑i^=1n∑j^=1mNi^,p(ξ)Mj^,q(η)wi^,j^,其中N和M为自由度．2 非连续等几何边界元法2.1　非连续几何和多种边界条件耦合的角点问题在边界元法中，求解边界积分方程须要用对外法向量n(x)的导数，因此要求法向量n(x)存在且唯一．若相邻单元节点处的法向量存在但不唯一，则会出现角点问题．角点问题可以归纳为以下两种形式．a. 几何角点：在几何边界上，当被共享的节点处法线或切线的变化不连续时，该点便被称为几何角点．几何角点通常有两种类型：一种是几何结构的表面不光滑，法向量有多个方向，如图2(a)所示；另一种是当用拉格朗日单元离散几何边界时，相邻单元在交界处C1不可导，如图2(b)所示．b. 物理角点：若相邻单元的边界条件类型不同，或者同类边界条件的参数值不同，被相邻单元共享的节点处便无法确定边界条件的取值，则该节点被称为物理角点，如图2(c)所示．双节点法是边界元法中处理角点问题较为常用的一种方法．如图3所示，在角点处配置两个节点，其几何坐标相同，但分别属于不同的单元，具有不同的边界条件．在这两个节点处分别建立边界积分方程，约束相应的边界条件，从而解决角点问题．当用双节点法处理角点问题时，不增加积分复杂性，并且对网格有很好的适应性．本研究结合双节点法和等几何边界元分析，在角点处生成重复的控制点，将配置点的位置移到单元内部，建立非连续等几何边界单元，在等几何分析的基础上发展能有效处理角点问题的非连续等几何边界元法．10.13245/j.hust.230251.F002图2角点问题10.13245/j.hust.230251.F003图3双节点法(法向量存在且唯一)2.2　非连续单元构建构造如图4所示的单元，在角点处生成两个控制点，其几何坐标相同，但分别属于不同的单元，具有不同的边界条件．在这两个节点处分别建立边界积分方程，分别约束所属单元上的边界条件．非连续单元内部在几何上可以保证C1连续，而且在配置点处的外法向向量n(x)存在且唯一．10.13245/j.hust.230251.F004图4等几何边界单元在非连续等几何边界元中，为了最终获得一个方形矩阵求解系统，配置点的个数和未知量的个数必须相同．为了确保合适数量的配置点，最直接的方法就是将其放在如图4(b)所示的单元内部．当部分或者所有的节点被放置在单元内部而不是单元边界时，这样的单元即被称为非连续单元．由于配置点移到了单元内部，相邻单元不再共享节点，因此当用非连续单元表征几何场和物理场时，不会出现同一个节点上的参数值有多个的情况，从根本上规避了几何角点和物理角点问题．由于配置点不再被相邻的单元所共享，因此这种独立性使得当非连续单元构建边界单元网格和实现自适应网格细化时更加灵活，并且可以实现并行计算．充分利用非连续单元独有的结构特点，将会大大提高计算效率．2.3　配置策略分析在非连续单元中，控制点可能是重合的，即同一个位置可能有多个控制点．若配置点与控制点重合则会产生相同的方程，从而出现秩亏方程系统，因此需要合适的配置方法来决定配置点的位置．在边界元法中，节点的配置策略是一个重要的研究课题．对此，不同研究者持不同的意见[9]：Brebbia等认为将配置点放在高斯点处可以得到最高计算精度；Atkinson和Marburg等认为最佳配置点的位置是正交函数的零点．其中，Marburg基于长杆模型研究了配置点的位置参数对计算精度的影响规律，指出计算频率会影响最佳配置点位置；而Banerjee等认为当配置点均匀放置在参数空间时，计算精度最高．本研究以二阶单元为例，研究配置点位置对于非连续等几何边界单元计算精度的影响，确定最佳配置点位置参数，为后续计算提供依据．图5展示了一个二阶非连续等几何边界单元，引入参数ai和bj配置点位置，为了简化计算，假设配置点对称分布，设a1=a2=a3=a4=ai，b1=b2=…=bj，则bj=(1-2ai)/(p-1)，因此一个参数ai即可决定非连续单元中所有节点的位置．10.13245/j.hust.230251.F005图5二阶非连续等几何边界单元利用一个立方体声腔模型，研究配置点的位置对计算精度的影响．该模型位于笛卡尔坐标系(x,y,z)∈[0,3]3中．假设一个波长为5 m的平面波，沿x方向穿过这个声腔结构．在这个算例中，参数空间里两个方向的初始节点向量为均匀节点向量γ={0,0,0,1,1,1}，p=2，n=3．在位于x=0面上加载狄利克雷边界条件ϕ¯=1，在位于x=3面上加载诺依曼边界条件q¯=-ksin(3k)+ikcos(3k)．诺依曼边界条件q¯=0被加在其他所有剩下的面上．对于该立方体算例，表面势能有解析解ϕ=cos(kr)+iksin(kr)．单纯考察一个点的误差无法全面评估计算结果，因此定义如下的面误差估算公式ε=ϕL2(Γ)-ϕrefL2(Γ) / ϕrefL2(Γ)，(9)式中：ϕref为参考值；ϕL2(Γ)为势能在边界上的L2范数形式，ϕL2(Γ)=∫Γϕ2dΓ．本研究计算均采用式(9)定义的面误差估算方法进行误差评估．以Δa=0.02为间隔，分别计算当频率为100，400，800，1 200 Hz时立方体的表面势能分布．以解析解为参考解，用面误差估算方法评估计算结果，得到误差随着配置点参数变化的关系曲线，如图6所示．10.13245/j.hust.230251.F006图6基于立方体模型的配置方法计算误差比较由图6可以看出：4种频率下的计算误差曲线随着配置点位置参数值的增加逐渐下降又上升，4条曲线的最小值出现的位置稍有不同，但是当a=0.17时计算误差都非常接近最小值．值得注意的是，当位置参数a=0.17时，配置点在参数空间近似呈均匀分布，该结论与Banerjee的观点一致．本节研究仅为后文计算提供合适的配置策略，即后文将依据a=0.17确定配置点位置．2.4　非连续等几何边界积分方程构建当场点x和源点s位于同一单元时，r=s-x趋近于零，被积函数在奇异点附近剧烈变化，高斯积分的计算结果误差较大，产生奇异问题．在式(6)中，等号右边第一项积分为弱奇异积分，第二项积分为强奇异积分．对于建立的非连续等几何边界单元，由于配置点均位于单元内部，因此在配置点处的角系数C(s)都满足C(s)=1/2．根据Liu的研究[10]，当源点s在边界上时，有以下积分特性∫Γ[∂G¯(s,x)/∂n]dΓ(x)=-1/2，(10)式中G¯(s,x)为拉普拉斯方程的格林函数，G¯(s,x)=1/(4πr)．将式(10)代入式(6)，有    -∫Γ[∂G¯(s,  x)/∂n]dΓ(x)ϕ(s)=∫Γ{G(s,  x)[∂ϕ(x)/∂n]-ϕ(x)[∂G(s,  x)/∂n]}dΓ(x). (11)整理式(11)，得到如下只含有弱奇异项的正则化边界积分方程    ∫Γ{[∂G(s,  x)/∂n]ϕ(x)-[∂G¯(s,  x)/∂n]ϕ(x)}dΓ(x)=∫ΓG(s,  x)[∂ϕ(x)/∂n]dΓ(x). (12)把式(7)和式(8)代入式(12)中，得到离散形式的非连续等几何边界积分方程    ∑e=1E∑i=1n∑j=1mPeij(s)ϕeij-∑e=1E∑i=1n∑j=1mP¯eij(s)ϕe^ij=∑e=1E∑i=1n∑j=1mQeij(s)∂ϕeij∂n; (13)    Peij=∫-11∫-11[∂G(ξ¯,η¯)/∂n]∙Reij(ξ¯,η¯)Jeij(ξ¯,η¯)dξ¯dη¯;    P¯eij=∫-11∫-11[∂G¯(ξ¯,η¯)/∂n]∙Re¯ij(ξ¯,η¯)Jeij(ξ¯,η¯)dξ¯dη¯;Qeij=∫-11∫-11G(ξ¯,η¯)Reij(ξ¯,η¯)Jeij(ξ¯,η¯)dξ¯dη¯,式中：ξ¯=ξ¯(u,v)；η¯=η¯(u,v)；e¯为源点s所在单元的索引号；Reij为NURBS基函数；Jeij为等几何空间映射的雅克比矩阵．方程(13)中仍然包含一个弱奇异项，利用Telles坐标转化法即可解决这一问题．依据上节中最佳位置参数确定配置点的位置，将每一个配置点轮流当作源点s，对式(13)进行移项整理，可以得到一个如下形式的方程Hu=Gq，式中：矩阵H包含了式(13)等号左边所有的积分项；矩阵G包含了式(13)等号右边所有的积分项；向量u和q分别包含了ϕ和∂ϕ/∂n的值．将所有未知量移到等号左边，已知量移到等号右边，最终可获得一个如下形式的方程Aμ=b，(14)式中：系数矩阵A为非对称稠密矩阵；向量μ包含了所有的未知量；向量b包含了所有的已知量．式(14)可以由直接法或者迭代法来求解，进而可得到结构边界上的声势能，将s点移到分析域Ω的内部，则C(s)=1，代入式(6)即可求解分析域Ω内任意一点处的声势能．3 数值算例在此利用一个微型客车车身模型验证非连续等几何边界元法计算辐射噪声的正确性．该车身结构如图7(a)所示，蓝色面板代表仪表板，灰色面板代表车窗玻璃，红色面板代表贴有衬里材料的车身面板．该车身的整体尺寸为3.47 m×1.63 m×1.36 m，车内辐射噪声问题可以用三组边界条件来模拟．10.13245/j.hust.230251.F007图7车身模型对于振动的面板，例如仪表板，该面上的边界条件可以用诺依曼条件来表达，即∂ϕ(x)/∂n=-iρ0ωv，(15)式中：ρ0为空气密度；ω=2πc/λ为角频率，c为声速，λ为波长；v为边界表面法向振速的幅值．对于具有完全反射特性的面板，例如车窗玻璃，该面上的边界条件可以用均匀的诺依曼条件来表达，即∂ϕ(x)/∂n=0．对于贴附有吸声材料的面板，例如车身的内部衬里材料，该面上的边界条件可以用罗宾条件来表达，即∂ϕ(x)/∂n=-iρ0ω[ϕ(x)/Z]，式中Z为声阻抗．依据Marburg等[11]的研究，基于实验数据声导纳值Yz和频率f有如下近似关系，Yz=f /f0，式中f0=2.8 kHz，进而由声阻抗和声导纳的倒数关系Yz=1/Z，即可计算出声阻抗Z．本算例中，将第一个诺依曼条件加载在蓝色面上，第二个诺依曼条件加载在灰色面上，罗宾条件被施加在红色面上．取声速c=340.29 m/s，计算频率f =113 Hz，蓝色仪表板的法向振动速度v=1.452 mm/s，空气密度ρ0=1.29 kg/m3．车身边界元网格离散模型和NURBS模型分别如图7(b)和图7(c)所示．对于NURBS车身模型，通过插入节点以扩大参数空间，分别将车身结构扩展为自由度N=500～4 000的模型．采用非连续等几何边界元法计算车内场点A处(1.00 m，0.60 m，1.03 m)的声势能，首先依据第2节研究的配置策略确定配置点位置，然后通过式(12)建立车内辐射噪声问题的边界积分方程，进而求解式(14)得到车身表面势能，最后依据式(6)求得场点A处的势能．将结果与连续边界元法和连续等几何边界元法的计算结果进行比较，验证本文提出的非连续等几何边界元法在计算精度和效率上的优越性．采用广义最小残差法计算最终线性方程组，迭代残差为1×10-3，并且三种方法采用了相同数量的高斯积分点．本研究基于Matlab平台，用“parfor”循环实现并行计算．程序在配置为Dual core CPU 2.2 GHz、内存为8 GiB的个人计算机上执行．以传统边界元精细模型(N=10 088)的收敛计算结果作为参考值，利用式(9)进行误差分析，计算误差和时间分别如表1和表2所示．10.13245/j.hust.230251.T001表1三种方法计算误差方法自由度5001 0001 5002 0002 5003 0003 500BEM77.9027.9014.308.996.224.553.43IGABEM24.007.332.971.681.050.710.52dIGABEM17.603.861.370.670.390.250.1710-310.13245/j.hust.230251.T002表2三种方法计算时间方法自由度5001 0001 5002 0002 5003 0003 500BEM0.120.862.605.7810.5017.5025.90IGABEM0.070.330.721.352.183.204.43dIGABEM0.060.190.370.580.831.111.42102 s由表1可知：非连续等几何边界元法明显比传统边界元法和连续等几何边界元法的计算精度更高，并且随着自由度的增加，收敛速度也更快．由表2可知：连续和不连续等几何边界元法的计算时间远少于边界元法；随着自由度的增加，非连续等几何边界元法在计算效率上的优势变得愈加明显．并行计算使得非连续等几何边界元法的计算时间明显短于连续等几何边界元法，而和传统边界元法相比，其计算效率有数量级的提高．通过构建配置点位于单元内部的非连续单元，在各个配置点上分别建立边界积分方程并约束所在单元的边界条件，非连续等几何边界元法规避了相邻单元边界上法向量或者边界条件的参数无法确定取值的问题，解决了车内辐射噪声计算中由非连续几何和多种边界条件耦合导致的角点问题，提高了计算精度．该非连续单元的构造方法简单，不改变NURBS基函数的形式，并且易于实现并行计算，提高了非连续等几何边界元法的计算效率，而本研究推导的仅含有弱奇异项的正则化边界积分方程则同步提高了算法的计算精度和稳定性．通过比较三种方法的计算精度和时间，证实了非连续等几何边界元理论推演的正确性，展示了非连续等几何边界元法计算车内辐射噪声问题的优势和潜力．4 实验验证本研究基于一个方形箱体设计声压实验，进一步验证文中方法计算结构辐射噪声的正确性．该箱体材料为Q235钢，尺寸为840 mm×700 mm×460 mm．将底板厚度设计为30 mm，四周围板板厚设计为20 mm，从而可以将地板和四周围板视为刚体，忽略其振动对测量结果的影响．上盖板的厚度为2 mm，其四周边长分别用横截面为20 mm×20 mm的正方形压条，通过密封垫片和螺栓紧固在箱体的壁板上．实验在截止频率为50 Hz、本底噪声为20 dB的半消声室内完成，实验设备包括学习管理系统(LMS)数据采集器、功率放大器、方形箱体和声压传感器阵列．选择箱体内坐标为(400 mm，300 mm，100 mm)的场点为声压观测点，用GRAS公司的46AE声压传感器测量该点处的声压级S为实验值；将上表面各节点处的振动速度代入式(15)中作为边界条件，计算观测点处声压级S为计算值，实验值和计算值的结果对比如图8所示．由图8可知：在分析频率范围内，箱体内观测点处的实验值与计算值的变化趋势比较一致，在关键点处两者之间的误差较小，再次证明了采用非连续等几何边界元法计算内部辐射噪声的准确性．10.13245/j.hust.230251.F008图8声压响应实验值与计算值对比5 结论本研究提出了一种车内辐射噪声计算的非连续等几何边界元法，不仅统一了几何模型和计算模型，简化了网格划分过程，而且解决了车内辐射噪声计算中由非连续几何边界和多种边界条件耦合导致的角点问题，提高了计算精度和效率，得到的主要结论如下．a. 利用“双节点”思路建立配置点位于单元内部的非连续单元，研究了配置点的位置对计算精度的影响，得出当采用二阶均匀节点向量时，a=0.17为最佳配置点位置参数．b. 非连续单元的构造方法简单，不改变NURBS基函数的形式，并且易于实现并行计算．利用非连续单元的积分特性建立了正则化边界积分方程，将原核函数的奇异性由强奇异降为弱奇异，提高了计算精度和稳定性．c. 数值算例和实验结果证明了非连续等几何边界元法计算车内辐射噪声的正确性，并且较连续边界元法和连续等几何边界元法计算精度和效率有明显提升，展示了该方法在汽车噪声与振动性能工程应用中的优势和潜力．