基于α有限元法的二维水下声散射计算

晋文超，岳智君，李 威，李耀飞，柴应彬

(1.海军装备研究院，北京 100036; 2.中国舰船研究设计中心，湖北 武汉 430064; 3.华中科技大学 船舶与海洋工程学院，湖北 武汉 430074)

基于α有限元法的二维水下声散射计算

晋文超1，岳智君2，李 威3，李耀飞3，柴应彬3

(1.海军装备研究院，北京 100036; 2.中国舰船研究设计中心，湖北 武汉 430064; 3.华中科技大学 船舶与海洋工程学院，湖北 武汉 430074)

众所周知，由于存在插值误差和污染误差的原因，高波数声学问题的有限元解是不可靠的。为了提高有限元法的求解精度，本文提出了一种新型的α有限元法来求解水下声散射问题。在α有限元法中，首先运用基于点的梯度光滑技术得到点光滑的梯度场，然后对点光滑有限元法进行推导。因此，α有限元模型既包含来自点光滑有限元模型的梯度成分又包含来自标准光滑有限元模型的梯度成分，充分利用了点光滑有限元模型的“过软”特性和标准有限元模型的“过刚”特性。为了处理外声场问题，本文运用了DtN无反射边界条件。数值计算结果表明：与标准有限元法相比，α有限元法在水下声散射计算中具备更高的计算精度和计算效率。

α有限元法；水下声散射；无界域；无反射边界条件

声学计算在航空、航海等领域有着广泛的应用前景；有许多数值方法可以用来求解声学问题，其中有限元法是各种求解声学问题的数值算法中运用最为广泛的算法之一。然而众所周知的是，标准有限元难以在高频域内提供较为可靠的计算结果[1-3]。

为了能够有效地计算高频域的声学问题，国内外学者在近年来提出了许多改进型的数值算法，例如伽辽金/最小二乘有限元法[4-5]、高阶有限元算法[6-7]，无网格法[8]和光滑有限元法[9-13]等，其中直接“软化”(降低模型的计算刚度)离散模型的光滑有限元法能够明显地提高标准有限元法的求解精度[13]。在常规有限元模型中，基于伽辽金法得到的标准有限元模型的刚度比实际的连续系统的刚度大，所得到的波数比实际波数要小，从而导致了数值离散误差。为了解决这一问题，刘桂荣[14]提出了梯度光滑技术并将之运用于节点光滑的点插值法(Node-based smoothed point interpolation method,NS-PIM)和节点光滑有限元法(Node based smoothed finite element method,NS-FEM)中。但是这两种方法又表现得“过软”(计算刚度比实际刚度小)，为此刘桂荣等人提出了α有限元法[15](Alpha finite element method,α-FEM)。这种方法通过使用参数α来结合“过刚”的标准有限元法和“过软”的节点光滑有限元法，从而得到较为准确的数值模型刚度矩阵；因此相比于标准有限元法，α有限元法能够在高频率的声散射计算中得到更准确的计算结果。

在二维声学问题的有限元法中，随着波数k的增高，标准有限元法的计算误差也随之增大，甚至会出现完全不可信的计算结果；针对此问题，本文根据文献[14-15]，把α有限元法推广到无界域的声学计算中，运用声压梯度光滑技术，结合标准有限元和光滑有限元的特点，在高频下得到更加准确的计算结果。在本文中，分别对标准有限元法和α有限元法对求解声散射问题进行了推导，采用了DtN无反射边界条件[16]来将无限域转化为有限域。通过两个典型的数值算例，将α有限元解、标准有限元解和解析解或参考解来对比，验证了α有限元法的高精度和高计算效率特性。

1 二维声学Helmholtz方程和DtN人工边界

理想介质中小振幅声波波动方程为：

考虑到小振幅声波的时谐特性，上式可表示为以下的Helmholtz方程：

式中：k表示波数。

图1 水下声散射问题计算模型示意Fig.1 The numerical model of the underwater acoustic scattering problems

声学计算中，通常考虑散射体的三种不同典型边界条件：即Dirichlet边界条件ΓDΓD，对应声学软边界；Neumann边界条件ΓNΓN，对应声学硬边界；Robin边界条件ΓRΓR，对应声学阻抗边界。本文算例中散射体为刚性体，则其边界条件为Neumann边界条件ΓNΓN，而问题(计算)域采用无反射边界条件是DtN边界条件，具体边界条件的表达式如下：

ΓN：p·n=-jρωvn

DtN：

式中：j表示复数单位，υnvn表示相应的边界上质点的法向振动速度，n表示相应边界上的单位法向量，AnAn表示导纳系数，ρ表示介质密度，M表示DtN算子，ω表示角频率。

在二维情况下，DtN算子的表达式为[16]：

采用伽辽金加权残值法，将式(2)所示的Helmholtz方程两边同时乘上一个加权函数w并且带入式(3)中所示的声学边界条件，同时运用分部积分和格林公式在整个问题域内进行积分，就得到以下声学问题控制方程的伽辽金“弱”形式：

式中：Γ为声学问题域的边界。

2 α有限元离散系统方程

在标准有限元中：声压p=NpP=∑NiPi=NP，Ni表示与节点i相关联的单元形函数，Pi表示节点i处的声压，将其代入式(5)可得：

其矩阵形式为：

式中

PT=[P1,P2,…,Pn] 节点声压矩阵

图2 基于节点的光滑域及其背景网格Fig.2 Node-based smoothing domain and background mesh

在节点光滑有限元中，所用的背景网格和标准有限元法是相同的，如图2所示：网格中共包含Ne个单元和Nn个节点。在二维情况下，通过连接一个节点的所有相邻三角单元的形心和三角形单元边的中点，形成与该节点相关联的光滑域Ωk；于是，网格中总共包含有Nn个光滑域。在NS-FEM中，使用的单元形函数与标准有限元相同，不同之处是N被光滑算子替代由节点梯度光滑技术算得[17-18]。NS-FEM中的声学刚度矩阵表达式为：

以上积分的计算均在基于节点的光滑域内，故而有如下形式：

式中：K(k)是与节点k有关的相关联的光滑域Ωk的单元刚度矩阵，计算式为：

式中：Ak为光滑域Ωk的面积，在二维情况下：

式中：Ni表示标准有限元中节点i的形函数，Γk表示光滑域的边界，np为积分域边界上的单位法向量。

图3 α-FEM单元离散示意Fig.3 Illustration of domain discretization for α-FEM

除此以外，本文计算的是外声场问题，采用了DtN边界条件，根据Givoli 和 Kaller[16]的推导，在DtN边界条件下的标准有限元计算中，刚度矩阵由两个部分组成：

式中：Kb是人工边界矩阵，它包含DtN算子M和有限元形函数：

将式(13)～(15)代入式(7)中，得到最终的离散系统方程矩阵形式如下：

应当指出的是：与标准FEM 相比，α-FEM既没有增加数学模型的单元数目和节点数目，也没有增加系统总的自由度数目；但α-FEM中的系统总体刚度矩阵的“带宽”却要大于FEM中系统总体刚度矩阵的“带宽”，也就是说α-FEM模型得到的刚度矩阵中的“非零”元素数目与标准FEM模型中的刚度矩阵中的“非零”元素数目相比要大；这是因为在装配标准FEM的单元刚度矩阵时，仅仅需要单个单元内部的节点信息，而在装配α-FEM的单元刚度矩阵时，除了单元内部的节点信息外，其相邻单元的节点信息也要参与单元刚度矩阵的装配过程，也就是说α-FEM的单元刚度矩阵考虑了更多的节点信息的影响；因此，α-FEM能比标准FEM得到更加精确的计算结果是不难理解的。

为了考察α-FEM求解声学问题的计算精度和计算效率，参考文献[1]的研究成果，本文定义了声学数值计算的误差因子，它能够较为客观地反映数值方法求解声学问题的求解精度，其表达式为：

式中：v表示声场中的质点振动速度，符号“～”表示相关变量的共轭复数，上标e和上撇号分别表示解析法和数值方法得到的计算结果。从式中可以看出，该误差因子能够较为客观地反映声学问题数值解的全局误差(也即求解精度)，该误差因子的值越小，则表明所用数值方法的求解精度越高；反之则越低。

3 数值算例

3.1无限长刚性圆柱在平面波入射下的声散射计算

以水中无限长刚性圆柱对平面波的散射为例进行计算分析。无限长刚性圆柱的半径为0.2 m，人工边界半径为1 m，有限元网格尺寸0.02 m，共有3 806个节点，7 356个单元。水中声速v=1 500 m/s，水密度ρ=1 000 kg/m3。无限长刚性圆柱在平面波入射下的声散射解析解为：

首先需确定α的值，本文中是通过计算对比来找到确定的α；图4所示的是在波数k=30时人工边界上的解析解和α-FEM解。α-FEM解的三条曲线分别对应三个α值：0.7、0.8、0.9。从图中可以看到，当α=0.8时数值解与解析解最为接近，故而后续计算均使用α=0.8。

图4 α分别为0.7、0.8、0.9时的α-FEM解和解析解在人工边界上的声压值(Pa)Fig.4 The acoustic pressure results (Pa) on the artificial boundary from analytical solution and α-FEM with different values of α: α=0.7，α=0.8，α=0.9

图5～图8显示的是α=0.8时，α-FEM解与FEM解在不同波数下的比较。从图中可以看出，当k=5时，两者的结果均非常接近解析解，当k=15、25、30时，α-FEM解仍然比较接近解析解，而FEM解误差逐渐增大，特别是在前向散射处最为明显。

图5 波数k=5时解析解、FEM解和α-FEM解在人工边界上声压值(Pa)Fig.5 The comparison of the acoustic pressure results (Pa) on the artificial boundary from the analytical solution，FEM and α-FEM for wavenumber k=5

图6 波数k=15时解析解、FEM解和α-FEM解在人工边界上声压值(Pa)Fig.6 The comparison of the acoustic pressure results (Pa) on the artificial boundary from the analytical solution，FEM and α-FEM for wavenumber k=15

图7 波数k=25时解析解、FEM解和α-FEM解在人工边界上声压值(Pa)Fig.7 The comparison of the acoustic pressure results (Pa) on the artificial boundary from the analytical solution，FEM and α-FEM for wavenumber k=25

图8 波数k=30时解析解、FEM解和α-FEM解在人工边界上声压值(Pa)Fig.8 The comparison of the acoustic pressure results (Pa) on the artificial boundary from the analytical solution，FEM and α-FEM for wavenumber k=30

图9给出的是α-FEM解、FEM解和解析解在人工边界上特定观察点的声压随波数的变化曲线。从图中可以看出，在低波数情况下，三者结果十分接近，随着波数的增高，α-FEM解和FEM解均会逐渐与解析解产生偏差，但α-FEM解的偏差远小于FEM解的偏差。

图9 人工边界上θ=0°，θ=30°，θ=45°，θ=60°四个观察点的α-FEM解、FEM解和解析解随波数变化曲线Fig.9 The acoustic pressure results from α-FEM， FEM and analytical method versus the increasing wavenumber at different observation points θ=0°，θ=30°，θ=45°，θ=60°

图10 FEM和α-FEM的求解时间结果比较Fig.10 The CPU time versus the numerical error indicator for both FEM and α-FEM

从上节内容可知，在同样的网格情况下，α-FEM的总体刚度矩阵“带宽”是要大于标准FEM的，这一现象导致的结果就是：在相同的条件下，运用α-FEM求解声学问题所需要的时间要比标准FEM要长，也就是说虽然α-FEM的求解精度提高了，但是其求解时间也变长了；为了进一步考察α-FEM的计算效率，即比较α-FEM和标准FEM在相同的求解精度前提下所需要的时间，运用了4套由疏到密的不同网格来计算上述的声散射问题；图10所示的就是α-FEM和标准FEM的计算效率结果比较，从图中可以看出：在相同的网格情况下，α-FEM的求解时间确实要比标准FEM的求解时间要长，但是在相同计算效率的前提下，α-FEM所需的求解时间明显要小于标准FEM，也就是说α-FEM的计算效率是要高于标准FEM的。

3.2刚性舵截面在平面波入射下的声散射计算

刚性舵截面的尺寸如图11，长2.1 m，宽0.4 m，对比用的α-FEM和FEM的计算网格尺寸是0.1 m，共有4 692个节点，9 108个单元。由于没有解析解，故将更细网格划分的FEM解作为参考解，其网格尺寸为0.05 m，共有12 795个节点，25 125个单元。

图12～图14对应的波数分别为k=0.5π、π和1.5π时，α-FEM解、FEM解和参考解在人工边界上的声压值。从图中可以看出，其结果与无限长刚性圆柱的情况类似，在低波数下α-FEM解、FEM解和参考解的结果十分接近，而当波数增高时，α-FEM解和FEM解均产生偏差，但FEM解的偏差更大，在前向散射范围内尤为明显。

图11 刚性舵截面尺寸Fig.11 Geometry of the rudder scatterer

图12 波数k=0.5π时α-FEM、FEM和参考解在人工边界上的声压值(Pa)Fig.12 Comparison of the acoustic pressure results (Pa) on the artificial boundary from α-FEM，FEM and reference solutions for the rudder at wavenumber k=0.5π

图13 波数k=π时α-FEM、FEM和参考解在人工边界上的声压值(Pa)Fig.13 Comparison of the acoustic pressure results (Pa) on the artificial boundary from α-FEM，FEM and reference solutions for the rudder at wavenumber k=π

图14 波数k=1.5π时α-FEM、FEM和参考解在人工边界上的声压值(Pa)Fig.14 Comparison of the acoustic pressure results (Pa) on the artificial boundary from α-FEM，FEM and reference solutions for the rudder at wavenumber k=1.5π

4 结 语

本文结合了α-FEM和DtN人工边界，计算了无界域中无限长刚性圆柱和刚性舵截面在平面波入射下的散射声场，通过对α-FEM解、FEM解、解析解和参考解进行对比和分析，得到以下几点结论：

1)参数α决定了光滑有限元刚度和标准有限元刚度所占的比例，影响α-FEM计算结果的准确性。

2)α-FEM与FEM使用的是相同的网格，所以α-FEM模型可以从FEM模型略微修改得到，带来计算上的方便。

3)随着波数的增高，标准有限元的误差逐渐增大，而α有限元的结果仍保有较高的计算精度。

4)通过使用梯度光滑技术，α-FEM“软化”系统刚度矩阵，减小了离散误差。数值计算结果表明：与标准有限元法相比，α有限元法在水下声散射计算中具备更高的计算精度和计算效率。

[1] IHLENBURG F,BABUKA I.Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation[J].International Journal for Numerical Methods in Engineering,1995,38(22): 3 745-3 774.

[3] DERAEMAEKER A,BABUKA I,BOUILLARD P.Dispersion and pollution of the FEM solution for the Helmholtz equation in one,two and three dimensions[J].International Journal for Numerical Methods in Engineering,1999,46(4): 471-499.

[4] HARARI I,HUGHES T J R.Galerkin/least-squares finite element methods for the reduced wave equation with non-reflecting boundary conditions in unbounded domains[J].Computer Methods in Applied Mechanics and Engineering,1992,98(3): 411-454.

[5] THOMPSON L L,PINSKY P M.A Galerkin least squares finite element method for the two-dimensional Helmholtz equation[J].International Journal for Numerical Methods in Engineering,1995 38(3): 371-397.

[6] PETERSEN S,DREYER D,VON ESTORFF O.Assessment of finite and spectral element shape function or efficient 1iterative simulations of interior acoustics[J].Computer Methods in Applied Mechanics and Engineering,2006,195(44): 6 463-6 478.

[7] BIERMANN J,VON ESTORFF O,PETERSEN S,et al.Higher order finite and infinite elements for the solution of Helmholtz problems[J].Computer Methods in Applied Mechanics and Engineering,2009,198(13): 1 171-1 188.

[8] BOUILLARD P,SULEAU S.Element-free Garlekin solutiong for Helmholtz problems: formulation and numerical assessment of the pollution effect[J].Computer Methods in Applied Mechanics and Engineering,1998,162(1): 317-335.

[9] LI W,CHAI Y B,LEI M,et al.Analysis of coupled structural-acoustic problems based on the smoothed finite element method(S-FEM)[J].Engineering Analysis with Boundary Elements,2014,42: 84-91.

[10] CHAI Y B,LI W,GONG Z X,et al.Hybrid smoothed finite element method for two dimensional acoustic radiation problems[J].Applied Acoustics,2016,103:90-101.

[11] CHAI Y B,LI W,LI TY,et al.Analysis of underwater acoustic scattering problems using stable node-based smoothed finite element method[J].Engineering Analysis with Boundary Elements,2016,72: 27-41.

[12] CHAI Y B,LI W,GONG Z X,et al.Hybrid smoothed finite element method for two-dimensional underwater acoustic scattering problems[J].Ocean Engineering,2016,116:129-141.

[13] LI W,CHAI Y B,LEI M,et al.Numerical investigation of the edge-based gradient smoothing technique for exterior Helmholtz equation in two dimensions[J].Computers and Structures,2017,182: 149-164.

[14] LIU G R.A generalized gradient smoothing technique and the smoothed bilinear form for Galerkin formulation of wide class of computational methods[J].International Journal of Computational Methods,2008,5(2): 199-236.

[15] LIU G R,NGUYEN T T,LAM K Y.A novel FEM by scaling the gradient of strains with factor α(αFEM)[J].Computational Mechanics,2009,43(3): 369-391.

[16] KELLER J B,GIVOLI D.Exact non-reflecting boundary conditions[J].Journal of Computational Physics,1989,82(1): 172-192.

[17] LIU G R,ZHANG G Y,DAI K Y,et al.A linearly conforming point interpolation method for 2D solid mechanics problems[J].International Journal of Computational Methods,2005,2(4): 645-665.

[18] ZHANG G Y,LIU G R,WANG Y Y,et al.A linearly conforming point interpolation method for three-dimensional elasticity problems[J].International Journal for Numerical Methods in Engineering,2007,72(13): 1 524-1 543.

Alpha finite element method for two-dimensional underwater acoustic scattering numerical computation

JIN Wenchao1,YUE Zhijun2,LI Wei3,LI Yaofei3,CHAI Yingbin3

(1.Naval Academy of Armament,Beijing,100036,China; 2.China Ship Development and Design Center,Wuhan,430064,China; 3.School of Naval Architecture and Ocean Engineering,Huazhong University of Science and Technology,Wuhan,430074,China)

It is well known that the finite element method (FEM) is unreliable to solve acoustic problems at high wavenumbers due to the interpolation and pollution error.In order to improve the performance of the standard FEM solutions for acoustic problems,a novel alpha finite element method (α-FEM) is presented to analyze the underwater acoustic scattering problems.In the α-FEM，the node-based gradient smoothing technique is used to obtain the smoothed node-based gradient field and then the node-based smoothed finite element (NS-FEM) is formulated.The α-FEM contains both the gradient components from the NS-FEM and the gradient components from the standard FEM，and takes advantage of the properties from the “overly-soft” NS-FEM and the “overly-stiff” FEM.In order to handle the exterior acoustic problems，the Dirichle-to-Neumann map is introduced.Several typical numerical examples are presented and it is verified that the α-FEM possesses higher computational efficiency and can provide more accurate solutions than standard FEM for underwater acoustic scattering problems.

alpha finite element method (α-FEM); underwater acoustic scattering; unbounded domain; non-reflecting boundary

P733.24；O427.2

A

10.16483/j.issn.1005-9865.2017.05.017

1005-9865(2017)05-0141-08

2016-10-04

国家自然科学基金资助项目(51579112)

晋文超(1985- )，男，安徽合肥人，工程师，从事舰船工程研究。E-mail: 812114744@qq.com

李 威。E-mail: hustliw@hust.edu.cn