一種預(yù)報(bào)水下聲輻射的無(wú)網(wǎng)格弱式徑向點(diǎn)插值方法

吳紹維，柯 磊，韓國(guó)文

(1. 重慶交通大學(xué) 航運(yùn)與船舶工程學(xué)院，重慶 400074；2. 船舶動(dòng)力工程技術(shù)交通運(yùn)輸行業(yè)重點(diǎn)實(shí)驗(yàn)室，武漢 430063)

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

用無(wú)窮傅里葉級(jí)數(shù)表示的DtN邊界條件是一種精確的非局部邊界條件[11-12]，解析構(gòu)建了Dirichlet變量與Neumann變量之間的關(guān)系。在計(jì)算時(shí)，需將無(wú)窮級(jí)數(shù)截?cái)酁橛邢揄?xiàng)數(shù)N，當(dāng)N不足夠大時(shí)，無(wú)法確保解的唯一性和可解性[13-14]。因此，截?cái)嘈问降腄tN邊界條件不再精確，雖然通過(guò)增加項(xiàng)數(shù)可解決該問(wèn)題，但對(duì)高頻或大尺度人工邊界計(jì)算量過(guò)大。為克服該缺陷，MDtN(modified Dirichlet-to-Neumann)邊界條件被提出，對(duì)于任意的N值，可確保聲場(chǎng)求解的唯一性和可解性[15]。與截?cái)嗟腄tN邊界條件相比，MDtN邊界條件使用少量的項(xiàng)數(shù)可獲得更高的精度。此外，在聲場(chǎng)計(jì)算中無(wú)法預(yù)先確定具體的項(xiàng)數(shù)以達(dá)到期望的精度。因此，MDtN邊界條件更具優(yōu)勢(shì)。

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

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

1 水下聲輻射計(jì)算理論背景

1.1 控制方程

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

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

(1)

式中： ?2為拉普拉斯算子；k=ω/c為聲波波數(shù)，ω為結(jié)構(gòu)振動(dòng)圓頻率，c為流體介質(zhì)中的聲傳播速度。結(jié)構(gòu)邊界可分為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為邊界外法向；ρ為流體介質(zhì)密度； i2=-1。

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

(4)

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

1.2 MDtN邊界條件

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

(5)

(6)

圖1 使用一組離散分布場(chǎng)點(diǎn)表示的有限計(jì)算域Fig.1 Finite computational domain represented by a set of discrete distributed field nodes

為完整起見(jiàn)，簡(jiǎn)要介紹DtN邊界條件理論，更為詳細(xì)的理論和推導(dǎo)可參考Keller等的研究。對(duì)二維聲問(wèn)題，采用極坐標(biāo)(r,θ)可解析求解出式(6)對(duì)應(yīng)問(wèn)題的聲壓為

(7)

(8)

其中，

(9)

mn(θ-θ′)=

(10)

在實(shí)際計(jì)算中，需將式(8)中的無(wú)窮級(jí)數(shù)截?cái)酁镹個(gè)有限項(xiàng)，則式(8)變?yōu)?/p>

(11)

式中，MN為截?cái)嗪蟮腄tN映射算子。從而截?cái)嗪蟮腄tN邊界條件不再精確，當(dāng)N取值較小時(shí)，無(wú)法確保問(wèn)題的可解性和數(shù)值解的唯一性。為克服該缺陷并提高計(jì)算精度，MDtN邊界條件被提出，以確保N取任意值時(shí)具有唯一的解。通過(guò)在式(8)右邊加上和減去Bp(R,θ)，同時(shí)將Bp(R,θ)求和并截?cái)?，得到如下MDtN邊界條件

(12)

(13)

則MDtN邊界條件可確定為

(14)

其中，

(15)

2 RPIM-MDtN計(jì)算理論

2.1 形函數(shù)構(gòu)造

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

(16)

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

(17)

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

(18)

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

(19)

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

(20)

式中：Ri(x)為RBF；pj(x)為基于空間坐標(biāo)x=[x,y,z]T表示的單項(xiàng)式；n和m分別為RBF基函數(shù)和多項(xiàng)式的項(xiàng)數(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，通過(guò)在計(jì)算點(diǎn)x的支持域Ωs中的n個(gè)場(chǎng)點(diǎn)處滿足式(20)來(lái)構(gòu)造n個(gè)方程，并結(jié)合如下m個(gè)附加方程

(21)

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

(22)

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

(23)

(24)

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

(25)

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

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

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

(26)

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

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

(27)

(28)

(29)

(30)

(31)

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

3 數(shù)值算例

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

3.1 浸沒(méi)周向簡(jiǎn)諧振動(dòng)無(wú)限長(zhǎng)圓柱聲輻射

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

(32)

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

(33)

(34)

圖2 浸沒(méi)周向簡(jiǎn)諧振動(dòng)無(wú)限長(zhǎng)圓柱聲輻射模型Fig.2 Submerged infinitely long cylinder with circumferential oscillation

3.1.1 影響域無(wú)量綱尺寸對(duì)精度的影響

對(duì)一類基準(zhǔn)問(wèn)題，通常需進(jìn)行數(shù)值試驗(yàn)來(lái)預(yù)先確定αs的可用取值，本節(jié)通過(guò)分析聲壓預(yù)報(bào)誤差隨αs的變化來(lái)研究αs對(duì)精度的影響規(guī)律。使用一組平均節(jié)點(diǎn)間距為d=0.023 m的不規(guī)則分布場(chǎng)點(diǎn)離散有限計(jì)算域，采用496個(gè)二次四邊形單元作為全局背景單元實(shí)施積分運(yùn)算(如圖3所示)，構(gòu)造形函數(shù)時(shí)在式(20)添加線性多項(xiàng)式。聲波波數(shù)滿足ka=π/2，使得a等于一個(gè)波長(zhǎng)。為了確保無(wú)網(wǎng)格全局弱式法的積分精度，對(duì)于二維問(wèn)題，積分點(diǎn)nq與場(chǎng)點(diǎn)nf的數(shù)量應(yīng)滿足nq>2nf/3。由于nf=1 171，為精確計(jì)算系統(tǒng)矩陣系數(shù)，在背景單元的每個(gè)方向上使用5個(gè)Gauss-Legendre積分點(diǎn)。圖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 使用任意分布的場(chǎng)點(diǎn)表示的計(jì)算域及用于積分的背景單元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

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

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

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

圖6 對(duì)應(yīng)MQ-RBF的εreal和εimag隨αc的變化Fig.6 εreal and εimag against αc corresponding to MQ-RBF

圖7 對(duì)應(yīng)MQ-RBF的εreal和εimag隨q的變化Fig.7 εreal and εimag against q corresponding to MQ-RBF

圖8 對(duì)應(yīng)EXP-RBF的εreal和εimag隨αc的變化Fig.8 εreal and εimag against αc corresponding to EXP-RBF

圖9 對(duì)應(yīng)TPS-RBF的εreal和εimag隨η的變化Fig.9 εreal and εimag against q corresponding to TPS-RBF

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

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

圖10 MDtN和DtN邊界條件對(duì)應(yīng)的誤差隨N的變化Fig.10 Errors obtained from MDtN and DtN maps versus N

3.1.4 收斂性和計(jì)算效率研究

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

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

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

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

3.2 浸沒(méi)無(wú)限長(zhǎng)圓柱局部聲輻射

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

(35)

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

圖13 浸沒(méi)無(wú)限長(zhǎng)圓柱局部聲輻射Fig.13 Non-uniform acoustic radiation by a sector of cylinder

圖14 當(dāng)k=4π時(shí)人工邊界處的聲壓分布Fig.14 Sound pressure on artificial boundary when k=4π

圖15 當(dāng)k=8π時(shí)人工邊界處的聲壓分布Fig.15 Sound pressure on artificial boundary when k=8π

3.3 舵形結(jié)構(gòu)聲輻射

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

圖16 剛性舵形結(jié)構(gòu)的幾何形狀與尺寸Fig.16 Shape and geometry of rigid rudder-shaped structure

圖17 當(dāng)kL=4π時(shí)人工邊界處的聲壓分布Fig.17 Sound pressure on artificial boundary when kL=4π

圖18 當(dāng)kL=8π時(shí)人工邊界處的聲壓分布Fig.18 Sound pressure on artificial boundary when kL=8π

3.4 潛艇結(jié)構(gòu)聲輻射

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

圖19 二維海底形結(jié)構(gòu)的幾何形狀和尺寸Fig.19 Size and shape of submarine-shaped structure

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

圖20 當(dāng)f=100 Hz時(shí)人工邊界上的聲壓分布Fig.20 Sound pressure on artificial boundary at f=100 Hz

圖21 當(dāng)f=400 Hz時(shí)人工邊界上的聲壓分布Fig.21 Sound pressure on artificial boundary at f=400 Hz

圖22 觀測(cè)點(diǎn)處的聲壓頻率響應(yīng)Fig.22 Sound pressure frequency response at observation point

4 結(jié) 論

基于徑向點(diǎn)插值法和修正的Dirichlet-to-Neumann邊界條件提出了一種預(yù)報(bào)水下結(jié)構(gòu)聲輻射的無(wú)網(wǎng)格弱式方法，通過(guò)算例對(duì)該方法進(jìn)行了研究，主要結(jié)論總結(jié)如下：

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

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

(3) 由于PRIM形狀函數(shù)構(gòu)造，在相同場(chǎng)場(chǎng)點(diǎn)配置下，RPIM-MDtN方法比標(biāo)準(zhǔn)的FEM需要更多的計(jì)算時(shí)間，當(dāng)要求高計(jì)算精度時(shí)，所提方法具有計(jì)算效率優(yōu)勢(shì)。

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

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