三维声腔的模态分析对压缩机和消声器等机械设备的声学设计有重要的指导意义，对于比较规则的三维声腔，如矩形、圆柱形和球形等，其特征解可以用解析法^[1]求得。有限元法(finite element method，FEM)^[2]被广泛应用于复杂结构内部声学的计算与分析。FEM中的网络离散过程舍弃了几何模型的精确信息，导致其精度下降，这种误差主要来源于简单的多项式形函数。Huttunen等^[3]进行了单位分解有限元(partition of unity FEM，PUFEM)与超弱变分(ultra-weak variational formulation，UWVF)2种方法的对比研究，通过分析波在二维管道中的传播，指出UWVF在高频范围具有更好的计算精度。广义有限元法^[4]通过引入特定的局部逼近形函数减弱了“污染效应”，其思想是在节点处引入广义自由度，随着自由度的增加，计算量逐渐增大。Grosu等^[5]采用间断富集法(discontinuous enrichment method, DEM)研究了二维区域的内部声学问题，DEM将非协调富集函数与FEM中的多项式形函数之和作为近似解，与PUFEM和UVWF相比具有更高的精度与计算效率^[6]。Desmet^[7]以Trefftz法为基础提出了波函数法(wave based method，WBM)，并进行了结构-声耦合系统的稳态动态分析。WBM在整个求解区域内采用满足控制方程的波函数描述场变量，与FEM相比，WBM不需要划分单元，具有更高的收敛率，但保证其收敛的条件是计算域为凸形域^[8]，这在一定程度上限制了WBM的应用范围。He等^[9]采用alpha有限元(α-FEM)研究了一维及二维声学问题的无色散分析；α-FEM通过改变参数α的值来控制数值模型的刚度，由于波数、网格等多种因素的影响，难以选取最佳的α值。Missaoui等^[10]提出了积分模态法，并研究了二维不规则声腔的声学特性；该方法将整个声腔分割成多个子声腔，通过增加子声腔的数量来提高结果的准确性。此外，一些学者采用级数展开法，如Taylor级数^[11]、Fourier级数^[12]、Chebyshev级数^[13]对声压进行计算，研究了封闭空间的声学特性。

在结构设计中，计算机辅助设计(computer aided design，CAD)与计算机辅助工程(computer aided engineering，CAE)已经密不可分，但两者采用不同的数学语言描述几何模型与计算模型，这种不一致性割裂了产品的设计与分析，还会引起几何误差。为了克服这些障碍，Hughes等^[14]提出了等几何分析(isogeometric analysis，IGA)，将非均匀有理B样条(non-uniform rational B-spline，NURBS)基函数作为FEM中的形函数，保证了几何模型的精确性，从源头消除了几何误差，实现了CAD和CAE数学描述的统一。近十几年来，等几何分析已成功应用于板壳振动^[15]、流体力学^[16]等多个领域。

结构形状对域内声场分布有明显影响，有限元法缺乏精确的几何模型，会导致求解失真，且网格划分过程十分冗繁，网络修改难度很大。NURBS基函数能灵活建立复杂结构的精确模型，IGA通过节点插入与基函数升阶来增加几何模型控制点的个数，避免了繁琐的网格划分过程，具有较高的潜在应用价值。目前等几何分析在声学领域^[17-18]的应用较少，等几何分析在声学领域的研究有待深入。

本文介绍了NURBS的基本概念和三维内部声场的等几何分析，以截面为椭圆或椭圆环声腔为例，计算了其波数和模态振型。将数值计算结果与解析解对比，验证了等几何分析法的正确性及收敛性，分析了几何参数对椭圆声腔固有特性的影响。

1 非均匀有理B样条

给定一个单调不减的节点矢量E = {ξ₁, ξ₂, …, ξ_i, …, ξ_n+p+1}, 其中ξ_i为第i个节点，n和p分别为基函数的个数与阶次，B样条基函数定义为：

$ \left\{\begin{array}{l} p=0: N_{i, 0}(\xi)=\left\{\begin{array}{ll} 1, & \xi_{i} \leqslant \xi \leqslant \xi_{i+1} \\ 0, & \text { 其他 } \end{array}\right. \\ p \geqslant 1: N_{i, p}(\xi)=\frac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}} N_{i, p-1}(\xi)+ \\ \ \ \ \ \ \ \frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}} N_{i+1, p-1}(\xi) \end{array}\right. $

(1)

图 1给出了定义在E={0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1}上的二次B样条基函数。一般地，B样条基函数仅在节点0与1处为插值函数。

对每个B样条基函数乘以对应的权因子ω_i, i = 1, 2, …, n, 并除以权函数W(ξ)，可以得到B样条基函数的有理形式，即NURBS基函数：

$ R_{i}^{p}(\xi)=\frac{\omega_{i} N_{i, p}(\xi)}{W(\xi)}=\frac{\omega_{i} N_{i, p}(\xi)}{\sum\limits_{i=1}^{n} \omega_{i} N_{i, p}(\xi)} $

(2)

张量积形式的二维及三维NURBS基函数为：

$ R_{i, j}^{p, q}(\xi, \eta)=\frac{\omega_{i, j} N_{i, p}(\xi) M_{j, q}(\eta)}{\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m} \omega_{i, j} N_{i, p}(\xi) M_{j, q}(\eta)} $

(3)

$ \begin{gathered} R_{i, j, k}^{p, q, r}(\xi, \eta, \zeta)= \\ \frac{\omega_{i, j, k} N_{i, p}(\xi) M_{j, q}(\eta) L_{k, r}(\zeta) }{\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m} \sum\limits_{k=1}^{l} \omega_{i, j, k} N_{i, p}(\xi) M_{j, q}(\eta) L_{k, r}(\zeta)} \end{gathered} $

(4)

式中：ω_{i, j}、ω_{i, j, k}是二维和三维权因子；M_{j, q}(η)、L_{k, r}(ζ)分别为定义在节点矢量H =(η₁, η₂, …, η_j, …, η_m+q+1)和Z =(ζ₁, ζ₂, …, ζ_k, …, ζ_l+r+1)上的B样条基函数。NURBS曲面和实体分别定义为：

$ \boldsymbol{S}(\xi, \eta)=\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m} R_{i, j}^{p, q}(\xi, \eta) \boldsymbol{B}_{i, j} $

(5)

$ \boldsymbol{S}(\xi, \eta, \zeta)=\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{m} \sum\limits_{k=1}^{l} R_{i, j, k}^{p, q, r}(\xi, \eta, \zeta) \boldsymbol{B}_{i, j, k} $

(6)

式中B_{i, j}和B_{i, j, k}是曲面和实体的控制点。

2 三维声场等几何分析 2.1 控制方程

考虑边界为刚性壁面的三维内部声场，其控制方程及边界条件为：

$ {\nabla ^2}p\left( {x,y,z} \right) + {k^2}p\left( {x,y,z} \right) = 0,区域\mathit{\Omega} 内 $

(7)

$ \partial p/\partial n = - {\rm{j}}{\rho _0}\omega {u_n},边界\mathit{\Gamma} 上 $

(8)

式中：$ {\nabla ^2} = {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2} + {\partial ^2}/\partial {z^2} $为拉普拉斯算子；p为声压；k=ω/c₀为波数；n为边界Г的外法线方向。

2.2 基于NURBS的离散方程

等几何分析的核心思想是将精确描述几何的NURBS基函数作为有限元中场变量的形函数。因此，式(7)中的声压可表示为：

$ p^{h}(\xi, \eta, \zeta)=\sum\limits_{A=1}^{m \times n \times l} R_{A}(\xi, \eta, \zeta) p_{A}=\boldsymbol{R}^{\mathrm{T}} \boldsymbol{p} $

(9)

式中：m、n及l分别代表ξ、η和ζ方向的基函数个数；p_A为控制点A处的声压；R=[R₁R₂…R_m×n×l]^T为基函数向量；p=[p₁p₂…p_m×n×l]^T为声压向量。

采用伽辽金法，将NURBS基函数乘以式(7)并在区域Ω上积分，得到：

$ \int_{\varOmega} \boldsymbol{R}\left(\nabla^{2} p^{h}+k^{2} p^{h}\right) \mathrm{d} \varOmega=0 $

(10)

应用格林公式，式(10)可以写为：

$ \int_{\varGamma} \boldsymbol{R} \frac{\partial p^{h}}{\partial n} \mathrm{~d} \varGamma-\int_{\varOmega} \nabla \boldsymbol{R} \cdot \nabla p^{h} \mathrm{~d} {\varOmega}+k^{2} \int_{\varOmega} \boldsymbol{R} p^{h} \mathrm{~d} \varOmega=0 $

(11)

将式(8)、(9)代入式(11)，可得声学问题的特征方程为：

$ \left(\boldsymbol{K}-k^{2} \boldsymbol{M}\right) \boldsymbol{p}=0 $

(12)

式中K、M为整体刚度与质量矩阵，表示为：

$ \boldsymbol{K}=\int_{\varOmega}(\nabla \boldsymbol{R}) \cdot(\nabla \boldsymbol{R})^{\mathrm{T}} \mathrm{d} \varOmega $

(13)

$ \boldsymbol{M}=\int_{\varOmega} \boldsymbol{R} \boldsymbol{R}^{\mathrm{T}} \mathrm{d} \varOmega $

(14)

3 数值结果与分析 3.1 椭圆台形声场特性分析

考虑截面连续变化的椭圆台形声场，如图 2所示，两端椭圆的长短半径分别为a₁、b₁、a₂、b₂，椭圆台高度为h。采用NURBS建立其精确的几何模型，初始节点矢量为E=H=(0, 0, 0, 1, 1, 1)，Z=(0, 0, 1, 1)，初始控制点数为3×3×2。当a₁=a₂，b₁=b₂时，声腔为椭圆柱，其特征值问题存在解析解。取尺寸参数a₁=a₂=0.15 m，b₁=b₂=0.09 m，h=0.4 m，声速c₀=1 m/s，前8阶波数的计算结果如表 1所示，随着基函数阶次的提高和控制点数的增加，等几何分析的计算结果收敛并与解析解^[19]吻合。基函数的阶次越高，前8阶波数收敛的速率越快。采用10×10×10(自由度为1 000)和12×12×12(自由度为1 728)的计算结果优于自由度为1 361和3 672时有限元^[19]得到的结果，因此等几何分析与传统有限元相比具有更高的准确性和计算效率。

接下来分析离心率对椭圆台声腔固有特性的影响，截面的连续变化增加了解析法求解的难度。取a₁=0.10 m，a₂=0.15 m，b₂=0.09 m，h=0.4 m，基函数阶次p=q=r=3，控制点10×10×10，不同离心率下椭圆台声腔的前8阶波数如表 2所示。可以看出，随着离心率逐渐增大，前8阶波数呈上升趋势，这是因为离心率的增大使b₁减小，椭圆台形声场的几何变化导致其波数增大。此外，离心率在0~0.4波数的增幅小于离心率在0.4~0.8波数的增幅，原因是离心率较小时，离心率对短半径值的影响不大。当e₁=0.4时，椭圆台形声场的模态振型如图 3所示。

3.2 同心椭圆环柱形声场特性分析

同心椭圆环柱形声场的固有特性，柱体截面如图 4所示。三维柱体的几何参数为a₁=0.10 m，b₁=0.06 m，a₂=0.15 m，b₂=0.09 m，h=0.4 m。采用NURBS建立精确的几何模型，初始节点矢量为E=(0, 0, 0, 0.25, 0.25, 0.5, 0.5, 0.75, 0.75, 1, 1, 1)，H=Z=(0, 0, 1, 1)，周向、径向和轴向的初始控制点个数分别为9、2、2，记为9×2×2。通过基函数升阶和节点插入后，不同网格下的计算结果如表 3所示，NURBS基函数阶次为p=q=r=3。可以看出随着控制点数的逐渐增加，计算结果逐渐收敛，当控制点个数达到25×7×23后，继续细化网格对结果的影响很小，所以该控制点网格将用于本节后续的计算与分析。图 5给出了椭圆环柱形声场的前8阶模态振型，可以看出第1阶与第6阶为轴向模态，第2阶与第3阶为周向模态。

接下来分析椭圆内环的离心率对波数的影响，取a₁=0.06 m，椭圆外环及柱形高度的尺寸参数同上，图 6给出了离心率e₁对椭圆环柱形声腔波数的影响。对于给定内环椭圆长半径的柱形声腔，内环椭圆短半径随着离心率e₁的增大而减小，第1阶与第5阶波数未出现明显变化，第2、3、7阶波数逐渐增大，第4、6、8阶波数逐渐减小。

3.3 椭圆偏心环柱形声场特性分析

最后分析椭圆偏心环柱形声场的固有特性，柱体截面如图 7所示，几何参数为内圆半径r=0.06 m，椭圆长半径a=0.15 m，椭圆短半径b=0.09 m，柱体高度h=0.4 m。采用NURBS基函数精确建立此声腔初始几何模型所用到的节点矢量与上节同心椭圆环柱相同。取NURBS基函数阶次为p=q=r=3，控制点个数为25×7×23，表 4给出了不同圆心位置下此声腔的前8阶波数。从表 4可以看出，改变内环圆心O₂的位置会引起声腔波数的变化。

在图 7所示的柱体截面所在直角坐标系中分析圆心位置的影响，当圆心O₂从原点(0, 0)向右沿着x轴方向移动到(0.02, 0)，第2、3、4、6、7阶波数逐渐增加；然后圆心O₂从(0.02, 0)向上移动到(0.02, 0.02)，再向左移动到(0.01, 0.02)，在此过程中第2、3、7阶波数逐渐减小，第4阶波数继续增加，第6阶与第8阶波数先增大后较小。不同圆心位置下椭圆偏心环柱形声场的第1阶与第3阶模态振型如图 8所示。由此可见，不同圆心位置下的第1阶模态振型相似；然而对于第3阶模态振型，内部偏心圆导致了几何模型上的不对称，故其声场分布存在局部变化。

4 结论

1) 对于椭圆台形声腔，其体积随着离心率增加而减小，导致前8阶波数逐渐增大。

2) 对于同心椭圆环柱形声腔，随着内环椭圆离心率增大，第1阶与第5阶波数变化不明显，第2、3、7阶波数逐渐增大，第4、6、8阶波数逐渐减小。

3) 由于几何的不对称，偏心圆所在位置对椭圆偏心环柱形声场的前8阶波数影响较为复杂。

4) 等几何分析具有良好的收敛性和计算精度，适用于复杂声场固有特性的预测与评估。