聲學(xué)邊界元擬奇異積分計(jì)算的自適應(yīng)方法

繆宇躍, 李天勻,3, 朱 翔, 張冠軍, 郭文杰

(1.華中科技大學(xué) 船舶與海洋工程學(xué)院, 武漢 430074; 2.船舶與海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室, 武漢 430074;3.船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室, 武漢 430064)

聲學(xué)邊界元擬奇異積分計(jì)算的自適應(yīng)方法

繆宇躍1,2, 李天勻1,2,3, 朱 翔1,2, 張冠軍1,2, 郭文杰1,2

(1.華中科技大學(xué) 船舶與海洋工程學(xué)院, 武漢 430074; 2.船舶與海洋水動(dòng)力湖北省重點(diǎn)實(shí)驗(yàn)室, 武漢 430074;3.船舶振動(dòng)噪聲重點(diǎn)實(shí)驗(yàn)室, 武漢 430064)

提出一種自適應(yīng)方法計(jì)算聲學(xué)邊界元中的擬奇異積分，通過(guò)單元分級(jí)細(xì)分將總積分轉(zhuǎn)移到子單元上以消除擬奇異性。在此方法基礎(chǔ)上深入研究擬奇異性，進(jìn)一步提出接近度的概念，其中臨界接近度可作為擬奇異積分計(jì)算的理論依據(jù)，并可用于預(yù)估擬奇異性是否存在。此方法的積分精度可調(diào)控，且不受場(chǎng)點(diǎn)位置限制，相比于已有方法更加靈活高效。數(shù)值分析表明擬奇異性強(qiáng)弱由場(chǎng)點(diǎn)與單元的相對(duì)位置決定,單元上遠(yuǎn)離場(chǎng)點(diǎn)的區(qū)域擬奇異性很弱，無(wú)需處理。研究結(jié)果為處理邊界元法中的擬奇異性問(wèn)題提供了新的選擇和參考。

邊界元；擬奇異性；自適應(yīng)；接近度

聲學(xué)邊界元法經(jīng)過(guò)幾十年的發(fā)展, 已經(jīng)成為求解結(jié)構(gòu)振動(dòng)聲輻射問(wèn)題的重要方法，得到了比有限元法更廣泛的應(yīng)用，尤其與快速多極算法結(jié)合后[1]，體現(xiàn)出較強(qiáng)的實(shí)用性。奇異積分問(wèn)題[2,3]在邊界元法中由來(lái)已久，已有多種處理方法，隨后擬奇異積分問(wèn)題也成為研究熱點(diǎn)。擬奇異性，又稱(chēng)為近奇異性[4]或幾乎奇異性[5]，是由于場(chǎng)點(diǎn)十分靠近結(jié)構(gòu)表面，導(dǎo)致場(chǎng)點(diǎn)與源點(diǎn)間的距離很微小。雖然其物理機(jī)理不存在奇異性，但在數(shù)值計(jì)算時(shí)，這種微距使常規(guī)的高斯積分法失效，產(chǎn)生類(lèi)似奇異性的收斂性差問(wèn)題。

為處理邊界元法的擬奇異性，國(guó)內(nèi)外研究者提出了多種方法。TANAKA等[6]提供了一個(gè)較詳實(shí)的關(guān)于擬奇異積分正則化方法的綜述文獻(xiàn)。ZHANG等[7]提出了一種半解析法，采用曲面單元，有效地處理了薄壁結(jié)構(gòu)的擬奇異性。針對(duì)二維彈性問(wèn)題，NIU等[8]提出一種解析的擬奇異積分算法，隨后這種方法被應(yīng)用到二維的正交各向異性勢(shì)問(wèn)題[9]和各向異性勢(shì)問(wèn)題[10]中。對(duì)于三維勢(shì)問(wèn)題，NIU等[11]將奇異性降階，并結(jié)合正則化方法能夠更容易地計(jì)算擬奇異積分。

自適應(yīng)方法源自自動(dòng)控制優(yōu)化領(lǐng)域，并在有限元法中得到發(fā)展，隨后這種自適應(yīng)的思想被引入邊界元法的研究中，但是主要被用來(lái)進(jìn)行整體或局部的單元細(xì)分以提高計(jì)算精度[12]，KITA等[13]對(duì)自適應(yīng)邊界元法的發(fā)展應(yīng)用做了較詳實(shí)的概括。根據(jù)自適應(yīng)單元細(xì)分的特性，研究者開(kāi)始將這種方法應(yīng)用于擬奇異積分的計(jì)算中。WEI等[14]基于CFD理論，結(jié)合有限元和邊界元法分析了水下航行體噪聲，其中邊界元的奇異性問(wèn)題采用自適應(yīng)單元平均細(xì)分法處理。CROAKER等[15]對(duì)氣流中結(jié)構(gòu)的誘導(dǎo)聲散射進(jìn)行預(yù)報(bào)研究，為了得到近場(chǎng)壓力和壓力梯度，對(duì)擬奇異性體積分和面積分采用了單元細(xì)分法。DONG等[16]提出非穩(wěn)態(tài)導(dǎo)熱問(wèn)題的邊界元算法，為避免小時(shí)間步長(zhǎng)的奇異性，采用坐標(biāo)變換結(jié)合體單元細(xì)分法來(lái)提高域積分精度。DAVIES等[17]提出三維彈塑性邊界元法研究復(fù)雜的地質(zhì)力學(xué)問(wèn)題，將域積分變換為邊界積分消除強(qiáng)奇異性，并運(yùn)用自適應(yīng)細(xì)分法優(yōu)化積分過(guò)程，其中單元細(xì)分過(guò)程依據(jù)經(jīng)驗(yàn)公式進(jìn)行。GAO等[18]進(jìn)一步將自適應(yīng)單元細(xì)分法和解析法相結(jié)合計(jì)算多種類(lèi)型的二維奇異積分并取得良好效果。LI等[19-22]同樣依據(jù)這種經(jīng)驗(yàn)公式進(jìn)行單元細(xì)分來(lái)處理三維問(wèn)題的擬奇異性，其中子單元尺寸以及子單元與源點(diǎn)距離需要通過(guò)反復(fù)迭代計(jì)算和誤差分析確定，而這種經(jīng)驗(yàn)公式未充分考慮場(chǎng)點(diǎn)與單元相對(duì)位置的不同所存在的差異。ZHANG等[23-24]利用邊界點(diǎn)法和邊界面法研究三維勢(shì)問(wèn)題，擬奇異積分通過(guò)自適應(yīng)單元細(xì)分實(shí)現(xiàn)，其細(xì)分準(zhǔn)則是：比較每個(gè)子單元的對(duì)角線長(zhǎng)度和近場(chǎng)點(diǎn)到該子單元中心的距離，若前者小于后者則認(rèn)為細(xì)分達(dá)到積分要求，否則繼續(xù)進(jìn)行細(xì)分以達(dá)到積分要求。依據(jù)這種“對(duì)角線準(zhǔn)則”的自適應(yīng)細(xì)分法能夠準(zhǔn)確計(jì)算擬奇異積分，但若運(yùn)用在聲學(xué)分析中則過(guò)于苛刻，增加不必要的計(jì)算量。

本文提出分級(jí)細(xì)分自適應(yīng)方法，制定合理易行的細(xì)分方案，將多種情況下的擬奇異積分計(jì)算統(tǒng)一起來(lái)。在深入研究場(chǎng)點(diǎn)與單元不同的相對(duì)位置所存在差異的基礎(chǔ)上提出接近度的概念，為擬奇異積分計(jì)算提供理論基礎(chǔ)，并預(yù)估擬奇異性程度以確定是否需要處理擬奇異性。數(shù)值分析證明了這種自適應(yīng)方法能夠很好的計(jì)算擬奇異積分，具有較好的靈活性和高效性，研究結(jié)果可為自適應(yīng)邊界元法的發(fā)展提供思路和參考。

1 聲學(xué)邊界元法基本理論

對(duì)于自由空間的聲輻射問(wèn)題，其Helmholtz邊界積分方程為：

C(z)p(z)=

(1)

將基本解中的指數(shù)項(xiàng)按Taylor級(jí)數(shù)展開(kāi)：

(2)

并設(shè):

(3)

則有:

(4)

(5)

(6)

可見(jiàn)基本解G的弱奇異性是由Gr引起的，Gr是與頻率無(wú)關(guān)的，而Gk是無(wú)奇異性的，可以進(jìn)行高斯積分。

場(chǎng)點(diǎn)在S上時(shí)，對(duì)于平面單元，r⊥ns，故?r/?ns=0，則?G/?ns=?G/?r·?r/?ns=0，只存在弱奇異性。場(chǎng)點(diǎn)不在S上時(shí)，r⊥ns這個(gè)條件不成立，故?G/?ns≠0，則同時(shí)存在弱奇異性和強(qiáng)奇異性。實(shí)際上弱奇異積分比強(qiáng)奇異積分更容易收斂，所以只需保證強(qiáng)奇異積分收斂，便能消除積分方程的奇異性。

2 自適應(yīng)細(xì)分方法

本文提出分級(jí)細(xì)分的自適應(yīng)方法，并與其他自適應(yīng)方法對(duì)比驗(yàn)證精度和收斂性。由于二維高斯積分是通過(guò)等參變換建立各種不規(guī)則單元與規(guī)則四邊形單元的映射關(guān)系，之后在一個(gè)局部坐標(biāo)系的正方形平面區(qū)域數(shù)值求解積分，所以文中統(tǒng)一采用頂點(diǎn)坐標(biāo)為[-1, -1]、[1, -1]、[1, 1]和[-1, 1]的標(biāo)準(zhǔn)正方形單元來(lái)說(shuō)明自適應(yīng)細(xì)分方法，如圖1所示。下文中所有省略的長(zhǎng)度單位都是米。圖1中黑點(diǎn)表示近場(chǎng)點(diǎn)在單元上的投影點(diǎn)，每行從左往右是自適應(yīng)細(xì)分方案，第(1)行是平均細(xì)分方案，文獻(xiàn)[14]在處理奇異積分時(shí)便是采取這種方案；第(2)～(5)行是投影點(diǎn)在不同位置時(shí)的分級(jí)細(xì)分方案。

平均細(xì)分方法是將原始單元平均細(xì)分為四個(gè)子單元，在下一次細(xì)分過(guò)程中每個(gè)子單元再平均細(xì)分為四個(gè)更小的子單元，如此遞推下去將原始單元上的積分轉(zhuǎn)移到各個(gè)子單元上去，直到全部子單元積分總和的誤差達(dá)到收斂容許值為止。

分級(jí)細(xì)分方法是圍繞近場(chǎng)點(diǎn)在單元上的投影點(diǎn)來(lái)細(xì)分，對(duì)原始單元進(jìn)行平均細(xì)分后，只將投影點(diǎn)所在的子單元再次平均細(xì)分，其他子單元不變，如此下去直到收斂。由于投影點(diǎn)位置不同，于是產(chǎn)生了(2)～(5)這幾種分級(jí)細(xì)分方案。

圖1 自適應(yīng)細(xì)分方案Fig.1. The adaptive subdivision processes

細(xì)分后，原始單元上的弱奇異積分和強(qiáng)奇異積分就變?yōu)樽訂卧系姆e分總和：

(7)

(8)

式中:s為原始單元，sj為子單元，N為子單元總數(shù)，Ng為高斯點(diǎn)數(shù)，Jj為子單元Jacobian，wu和wv為權(quán)系數(shù)。

設(shè)num為細(xì)分次數(shù)和子單元級(jí)別，則圖1中五種細(xì)分方案的子單元總數(shù)為

(9)

由于平面單元的X、Y和Z坐標(biāo)都是線性變化的，以X坐標(biāo)為例說(shuō)明如何根據(jù)自適應(yīng)算法得到子單元坐標(biāo)。對(duì)于平均細(xì)分方案，如圖1中的第(1)行， 第一級(jí)子單元與原始單元的坐標(biāo)關(guān)系如下：

(10)

其中T1～T4是坐標(biāo)傳遞矩陣，X0是原始單元坐標(biāo)：

(11)

于是通過(guò)同樣的坐標(biāo)傳遞，第二級(jí)子單元坐標(biāo)為

(12)

如此繼續(xù)下去，細(xì)分num次獲得的所有子單元坐標(biāo)為

(13)

式中:m1=1～4，m2=1～4。

對(duì)于分級(jí)細(xì)分方案，如圖1中第(2)行，只對(duì)投影點(diǎn)所在單元進(jìn)行平均細(xì)分，每次細(xì)分產(chǎn)生三個(gè)不含投影點(diǎn)的子單元，最后剩下一個(gè)投影點(diǎn)所在單元，故第num級(jí)子單元坐標(biāo)為

(14)

圖1中第(3)～(5)行的分級(jí)細(xì)分方案具有很高相似性，其子單元坐標(biāo)的自適應(yīng)算法都可參照第(2)行。由此便可清晰地構(gòu)建各種情況下的自適應(yīng)細(xì)分方法，不需要迭代和判斷，只需要設(shè)定細(xì)分次數(shù)。

3 近場(chǎng)點(diǎn)的接近度

在進(jìn)行收斂性分析之前，先要提出接近度的概念，并根據(jù)接近度決定單元?jiǎng)澐值拇螖?shù)。由于近場(chǎng)點(diǎn)與單元的相對(duì)位置并不固定，不同的接近程度對(duì)擬奇異性的影響是有差別的。文獻(xiàn)[25]初步指出擬奇異性是近場(chǎng)點(diǎn)與高斯積分點(diǎn)相互作用的結(jié)果，接近高斯點(diǎn)的區(qū)域會(huì)產(chǎn)生較強(qiáng)的奇異性，其他區(qū)域則很弱。為了深入發(fā)掘接近度對(duì)擬奇異性的影響，本文在標(biāo)準(zhǔn)單元周?chē)×巳舾山鼒?chǎng)點(diǎn)進(jìn)行研究，如圖2所示。由于正方形的對(duì)稱(chēng)性，只在八分之一直角三角形OAB附近取十個(gè)近場(chǎng)點(diǎn)，采用十點(diǎn)高斯積分，①～⑥號(hào)點(diǎn)投影在OAB的邊上，⑦～⑩號(hào)點(diǎn)在平面OAB上，坐標(biāo)如表1所示，其中t為場(chǎng)點(diǎn)到單元的距離。使每一近場(chǎng)點(diǎn)沿直線不斷接近單元并考察收斂性：若頻率無(wú)關(guān)的強(qiáng)奇異積分∫s-1/r2·?r/?nsds在三次細(xì)分以?xún)?nèi)收斂并且不細(xì)分的積分值與細(xì)分收斂后的積分值相對(duì)誤差不超過(guò)1%，那么此時(shí)的距離稱(chēng)為“近遠(yuǎn)場(chǎng)”的臨界值，超過(guò)這個(gè)臨界距離就不存在擬奇異性，場(chǎng)點(diǎn)不再是近場(chǎng)點(diǎn)，單元也不需要細(xì)分就可正常積分。本文定義這種場(chǎng)點(diǎn)到單元的距離t與單元邊長(zhǎng)l的比值稱(chēng)為近場(chǎng)點(diǎn)的接近度E，即E=t/l。通過(guò)逼近計(jì)算得到場(chǎng)點(diǎn)到單元的臨界距離t0和臨界接近度E0如表2所示，其中⑦～⑩號(hào)點(diǎn)的?r/?ns=0，故這些點(diǎn)的強(qiáng)奇異積分以∫s-1/r2ds表示。

圖2 近場(chǎng)點(diǎn)與標(biāo)準(zhǔn)單元的相對(duì)位置Fig.2 The relative positions of near-field points to the standard element

從表2看出①～⑥號(hào)點(diǎn)和⑦～⑩號(hào)點(diǎn)的臨界距離t0和臨界接近度E0各不相同，明顯呈現(xiàn)出隨著與O點(diǎn)距離增大而減小的趨勢(shì)，表明場(chǎng)點(diǎn)離高斯積分區(qū)域越遠(yuǎn)，相互作用越弱，擬奇異性程度越低。

表1 近場(chǎng)點(diǎn)坐標(biāo)

表2 臨界距離和臨界接近度

文獻(xiàn)[17-22]采用的經(jīng)驗(yàn)公式為

(15)

式中:A1為近場(chǎng)點(diǎn)到積分單元最小距離，A2為單元在積分方向上的長(zhǎng)度，A3為該積分方向上的高斯點(diǎn)數(shù)，A4為奇異性階次，A5為該積分方向上的高斯積分誤差。對(duì)于標(biāo)準(zhǔn)單元取A2=2，A3=10，A4=2，A5=0.01，計(jì)算得到臨界距離為A1=0.168 6，并不能滿足任意位置近場(chǎng)點(diǎn)的誤差要求，未體現(xiàn)出相對(duì)位置變化所帶來(lái)的差異。

由于分級(jí)細(xì)分方法是圍繞近場(chǎng)點(diǎn)在單元上的投影點(diǎn)來(lái)細(xì)分，實(shí)際造成場(chǎng)點(diǎn)與最小子單元的位置關(guān)系等同于⑤號(hào)點(diǎn)與標(biāo)準(zhǔn)單元的位置關(guān)系，所以采用十點(diǎn)高斯積分并且誤差閾值為1%時(shí)，分級(jí)細(xì)分次數(shù)num可根據(jù)接近度(場(chǎng)點(diǎn)到單元距離與最小子單元邊長(zhǎng)之比)來(lái)決定，即t/(l/2num)≥0.04，于是得到：

(16)

式中l(wèi)為原始單元平均邊長(zhǎng)。式(16)可作為擬奇異積分的單元細(xì)分準(zhǔn)則。

4 收斂性分析

首先以標(biāo)準(zhǔn)單元(l=2)的收斂性分析驗(yàn)證分級(jí)自適應(yīng)方法和單元細(xì)分準(zhǔn)則的正確性。采用十點(diǎn)高斯積分，k=1 rad/m，取四個(gè)場(chǎng)點(diǎn)(編號(hào)1～4)，坐標(biāo)分別為[0, 0, 0.01]、[1, 1, 0.01]、[0, 0, 0.001]和[1, 1, 0.001]，四個(gè)場(chǎng)點(diǎn)對(duì)應(yīng)的強(qiáng)奇異積分分別表示為IR1、IR2、IR3和IR4，每次細(xì)分的積分值與最終收斂值的誤差分別表示為er1、er2、er3和er4，收斂性曲線如圖3所示。

圖3 標(biāo)準(zhǔn)單元的積分收斂性曲線Fig.3 The convergence curves for the standard element

圖3中只展示了IR實(shí)部，因?yàn)镮R虛部幾乎沒(méi)有變化，可見(jiàn)強(qiáng)奇異性由積分的實(shí)部來(lái)體現(xiàn)，而弱奇異積分收斂非?？?，同樣不予展示。通過(guò)誤差曲線er1和er2看出經(jīng)過(guò)3次以上細(xì)分，誤差就能降到1%以下，而誤差曲線er3和er4顯示經(jīng)過(guò)7次以上細(xì)分，誤差也降到1%以下。可見(jiàn)對(duì)于t分別為0.01和0.001的場(chǎng)點(diǎn)，收斂情況完全符合式(16)的細(xì)分準(zhǔn)則。

若按照式(15)計(jì)算得到離場(chǎng)點(diǎn)最近子單元長(zhǎng)度分別為0.118 6和0.011 9，需要細(xì)分5次和8次；若按照文獻(xiàn)[23-24]中的“對(duì)角線準(zhǔn)則”計(jì)算得到離場(chǎng)點(diǎn)最近子單元長(zhǎng)度分別為0.008 2和0.000 82，需要細(xì)分8次和12次，都增加了額外的計(jì)算量。另外，式(15)和“對(duì)角線準(zhǔn)則”需要考察每個(gè)子單元是否滿足要求，單元的細(xì)分是一種反復(fù)迭代判斷過(guò)程，而根據(jù)本文的單元細(xì)分依據(jù)式(16)可以進(jìn)行預(yù)判，然后按本文的自適應(yīng)方法一次性劃分子單元和計(jì)算坐標(biāo)，更加簡(jiǎn)捷實(shí)用。

再以寬為2，長(zhǎng)分別4和8的兩種矩形單元(分別簡(jiǎn)稱(chēng)為“矩四”和“矩八”)進(jìn)行驗(yàn)證。令t=0.01，近場(chǎng)點(diǎn)5和6分別投影于“矩四”的中心和頂點(diǎn)；近場(chǎng)點(diǎn)7和8分別投影于“矩八”的中心和頂點(diǎn)。相應(yīng)的奇異積分分別表示為IR5～8，每次細(xì)分的積分值與最終收斂值的誤差分別表示為er5～8，收斂性曲線如圖4所示。

圖4 矩形單元的積分收斂性曲線Fig.4 The convergence curves for the rectangular elements

圖4顯示對(duì)于t為0.01的場(chǎng)點(diǎn)，分別經(jīng)過(guò)4次和5次以上細(xì)分，誤差也降到1%以下，收斂情況同樣完全符合式(16)的細(xì)分要求。若按照式(15)的要求，需要細(xì)分6次和9次；若按照“對(duì)角線準(zhǔn)則”，需要細(xì)分11次和13次，這樣也增加了額外計(jì)算量。

值得注意的是，采用分級(jí)細(xì)分方法和平均細(xì)分方法得到的計(jì)算結(jié)果并無(wú)差別，這充分說(shuō)明了分級(jí)細(xì)分方法的準(zhǔn)確性，同時(shí)也表明單元上不同區(qū)域?qū)M奇異性的影響不同：越靠近場(chǎng)點(diǎn)的區(qū)域?qū)ζ娈愋缘呢暙I(xiàn)越大，而遠(yuǎn)離節(jié)點(diǎn)的區(qū)域?qū)ζ娈愋缘呢暙I(xiàn)迅速減小。常規(guī)高斯積分未考慮這種差別導(dǎo)致積分不準(zhǔn)，平均細(xì)分方法平均化處理整個(gè)單元會(huì)增加很多不必要得計(jì)算量，而分級(jí)細(xì)分方法梯度化地對(duì)待遠(yuǎn)近區(qū)域，符合奇異性產(chǎn)生的物理本質(zhì)，既保證了精度也節(jié)省了計(jì)算量。根據(jù)式(9)，兩種細(xì)分方法產(chǎn)生的子單元總數(shù)如表3所示。

表3 兩種細(xì)分方法產(chǎn)生的子單元總數(shù)

由以上積分收斂性分析可見(jiàn)，自適應(yīng)分級(jí)細(xì)分方法能克服弱奇異性和強(qiáng)奇異性，有效計(jì)算擬奇異積分。同時(shí)依據(jù)本文所提的細(xì)分準(zhǔn)則能夠準(zhǔn)確預(yù)判，不需要反復(fù)迭代判斷，產(chǎn)生子單元數(shù)量更少。

另外，對(duì)于邊長(zhǎng)比例不協(xié)調(diào)、形狀不規(guī)則的單元，以及場(chǎng)點(diǎn)靠近單元邊角的情況，將遠(yuǎn)離場(chǎng)點(diǎn)的區(qū)域做細(xì)分是不必要的。為了很快地收斂，只需將投影點(diǎn)所在區(qū)域劃分出來(lái)，得到一個(gè)比例較為規(guī)則的小塊(陰影部分)，然后按照?qǐng)D1中分級(jí)細(xì)分方案進(jìn)行細(xì)分，小塊以外區(qū)域保持不變，如圖5所示。

圖5 特殊情況下的細(xì)分方法Fig.5 Some especial cases for subdivision

5 數(shù)值算例

5.1 預(yù)估擬奇異性

對(duì)于任意不規(guī)則單元，可用表2中臨界接近度預(yù)估此種情況下是否存在擬奇異性。單元頂點(diǎn)和近場(chǎng)點(diǎn)的投影點(diǎn)坐標(biāo)如圖6所示，投影點(diǎn)落在不規(guī)則單元的一頭。

圖6 任意不規(guī)則單元和近場(chǎng)點(diǎn)的投影點(diǎn)Fig.6 The arbitrary irregular element and the subpoint

首先根據(jù)表2中②號(hào)點(diǎn)的臨界接近度E0作判斷：

t/l≥0.130

(17)

得到臨界距離為0.333 0，其中l(wèi)取為長(zhǎng)短邊的平均值，那么當(dāng)近場(chǎng)點(diǎn)到單元距離t不小于0.333 0時(shí)，一定不存在擬奇異性。當(dāng)t小于0.333 0時(shí)，將投影點(diǎn)所在部分劃分出來(lái)(虛線將Y軸右半部分單元平分)，再根據(jù)表2中①號(hào)點(diǎn)的臨界接近度E0作判斷：

t/l≤0.150 5

(18)

得到臨界距離為0.152 8，其中l(wèi)取為投影點(diǎn)所在子單元的長(zhǎng)短邊平均值，那么當(dāng)近場(chǎng)點(diǎn)到單元距離t不超過(guò)0.152 8時(shí)，一定存在擬奇異性。采用十點(diǎn)高斯積分，k=1 rad/m，t分別為0.333 0、0.3、0.2和0.152 8時(shí)，對(duì)應(yīng)的擬奇異積分(分別為IR-1、IR-2、IR-3、IR-4)和誤差(分別為Er-1、Er-2、Er-3、Er-4)隨細(xì)分次數(shù)的收斂性曲線如圖7所示。

圖7 任意單元和近場(chǎng)點(diǎn)的擬奇異積分收斂性曲線Fig.7 The convergence curves of the nearly singular integrals on the arbitrary element

從圖7看出，不細(xì)分時(shí)IR-1和IR-2是收斂的，IR-3和IR-4是不收斂的，表明t為0.333 0和0.3時(shí)不存在擬奇異性，t為0.2和0.152 8時(shí)存在擬奇異性。這符合前文中t不小于0.333 0時(shí)一定不存在擬奇異性和t不超過(guò)0.152 8時(shí)一定存在擬奇異性的預(yù)計(jì)，也顯示了0.152 8

對(duì)于某些典型問(wèn)題，利用臨界接近度能很快判斷是否存在擬奇異性。例如薄板聲輻射計(jì)算中，厚度為多少才算是“薄板”值得探究。實(shí)際聲學(xué)意義上的薄板并不是由板的長(zhǎng)寬厚相對(duì)比例決定的，而是依據(jù)計(jì)算時(shí)單元尺寸與厚度之比。將圖8中的矩形平板劃分一致的矩形單元(不要求單元都為正方形)，那么上下表面的單元是重合的，節(jié)點(diǎn)作為近場(chǎng)點(diǎn)，其相對(duì)位置固定。采用十點(diǎn)高斯積分，根據(jù)計(jì)算時(shí)使用單元的類(lèi)型和臨界接近度得出聲學(xué)薄板的臨界厚度如表4所示，其收斂性已驗(yàn)證，l為單元長(zhǎng)短邊的平均值。厚度低于此臨界值的板為聲學(xué)薄板，需要處理擬奇異性。

圖8 矩形薄板網(wǎng)格Fig.8 The rectangular thin plate grid

單元類(lèi)型臨界厚度/m常單元0．1505l線性單元0．04l二次單元0．1325l不連續(xù)線性單元0．11l不連續(xù)二次單元0．13l

5.2 脈動(dòng)球聲輻射

以水中脈動(dòng)球聲輻射為例驗(yàn)證方法的正確性，將球體劃分242個(gè)四邊形單元，如圖9所示。半徑a=1 m，c=1 500 m/s，ρ=1 000 kg/m3，vn=1 m/s，k=1 rad/m。沿半徑方向距離球心R=1.001～1.1 m平均分布若干近場(chǎng)點(diǎn)，近場(chǎng)點(diǎn)與球體相對(duì)距離為D=(R-a)/a，與球心距離為R的場(chǎng)點(diǎn)聲壓解析解為

(19)

圖9 球體網(wǎng)格模型Fig.9 The spherical grid

圖10 兩種邊界元的計(jì)算對(duì)比曲線Fig.10. The results and errors with two methods

從圖10看出，聲壓幅值曲線p2與p0十分吻合，p1與p0有一定差距，而且D越小差距越大，這一點(diǎn)在誤差曲線Er1和Er2上可以更清晰地看到。在D的整個(gè)區(qū)間上，Er2都是小于Er1的，Er1變化較劇烈，而Er2相對(duì)平穩(wěn)；Er1最大值超過(guò)18%，而Er2一直在1%左右。這些情況說(shuō)明常規(guī)邊界元適用于計(jì)算遠(yuǎn)場(chǎng)聲壓，計(jì)算近場(chǎng)聲壓則難以保證精度，本文的分級(jí)細(xì)分自適應(yīng)邊界元法能有效計(jì)算近場(chǎng)聲壓。二者計(jì)算時(shí)間基本相同，因?yàn)閷?duì)靠近場(chǎng)點(diǎn)的單元進(jìn)行細(xì)分的計(jì)算量相對(duì)于總體計(jì)算量十分微小。對(duì)于聲學(xué)分析，采用平均細(xì)分法、式(15)或“對(duì)角線準(zhǔn)則”所得到的場(chǎng)點(diǎn)聲壓與本文分級(jí)細(xì)分法計(jì)算結(jié)果并無(wú)差別，卻會(huì)產(chǎn)生更多子單元和額外計(jì)算量。

5.3 回轉(zhuǎn)殼聲輻射

如圖11所示，回轉(zhuǎn)殼由半球殼、圓柱殼和圓錐殼組成，半球體半徑為1 m，圓柱體長(zhǎng)7 m，圓錐體高2 m，上端面半徑為0.16 m，表面統(tǒng)一振速為vn=1×10-5m/s，其它參數(shù)同上。分別沿半球、圓柱和圓錐表面單元中心法線設(shè)置3條近場(chǎng)點(diǎn)鏈，距離R為0.001～2 m。三個(gè)單元平均邊長(zhǎng)l分別為0.29、0.39和0.26 m，故滿足式(16)的num皆為4。分別用常規(guī)高斯積分方法和自適應(yīng)積分方法計(jì)算這些近場(chǎng)點(diǎn)聲壓，結(jié)果如圖12所示，其中eR是常規(guī)方法計(jì)算結(jié)果與自適應(yīng)方法計(jì)算結(jié)果的誤差絕對(duì)值。

自適應(yīng)方法中num取4和6的計(jì)算結(jié)果相等，說(shuō)明num為4時(shí)結(jié)果已經(jīng)收斂，圖12中僅展示num取4的曲線?？煽闯鰞煞N積分方法的計(jì)算結(jié)果曲線的近場(chǎng)差別比脈動(dòng)球算例更加顯著， 而且eR降到約1%所需的R也更大，這個(gè)現(xiàn)象說(shuō)明常規(guī)積分方法的近場(chǎng)不精確程度也會(huì)受到結(jié)構(gòu)尺寸的影響。

圖11 回轉(zhuǎn)殼網(wǎng)格模型Fig.11 The rotative shell grid

圖12 3條近場(chǎng)點(diǎn)鏈聲壓和誤差曲線Fig.12 The sound pressures and absolute errors of the 3 chains of near-field points

6 結(jié) 論

對(duì)于邊界元中的擬奇異積分問(wèn)題，本文鑒于自適應(yīng)精確積分的思想提出一種自適應(yīng)分級(jí)細(xì)分方法，不需要迭代和判斷就能夠準(zhǔn)確有效地計(jì)算任意近場(chǎng)點(diǎn)的擬奇異積分，所需子單元更少，且精度可以控制，其靈活性和高效性是一些已有方法所不具備的。接近度的研究為擬奇異積分的計(jì)算和預(yù)估擬奇異性程度提供理論基礎(chǔ)，具有一定的實(shí)用價(jià)值。

收斂性分析發(fā)現(xiàn)：?jiǎn)卧峡拷鼒?chǎng)點(diǎn)的區(qū)域?qū)M奇異性的影響遠(yuǎn)遠(yuǎn)大于遠(yuǎn)離場(chǎng)點(diǎn)的區(qū)域，梯度化地對(duì)待遠(yuǎn)近區(qū)域能較好地考慮奇異性的機(jī)理，節(jié)省計(jì)算量。通過(guò)對(duì)接近度的研究發(fā)現(xiàn)：擬奇異性強(qiáng)弱由場(chǎng)點(diǎn)與單元的相對(duì)位置決定，并非單純由近場(chǎng)點(diǎn)到單元距離決定。

數(shù)值分析表明：在預(yù)估任意不規(guī)則單元和近場(chǎng)點(diǎn)是否存在擬奇異性時(shí)，臨界接近度可作為判斷擬奇異性存在與否的充分條件而非必要條件。對(duì)于薄板聲輻射問(wèn)題，聲學(xué)意義上的薄與厚并不是由板的尺寸相對(duì)比例決定的，而要以接近度為參考，不同類(lèi)型的計(jì)算單元所對(duì)應(yīng)的臨界厚度也不相同。常規(guī)邊界元積分方法近場(chǎng)計(jì)算精度較低，其不精確程度也會(huì)受到結(jié)構(gòu)尺寸的影響。

[1] 李善德, 黃其柏, 李天勻. 新的對(duì)角形式快速多極邊界元法求解聲學(xué) Helmholtz 方程[J]. 物理學(xué)報(bào), 2012,61(6): 64301. LI Shande, HUANG Qibai, LI Tianyun. A new diagonal form fast multipole boundary element method for solving acoustic Helmholtz equation[J]. Acta Physica Sinica, 2012, 61(6): 64301.

[2] 白楊, 汪鴻振. 聲學(xué)-結(jié)構(gòu)設(shè)計(jì)靈敏度分析[J]. 振動(dòng)與沖擊,2003, 22(3):43-45. BAI Yang, WANG Hongzhen. Acoustic-structural design sensitivity analysis[J]. Journal of Vibration and Shock, 2003,22(3):43-45.

[3] 楊迎春, 周其斗. 運(yùn)動(dòng)介質(zhì)中奇異邊界元積分式的精確求解[J]. 振動(dòng)與沖擊, 2011, 30(3): 242-245. YANG Yingchun, ZHOU Qidou. Accurate calculation of weak singular integrals in boundary element method in moving flows[J]. Journal of Vibration and Shock, 2011, 30(3): 242-245.

[4] 張耀明, 孫翠蓮, 谷巖. 邊界積分方程中近奇異積分計(jì)算的一種變量替換法[J]. 力學(xué)學(xué)報(bào), 2008, 40(2): 207-214. ZHANG Yaoming, SUN Cuilian, GU Yan. The evaluation of nearly singular integrals in the boundary integral equations with variable transformation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2008, 40(2): 207-214.

[5] 牛忠榮, 張晨利, 王秀喜. 邊界元法分析狹長(zhǎng)體結(jié)構(gòu)[J]. 計(jì)算力學(xué)學(xué)報(bào), 2003, 20(4): 391-396. NIU Zhongrong, ZHANG Chenli, WANG Xiuxi. Boundary element analysis of the thin-wall structures[J]. Chinese Journal of Computational Mechanics, 2003, 20(4): 391-396.

[6] TANAKA M, SLADEK V, SLADEK J. Regularization techniques applied to boundary element methods[J]. Applied Mechanics Reviews,1994, 47(10): 457-499.

[7] ZHANG Y, GU Y, CHEN J. Stress analysis for multilayered coating systems using semi-analytical BEM with geometric non-linearities[J]. Computational Mechanics, 2011, 47(5): 493-504.

[8] NIU Z, CHENG C, ZHOU H, et al. Analytic formulations for calculating nearly singular integrals in two-dimensional BEM[J]. Engineering Analysis with Boundary Elements, 2007, 31(12): 949-964.

[9] ZHOU H, NIU Z, CHENG C, et al. Analytical integral algorithm in the BEM for orthotropic potential problems of thin bodies[J]. Engineering Analysis with Boundary Elements, 2007, 31(9):739-748.

[10] ZHOU H, NIU Z, CHENG C, et al. Analytical integral algorithm applied to boundary layer effect and thin body effect in BEM for an isotropicpotential problems[J]. Computers & Structures, 2008, 86(1516):1656-1671.

[11] NIU Z, ZHOU H. The natural boundary integral equation in potential problems and regularization of the hypersingular integral[J]. Computers & Structures, 2004, 82(2):315-323.

[12] DAVIES T G, GAO X. Three-dimensional elasto-plastic analysis via the boundary element method[J]. Computers and Geotechnics,2006, 33(3):145-154.

[13] KITA E, KAMIGA N. Error estimation and adaptive mesh refinement in boundary element method, an overview[J]. Engineering Analysis with Boundary Elements, 2001,25(7):479-495.

[14] WEI Y, WANG Y, CHANG S, et al. Numerical prediction of propeller excited acoustic response of submarine structure based on CFD, FEM and BEM[J]. Journal of Hydrodynamics, 2012, 24(2): 207-216.

[15] CROAKER P, KESSISSOGLOU N, MARBURG S. Strongly singular and hypersingular integrals for aeroacoustic incident fields[J]. International Journal for Numerical Methods in Fluids, 2015, 77(5): 274-318.

[16] DONG Y, ZHANG J, XIE G, et al. A general algorithm for the numerical evaluation of domain integrals in 3D boundary element method for transient heat conduction[J]. Engineering Analysis with Boundary Elements, 2015, 51: 30-36.

[17] DAVIES T G, GAO X. Three-dimensional elasto-plastic analysis via the boundary element method[J]. Computers and Geotechnics, 2006, 33(3):145-154.

[18] GAO X W, YANG K, WANG J. An adaptive element subdivision technique for evaluation of various 2D singular boundary integrals[J]. Engineering Analysis with Boundary Elements, 2008, 32(8): 692-696.

[19] LI Y, YE J H. The interaction between membrane structure and wind based on the discontinuous boundary element[J]. Science China Technological Sciences, 2010, 53(2): 486-501.

[20] 王靜, 高效偉. 基于單元子分法的結(jié)構(gòu)多尺度邊界單元[J].計(jì)算力學(xué)學(xué)報(bào), 2010, 27(2): 258-263. WANG Jing, GAO Xiaowei. Structural multi-scale boundary element method based on element subdivision technique[J]. Chinese Journal of Computational Mechanics, 2010, 27(2): 258-263.

[21] QU S, CHEN H, LI S. Adaptive integration method based on sub-division technique for nearly singular integrals in near-field acoustics boundary element analysis[J]. Journal of Low Frequency Noise, Vibration and Active Control, 2014,33(1): 27-46.

[22] GAO X W,DAVIES T G. Adaptive integration in elasto-plastic boundary element analysis[J]. Journal of the Chinese Institute of Engineers, 2000, 23(3):349-356.

[23] ZHANG J, TANAKA M, MATSUMOTO T. Meshless analysis of potential problems in three dimensions with the hybrid boundary node method[J]. International Journal for Numerical Methods in Engineering, 2004, 59(9): 1147-1166.

[24] ZHANG J, QIN X, HAN X, et al. A boundary face method for potential problems in three dimensions[J]. International Journal for Numerical Methods in Engineering, 2009, 80(3): 320-337.

[25] KANG S C, IH J G. On the accuracy of nearfield pressure predicted by the acoustic boundary element method[J]. Journal of Sound and Vibration, 2000, 233(2): 353-358.

Adaptive method for calculating nearly singular integrals in acoustic boundary elements

MIAO Yuyue1,2, LI Tianyun1,2,3, ZHU Xiang1,2, ZHANG Guanjun1,2, GUO Wenjie1,2

(1. School of Naval Architecture and Ocean Engineering，Huazhong University of Science & Technology, Wuhan 430074,China;2.Hubei Key Laboratory of Naval Architecture & Ocean Engineering Hydrodynamics, Wuhan 430074,China;3.National Key Laboratory on Ship Vibration & Noise, Wuhan 430064,China)

An adaptive method was presented to calculate nearly singular integrals in acoustic boundary elements. The element was subdivided into subelements hierarchically so as to remove the near singularity by transforming the integral over the initial element to Gauss integrals over subelements. Based on this method, the near singularity was thoroughly studied and the concept of proximity was described. The critical proximity can be used as the criterion for calculating nearly singular integrals and to predict the existence of near singularity. Compared with previous methods, the adaptive method is more flexible and efficient. The integral precision can be adjusted and is not restricted by the positions of near-field points. It is found that the positional relationships between near-field points and the element have important effect on nearly singular integrals. The work provides more understanding and choice to the problem of near singularity in acoustic boundary elements.

boundary element; nearly singular intergral; adaptivity; proximity

國(guó)家自然科學(xué)基金(51379083;51479079);高等學(xué)校博士學(xué)科點(diǎn)專(zhuān)項(xiàng)科研基金(20120142110051)

2015-09-09 修改稿收到日期：2015-12-05

繆宇躍 男，博士生，1988年生

李天勻 男，教授，博士生導(dǎo)師，1969年生 E-mail:ltyz801@mail.hust.edu.cn

TB532

A

10.13465/j.cnki.jvs.2017.02.004