吳紹維，柯 磊，韓國文

(1. 重慶交通大學 航運與船舶工程學院，重慶 400074；2. 船舶動力工程技術交通運輸行業(yè)重點實驗室，武漢 430063)

水下聲輻射問題是近年來非?；钴S的研究領域，基于網(wǎng)格的方法是現(xiàn)階段計算水下聲輻射的主要手段，如邊界元法(boundary element method, BEM)和有限元法(finite element method, FEM)。BEM和FEM均可直接求解內部聲學問題，對于外部聲問題，BEM在無窮遠處滿足Sommerfeld輻射條件，是求解無限域聲學問題的有效手段[1]。但BEM在處理大型問題時，由于相關系數(shù)矩陣為非對稱和非稀疏矩陣，導致計算效率變低[2]。FEM是求解大型問題非常有效的方法[3]，然而標準的FEM不能直接用于求解無限域中的外部聲學問題，需對問題域進行處理，常采用人工邊界截斷無限域獲取有限計算域，為確保聲波向外傳播過程中的自由衰減[4]，需在人工邊界處施加無反射邊界條件，如DtN(Dirichlet-to-Neumann)邊界條件、完美匹配層、無限元和吸收邊界條件等。在標準的FEM中，存在數(shù)值色散誤差[5-6]，細化網(wǎng)格不能有效降低該誤差[7]。近年來提出了一些基于光滑數(shù)值和無網(wǎng)格的方法來軟化有限元模型“過硬”的剛度[8-10]。借助這些方法，色散誤差得以減小。

用無窮傅里葉級數(shù)表示的DtN邊界條件是一種精確的非局部邊界條件[11-12]，解析構建了Dirichlet變量與Neumann變量之間的關系。在計算時，需將無窮級數(shù)截斷為有限項數(shù)N，當N不足夠大時，無法確保解的唯一性和可解性[13-14]。因此，截斷形式的DtN邊界條件不再精確，雖然通過增加項數(shù)可解決該問題，但對高頻或大尺度人工邊界計算量過大。為克服該缺陷，MDtN(modified Dirichlet-to-Neumann)邊界條件被提出，對于任意的N值，可確保聲場求解的唯一性和可解性[15]。與截斷的DtN邊界條件相比，MDtN邊界條件使用少量的項數(shù)可獲得更高的精度。此外，在聲場計算中無法預先確定具體的項數(shù)以達到期望的精度。因此，MDtN邊界條件更具優(yōu)勢。

與基于網(wǎng)格的方法相比，新興的無網(wǎng)格法無需使用節(jié)點連接信息或網(wǎng)格進行變量插值或近似[16]，已被廣泛應用于處理一些具有挑戰(zhàn)性的數(shù)值和工程問題。無網(wǎng)格點插值法是一種全局弱式法[17]，該方法采用多項式點插值法(polynomial point interpolation method, PPIM)構造形函數(shù)，使之具有Kronecker-delta函數(shù)性質，可以像有FEM一樣直接施加本質邊界條件。然而，PPIM形成的矩矩陣可能是奇異的。使用兩階矩陣三角化算法可消除該奇異問題[18]，但由于PPIM形函數(shù)的不相容性，使得該方法對不規(guī)則分布節(jié)點的魯棒性較差，且局部支持域中含有過多節(jié)點會形成過高階的插值多項式，導致PPIM形函數(shù)劇烈振蕩。為消除奇異性，基于徑向基函數(shù)(radial basis function，RBF)的徑向點插值法(radial point interpolation method, RPIM)被提出[19]，其對任意分布的節(jié)點具有良好的穩(wěn)定性和魯棒性，具有比FEM更快的收斂速度。

RPIM已被用于內部聲學問題，能獲得比FEM更精確的結果[20]。作為一種基于域離散的方法，RPIM需與其他數(shù)值方法聯(lián)合處理無限域聲問題[21]。在水下聲輻射計算領域，為抑制中頻數(shù)值色散誤差，基于光滑有限元耦合DtN邊界的方法被提出，包括混合光滑有限元法、穩(wěn)定點基光滑有限元法[22]及邊基梯度光滑法[23]。這些耦合方法有效提高了水下中頻段聲輻射和聲散射計算精度。在光滑有限元框架下，Xu等[24]提出了單元基光滑徑向基點插值耦合DtN邊界條件方法計算水下聲輻射，將RPIM引入水下噪聲計算，與傳統(tǒng)的FEM相比，該方法在計算精度、收斂速度和抑制數(shù)值色散誤差方面具有優(yōu)勢?；赗PIM和MDtN邊界條件的各自優(yōu)點，本文提出了一種基于RPIM和MDtN邊界條件相耦合的無網(wǎng)格弱式方法來精確預報水下聲輻射。研究了無網(wǎng)格插值中的尺寸參數(shù)和形參數(shù)對聲場預報精度的影響規(guī)律，確定了可用取值范圍。結果表明，與FEM相比，所提方法在計算精度、收斂速度以及計算效率方面具有優(yōu)勢，對聲波數(shù)的敏感度顯著降低。

1 水下聲輻射計算理論背景

1.1 控制方程

在無限大自由場V中，穩(wěn)態(tài)簡諧結構振動輻射聲波的聲壓p(x)滿足如下Helmholtz方程

?2p(x)+k2p(x)=0,x∈V

(1)

式中： ?2為拉普拉斯算子；k=ω/c為聲波波數(shù)，ω為結構振動圓頻率，c為流體介質中的聲傳播速度。結構邊界可分為Dirichlet型邊界ΓD和Neumann型邊界ΓN，ΓD和ΓN滿足ΓD∪ΓN=Γ及ΓD∩ΓN=?，其定義為

p=pD, onΓD

(2)

?p(x)·nb=-ikρcv(x), onΓN

(3)

式中：p=pD為ΓD上的已知聲壓；v(x)為ΓN上給定的法向振速；nb為邊界外法向；ρ為流體介質密度； i2=-1。

對exp(iωt)形式的簡諧時間因子，為確保聲波向外傳播過程中的自由衰減和聲場求解的唯一性要求，需滿足如下Sommerfeld輻射條件[25]

(4)

式中：r為場點到坐標原點的距離；β為維度參數(shù)，對二維和三維問題分別取1和2。若p(x,t)=p(x)exp(-iωt)，則?p/?r-ikp=o(r-β)。通常只用第二項表示Sommerfeld 輻射條件，在無窮遠處，此輻射條件與ρc阻抗條件等效，即?p·n∞+ikp=0。

1.2 MDtN邊界條件

為獲取有限計算域，采用人工邊界ΓR(對2D和3D問題分別為半徑為R的圓和球)截斷無限問題域，有限計算域以Γ為內邊界，ΓR為外邊界，如圖1所示。在人工邊界上施加DtN邊界條件?p(x)/?n=-Mp(x)，式中M為DtN映射算子，其作用是在ΓR上建立Dirichlet量聲壓p與Neumann量?p/?n的關系。從而，最初的聲場求解問題可等效為如下兩個聲場求解問題。

(5)

(6)

圖1 使用一組離散分布場點表示的有限計算域Fig.1 Finite computational domain represented by a set of discrete distributed field nodes

為完整起見，簡要介紹DtN邊界條件理論，更為詳細的理論和推導可參考Keller等的研究。對二維聲問題，采用極坐標(r,θ)可解析求解出式(6)對應問題的聲壓為

(7)

(8)

其中，

(9)

mn(θ-θ′)=

(10)

在實際計算中，需將式(8)中的無窮級數(shù)截斷為N個有限項，則式(8)變?yōu)?/p>

(11)

式中，MN為截斷后的DtN映射算子。從而截斷后的DtN邊界條件不再精確，當N取值較小時，無法確保問題的可解性和數(shù)值解的唯一性。為克服該缺陷并提高計算精度，MDtN邊界條件被提出，以確保N取任意值時具有唯一的解。通過在式(8)右邊加上和減去Bp(R,θ)，同時將Bp(R,θ)求和并截斷，得到如下MDtN邊界條件

(12)

(13)

則MDtN邊界條件可確定為

(14)

其中，

(15)

2 RPIM-MDtN計算理論

2.1 形函數(shù)構造

采用ph(x)表示聲壓p(x)的試函數(shù)，W(x)表示檢驗函數(shù)，式(1)兩邊乘以W(x)，并在VI上積分得到

(16)

根據(jù)格林第一公式，得到如下方程

(17)

式中：n為包裹VI的邊界外法向；在Γ上，n=-nb。將式(17)代入式(16)得到

(18)

將式(3)和式(12)代入式(18)得到如下弱式方程

(19)

用散布在VI及其邊界上的一組離散場點表示該有限計算域及其邊界，VI中任意計算點x處的聲壓可使用場點聲壓值插值確定。

(20)

式中：Ri(x)為RBF；pj(x)為基于空間坐標x=[x,y,z]T表示的單項式；n和m分別為RBF基函數(shù)和多項式的項數(shù)；a=[a1,a2,…,an]T和b=[b1,b2, …,bm]T為未知的常數(shù)系數(shù)向量；RT=[R1(x),R2(x),…,Rn(x)]，pT=[p1(x),p2(x), …,pm(x)]。

為求解式(20)中的未知系數(shù)a和b，通過在計算點x的支持域Ωs中的n個場點處滿足式(20)來構造n個方程，并結合如下m個附加方程

(21)

得到如下系統(tǒng)方程

(22)

式中：ps=[p1,p2, …,pn]T為聲壓值向量；

(23)

(24)

對任意分布的場點，R-1通常存在，且式(20)通常使用低階多項式，則RPIM中不會出現(xiàn)奇異問題。求解式(22)得到系數(shù)a和b，將其代入式(20)得到

(25)

式中，NT(x)=[N1(x),N2(x), …,Nn(x)]為n個場點的RPIM形函數(shù)。在RPIM形函數(shù)的構造中，需確定局部支持域，通常采用影響域來確定計算點x的支持域中包含的場點，若計算點位于某場點的影響域中，則該場點參與構造該計算點的形函數(shù)。在后續(xù)部分，為了簡化起見，采用圓形影響域。根據(jù)建議，采用rw=αsd確定影響域半徑rw，式中αs表示場點xi的影響域的無量綱尺寸，該參數(shù)對計算精度至關重要，對于給定類型的問題，通常采用數(shù)值試驗來確定αs取值。

2.2 RPIM-MDtN離散系統(tǒng)方程

在加權殘差公式中，采用RPIM構造的形函數(shù)N(x)作為檢驗函數(shù)，則式(19)可表示為

(26)

經(jīng)離散化處理，建立如下矩陣形式的離散系統(tǒng)方程

[K-k2M+Kb]{p}=F

(27)

(28)

(29)

(30)

(31)

式中： {p}為聲壓值的向量；Cn(x)=[cos(nθ),sin(nθ)]。以上這些矩陣是基于場點i和j的全局索引組裝，計算這些系統(tǒng)矩陣需用到一組全局背景單元，其不依賴于場點，不參與變量插值，對插值精度無影響。只有當場點i和j同時屬于至少一個積分點的支持域時，Kij≠0和Mij≠0，否則，Kij=0且Mij=0。因此，全局系數(shù)矩陣K和M為稀疏矩陣。

3 數(shù)值算例

采用數(shù)值算例對所提方法進行驗證。通過使用具有解析解的浸沒無限長圓柱數(shù)值試驗確定用于水下聲輻射計算的影響域無量綱尺寸和形狀參數(shù)的取值，并對該方法的精度、收斂性以及計算效率開展研究。進一步通過浸沒無限長圓柱局部聲輻射模型、舵形結構和二維潛艇結構的聲輻射算例，對所提方法進行檢驗。

3.1 浸沒周向簡諧振動無限長圓柱聲輻射

無限長圓柱半徑為a=0.25 m，聲介質為水，密度取ρ=1 000 kg/m3，聲速為c=1 500 m/s，以坐標原點為中心構建R=0.5 m圓形人工邊界，如圖2所示。在圓柱的濕表面上(r=a)指定周向Dirichlet邊界條件p(θ)=cos(nθ)，對應的歸一化解析解為

(32)

采用εreal和εimag分別表示聲壓實部和虛部數(shù)值均方根誤差

(33)

(34)

圖2 浸沒周向簡諧振動無限長圓柱聲輻射模型Fig.2 Submerged infinitely long cylinder with circumferential oscillation

3.1.1 影響域無量綱尺寸對精度的影響

對一類基準問題，通常需進行數(shù)值試驗來預先確定αs的可用取值，本節(jié)通過分析聲壓預報誤差隨αs的變化來研究αs對精度的影響規(guī)律。使用一組平均節(jié)點間距為d=0.023 m的不規(guī)則分布場點離散有限計算域，采用496個二次四邊形單元作為全局背景單元實施積分運算(如圖3所示)，構造形函數(shù)時在式(20)添加線性多項式。聲波波數(shù)滿足ka=π/2，使得a等于一個波長。為了確保無網(wǎng)格全局弱式法的積分精度，對于二維問題，積分點nq與場點nf的數(shù)量應滿足nq>2nf/3。由于nf=1 171，為精確計算系統(tǒng)矩陣系數(shù)，在背景單元的每個方向上使用5個Gauss-Legendre積分點。圖4和圖5分別為第二階(n=1)和第三階(n=2)輻射模式下人工邊界處的εreal和εimag隨αs的變化曲線，其中ac=0.4和q=1.03(MQ)，ac=2(EXP)，η=3(TPS)，N=7(MDtN)。

圖3 使用任意分布的場點表示的計算域及用于積分的背景單元Fig.3 Representing computational domain using arbitrarily distributed nodes and performing integration with the use of background cells

圖4 第二階輻射模式下εreal和εimag隨αs的變化Fig.4 εreal and εimag against αs for second circumferential mode

圖5 第三階輻射模式下εreal和εimag隨αs的變化Fig.5 εreal and εimag against αs for third circumferential mode

由圖中曲線可見，聲壓預報誤差隨αs改變而變化，αs取值對所提方法的精度具有顯著的影響。為了獲得精確的計算結果，需合理選取αs。當影響域尺寸過小時(αs≤1.5)，矩陣G為奇異矩陣，導致計算失敗。因此，圖中沒有給出對應αs≤1.5的計算結果。采用EXP-RBF作為基函數(shù)時，用所提方法獲得的結果在2≤αs≤6范圍可接受；當αs<2或αs>6時，精度降低。對于MQ-RBF和TPS-RBF，在αs≥2范圍內，所提方法可以產(chǎn)生精確的結果，但隨著αs的增大，并不能有效提高計算精度，其原因可能是αs與形狀參數(shù)不匹配，對于該問題還需作進一步研究。然而，較大的影響域將顯著增加計算時間，尤其是對于大規(guī)模問題。因此，需保持計算精度和計算耗時之間的平衡。結果還表明，相比EXP-RBF，采用MQ-TPF和TPS-RBF作為基函數(shù)具有更高的精度。為了不增加計算耗時并獲得較高的計算精度，在后續(xù)算例中，對EXP-RBF取αs=3，對MQ-RBF和TPS-RBF取αs=4.5。

3.1.2 形狀參數(shù)對精度的影響

為研究形狀參數(shù)對所提方法精度的影響，圖6～圖9給出了人工邊界處的εreal和εimag隨形參數(shù)的變化。由圖中曲線可知，計算精度受形狀參數(shù)的影響，需合理選取形狀參數(shù)的值以獲得精確的計算結果。當采用MQ-RBF時，由圖6可知，αc可取值范圍為0≤αc≤5；圖7中的曲線表明，形狀參數(shù)q對計算精度具有顯著影響，q=j(j=1,2,3,…)時，RPIM中的G矩陣為奇異矩陣，導致錯誤的計算結果。對于EXP-RBF，由圖8可見，αc在0.07～2.50內可產(chǎn)生較精確的計算結果。當使用TPS-RBF時，從圖9給出的結果可見，η應在1≤η≤9范圍取值，需注意，當η取值接近2j(j=1,2,3,…)時，矩陣G的條件數(shù)變大，導致計算精度急劇變差。根據(jù)以上結果，在后續(xù)研究中，基函數(shù)的形狀參數(shù)取值為：對MQ-RBF，取αc=0.5和q=1.1；對EXP-RBF，取αc=1.9；對TPS-RBF，取η=3.5。

圖6 對應MQ-RBF的εreal和εimag隨αc的變化Fig.6 εreal and εimag against αc corresponding to MQ-RBF

圖7 對應MQ-RBF的εreal和εimag隨q的變化Fig.7 εreal and εimag against q corresponding to MQ-RBF

圖8 對應EXP-RBF的εreal和εimag隨αc的變化Fig.8 εreal and εimag against αc corresponding to EXP-RBF

圖9 對應TPS-RBF的εreal和εimag隨η的變化Fig.9 εreal and εimag against q corresponding to TPS-RBF

3.1.3 有限項數(shù)N的影響

為研究MDtN邊界條件中的截斷項數(shù)對所提方法精度的影響規(guī)律，圖10給出了在第五階輻射模式下人工邊界處的誤差隨有限項數(shù)N的變化，同時還給出了使用DtN邊界條件產(chǎn)生的數(shù)值誤差，其中波數(shù)k滿足ka=2π。結果表明：與截斷的DtN邊界條件相比，使用MDtN邊界條件可獲得更高的計算精度，尤其是對于較小的N取值，計算精度提升明顯。對于截斷的DtN邊界條件，通過設定N≥kR，雖能夠克服適定性問題，但由于通常無法預先確定所需項數(shù)來達到期望的精度，計算結果仍可能嚴重偏離真實值。與截斷的DtN邊界條件不同，MDtN邊界條件對積分算子內核中包含的低階模態(tài)是精確的，而局部算子包含了與高階模態(tài)對應的近似非反射邊界條件。因此，MDtN邊界條件對應的誤差在N=5時快速收斂到穩(wěn)定值，在后續(xù)算例中，取N=6。

圖10 MDtN和DtN邊界條件對應的誤差隨N的變化Fig.10 Errors obtained from MDtN and DtN maps versus N

3.1.4 收斂性和計算效率研究

為研究收斂性，采用四種平均場點間距(d=0.011, 0.020, 0.032,0.043)進行數(shù)值試驗。圖11給出在第三階輻射模式下所提方法和FEM的收斂曲線，其中k滿足ka=2π，兩種方法采用相同的場點配置，使用二次有限元單元和和對應的節(jié)點。圖11中的曲線表明，在相同的配置下，與FEM相比，所提出的方法具有更高的收斂率，能夠產(chǎn)生更精確的結果；當使用MQ-EBF或EXP-RBF時，收斂過程不穩(wěn)定，一個可能的原因是無量綱尺寸與形狀參數(shù)不匹配。到目前為止，還沒有可用的方法從理論上確定無量綱尺寸和形狀參數(shù)的取值，只能通過數(shù)值試驗獲得，需進一步從理論上確定其最優(yōu)值。

在相同配置下，通過對比所提方法與有限元法的計算耗時來開展計算效率研究，所有方案及程序均在相同的計算環(huán)境中運行(Intel?Xeon?W-2255 CPU 3.7 GHz和RAM 64 GB)。在效率研究中，采用TPS-RBF作為基函數(shù)。圖12(a)給出了在ka=2π時，CPU耗時隨d的變化。結果表明：在相同配置下，所提方法需更多的計算時間。其原因是，對所提方法，必須為所有積分點支持域內的場點構造形函數(shù)，而在FEM中，已預先確定覆蓋積分點的單元的FEM形函數(shù)；此外，與FEM相比，所提方法需使用更多的場點信息來組裝系統(tǒng)矩陣，導致系統(tǒng)矩陣的帶寬較FEM大。眾所周知，有效的數(shù)值方法應以短的計算耗時獲取高計算精度，因此將計算耗時與精度結合能合理評計算效率，圖12(b)給出了CPU計算時間隨聲壓實部誤差的變化。由結果可見，當εreal=17.5%時，該方法的計算耗時約為FEM的39%，從而當要求高精度時，所提方法具有計算效率優(yōu)勢。

圖11 所提方法與FEM的收斂性比較Fig.11 Comparison of convergence curves between the proposed scheme and FEM scheme

圖12 所提方法與FEM的計算效率比較Fig.12 Comparison of computational efficiency between the proposed scheme and FEM scheme

3.2 浸沒無限長圓柱局部聲輻射

為進一步檢驗采用所提方法預報聲場的準確性，采用浸沒無限長圓柱的局部聲輻射模式進行數(shù)值試驗，對于該算例，在圓柱邊界局部范圍(-φ≤θ≤φ,φ=π/6)指定p=1，在其余邊界處p=0，如圖13所示。幾何尺寸、相關參數(shù)和平均場點間距與上述算例相同，對應的歸一化解析解為

(35)

圖14和圖15分別為k=4π和k=8π，聲壓沿人工邊界的分布。由結果可見，數(shù)值解與解析解總體上吻合良好，高聲波數(shù)對應的聲壓值與解析解略有偏離，但計算結果仍可接受，與EXP-RBF相比，使用MQ-RBF和TPS-RBF作為基函數(shù)可獲得更精確的結果。

圖13 浸沒無限長圓柱局部聲輻射Fig.13 Non-uniform acoustic radiation by a sector of cylinder

圖14 當k=4π時人工邊界處的聲壓分布Fig.14 Sound pressure on artificial boundary when k=4π

圖15 當k=8π時人工邊界處的聲壓分布Fig.15 Sound pressure on artificial boundary when k=8π

3.3 舵形結構聲輻射

浸沒于水中的剛性舵形結構的形狀和尺寸如圖16所示，水介質密度ρ=1 000 kg/m3，聲速c=1 500 m/s，在結構的濕表面施加vn=10-5m/s的Neumann邊界條件，以半徑為R=4 m的圓作為人工邊界獲取有限計算域，并采用1 193個雙線性四節(jié)點單元對其進行離散(含1 277個節(jié)點)。為進行對比，在相同場點配置下，給出FEM計算結果。由于該問題不存在解析解，采用細化網(wǎng)格(68 677個四節(jié)點單元，69 357個節(jié)點)對應的FEM結果作為參考值。圖17和圖18和分別給出kL=4π和kL=8π時聲壓沿人工邊界的分布，其中RPIM采用TPS-RBF作為基函數(shù)進行插值。由圖中結果可見，在給定頻率處，RPIM-MDtN獲得的聲壓分布與參考值吻合良好。對于較小的波數(shù)，F(xiàn)EM-MDtN計算結果與參考值具有良好的一致性，但隨著聲波數(shù)的增加，在kL=8π時，F(xiàn)EM預報的聲壓嚴重偏離參考值。結果表明，所提方法對聲波數(shù)敏感度低，在水下聲輻射計算中具有更高的精度。

圖16 剛性舵形結構的幾何形狀與尺寸Fig.16 Shape and geometry of rigid rudder-shaped structure

圖17 當kL=4π時人工邊界處的聲壓分布Fig.17 Sound pressure on artificial boundary when kL=4π

圖18 當kL=8π時人工邊界處的聲壓分布Fig.18 Sound pressure on artificial boundary when kL=8π

3.4 潛艇結構聲輻射

本節(jié)采用浸沒于水中的二維潛艇結構的聲輻射試驗來檢驗所提方法的實用性，水的密度ρ=1 000 kg/m3，聲速c=1 500 m/s，潛艇形狀和尺寸如圖19所示。對于潛艇，需評估其輻射噪聲水平，敵人聲納通常距離被探測物體較遠，在所提RPIM-MDtN方法中，在確定人工邊界處的聲壓值后利用式(7)可計算任意遠場點處的聲壓。在實際預報潛艇噪聲時，潛艇的艉部和艏部船體的振動模態(tài)取決于分析頻率[26]，為簡單起見，僅考慮位于潛艇尾部的動力裝置引起的船尾振動。在此算例中，在尾部施加如圖19所示的振動速度為vn=10-5m/s的Neumann邊界條件。以原點為中心，半徑為R=50 m的圓作為人工邊界截斷無限問題域，采用10 545個雙線性四節(jié)點單元(10 841個節(jié)點)離散計算域，并在人工邊界處選取觀測點(θ=π/2)。由于該聲輻射問題不存在解析解，使用細化的網(wǎng)格(177 820個四節(jié)點單元，179 005個節(jié)點)對應的FEM結果作為參考值。

圖19 二維海底形結構的幾何形狀和尺寸Fig.19 Size and shape of submarine-shaped structure

圖20～圖21給出了f=100 Hz及400 Hz時，聲壓沿人工邊界的分布，其中采用TPS-RBF構造形函數(shù)。為進行對比，在相同配置下給出FEM的計算結果(10 545個雙線性四節(jié)點單元，10 841個節(jié)點)。由圖中結果可見，在低頻f=100 Hz時，RPIM-MDtN和FEM-MDtN方案均可產(chǎn)生良好的結果，隨著頻率的增加，RPIM-MDtN預報的聲壓依然非常接近參考值，而FEM-MDtN預報的聲壓在f=400 Hz時明顯偏離了參考值。為進一步檢驗所提方法的精度，研究觀測點處的10～400 Hz的聲壓頻率響應，其中頻率間隔為1.0 Hz。RPIM-MDtN和FEM-MDtN預報的聲壓及參考值繪于圖22。由結果可見，相比FEM-MDtN，RPIM-MDtN在全頻率范圍內具有更高的精度且在高頻處預報的結果與參考值仍然具有良好的一致性，而FEM-MDtN在高頻處預報的結果嚴重偏離參考值。結果再次表明：與FEM相比，所提方法在水下聲輻射預報中能夠獲得更高的計算精度，對聲波數(shù)的敏感度顯著降低，能夠有效抑制中頻色散誤差。

圖20 當f=100 Hz時人工邊界上的聲壓分布Fig.20 Sound pressure on artificial boundary at f=100 Hz

圖21 當f=400 Hz時人工邊界上的聲壓分布Fig.21 Sound pressure on artificial boundary at f=400 Hz

圖22 觀測點處的聲壓頻率響應Fig.22 Sound pressure frequency response at observation point

4 結 論

基于徑向點插值法和修正的Dirichlet-to-Neumann邊界條件提出了一種預報水下結構聲輻射的無網(wǎng)格弱式方法，通過算例對該方法進行了研究，主要結論總結如下：

(1) RPIM-MDtN能夠精確且穩(wěn)定地計算無限域水下聲輻射，無需使用網(wǎng)格或節(jié)點連接信息來進行場變量插值。

(2) 與傳統(tǒng)的FEM相比，RPIM-MDtN能夠獲得更精確的聲場預報結果，在收斂性方面表現(xiàn)得更好，對聲波數(shù)的敏感度顯著降低，能夠有效抑制中頻色散誤差。

(3) 由于PRIM形狀函數(shù)構造，在相同場場點配置下，RPIM-MDtN方法比標準的FEM需要更多的計算時間，當要求高計算精度時，所提方法具有計算效率優(yōu)勢。

(4) 相比EXP-RBF，MQ-RBF和TPS-RBF具有更高的精度；TPS-RBF具有穩(wěn)定的收斂性，而EXP-RBF和MQ-RBF的收斂過程不穩(wěn)定，對于該問題還需作進一步的研究。

(5) 在所提方法中，影響域無量綱尺寸和形參數(shù)取值均采用數(shù)值試驗獲取，可能存在無量綱尺寸與形狀參數(shù)不匹配的問題，需從理論上進一步研究其優(yōu)化取值。