基于α有限元法的二維水下聲散射計算

晉文超，岳智君，李 威，李耀飛，柴應彬

(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