雙層介質(zhì)聲傳播問題的標準解及有限元解

周益清，駱文于，吳雙林

（1. 中國科學(xué)院聲學(xué)研究所聲場聲信息國家重點實驗室，北京 100190；2. 中國科學(xué)院大學(xué)，北京 100049）

0 引 言

自20世紀60年代以來，國內(nèi)外學(xué)者對聲場計算方法進行了大量的研究，先后出現(xiàn)了許多種針對水下聲傳播的計算模型。目前比較流行的幾種聲場計算方法有：射線法[1-2]、簡正波法[3-4]、拋物方程法[5-6]、波數(shù)積分法[7-8]，研究人員使用這些方法建立了許多聲場計算模型，但是目前還沒有非常流行的有限元算法模型。可見，關(guān)于有限元方法在水下聲場計算中的研究還不是非常充分。

有限元方法是求解偏微分方程復(fù)雜邊值問題的一種通用數(shù)值解法。它首先將物理域離散成有限數(shù)量的單元，各單元以節(jié)點相連，然后利用各單元間的關(guān)聯(lián)性形成一系列的有限元方程組，求得單元內(nèi)的精確解或近似解，從而把復(fù)雜的偏微分方程邊值問題轉(zhuǎn)化為大型方程組的求解問題[9]。在水聲學(xué)中，目前有限元方法主要被用來分析換能器的結(jié)構(gòu)振動、解決聲散射問題[10]，而關(guān)于聲傳播的研究[11-12]相對較少。

由于過去的計算機性能較現(xiàn)在落后很多，且我們常關(guān)心大尺度的聲傳播問題對計算機的計算速度和內(nèi)存的要求都很高，所以過去的水下聲傳播方向的學(xué)者對有限元方法的研究較少。但是，近些年，隨著計算機軟硬件的快速發(fā)展，已經(jīng)有充分的條件使用有限元方法進行水下聲傳播的研究。因此，希望開發(fā)一個基于有限元方法的聲場計算模型，以處理復(fù)雜邊界下的實際水下聲傳播問題。

本文主要研究雙層介質(zhì)中的聲傳播問題。首先，討論有限深度的雙層介質(zhì)中的聲場計算，使用波數(shù)積分方法計算出標準解后，再使用有限元方法計算相同問題的聲場，對于出射邊界，我們分別選擇傳統(tǒng)輻射條件、Dirichlct to Ncumann （DtN）非局部算子和完美匹配層來處理，經(jīng)過對比，發(fā)現(xiàn)使用完美匹配層方法處理出射邊界精度更高。然后，討論無限深度的雙層介質(zhì)，即 Pckcris波導(dǎo)中的聲場計算，這里同樣使用波數(shù)積分方法計算標準解，再與結(jié)合了完美匹配層技術(shù)的有限元方法以及常用的簡正波模型KRAKEN進行對比。我們發(fā)現(xiàn)，當某號簡正波的本征值非常接近割線時，KRAKEN可能無法對該本征值進行準確求解，從而導(dǎo)致較大的計算誤差。但是有限元方法則不受此類限制，所以相比KRAKEN而言，有限元方法的普適性更優(yōu)。

1 問題描述

1.1 有限深度的雙層介質(zhì)問題

如圖1所示為有限深度的雙層介質(zhì)，上層介質(zhì)的下邊界深度為D1，聲速和密度分別為c1，ρ1，下層介質(zhì)的下邊界深度為D2，聲速和密度分別為c2，ρ2。在z=0處為絕對軟的上邊界，聲壓為0；在z=D2處為絕對硬的下邊界，質(zhì)點法向振速為0。對于無限大三維介質(zhì)中的無限長線聲源，可以將它簡化為圖1所示的二維問題，假設(shè)聲源位于第一層介質(zhì)中，深度為zs。

圖1 有限深度的雙層介質(zhì)環(huán)境示意圖Fig.1 Schematic diagram of a two-layer medium with bothlimited depths

1.2 無限深度的雙層介質(zhì)問題

如圖2所示為無限深度的分層介質(zhì)，即Pckcris波導(dǎo)。上層介質(zhì)的下邊界深度為D1，聲速和密度分別為c1，ρ1，下層介質(zhì)厚度無限大，聲速和密度分別為c2，ρ2。在z=0 處為絕對軟的上邊界，聲壓為0。對于無限大三維介質(zhì)中的無限長線聲源，同樣可以將它簡化為圖2所示的二維問題，聲源深度為zs，假設(shè)位于第一層介質(zhì)中。

圖2 無限深度的雙層介質(zhì)環(huán)境示意圖Fig.2 Schematic diagram of a two-layer medium with an unlimited depth

由于這兩個問題都不存在解析解，所以我們使用波數(shù)積分方法推導(dǎo)出標準解，用于與其他方法進行比對。

2 使用波數(shù)積分法計算標準解

2.1 有限深度的雙層介質(zhì)的標準解

對于圖1所示的無限長線源問題，選用直角坐標系（x，z），z軸通過聲源且垂直向下，x軸平行于海面，此時線源的Hclmholtz方程表示為[9]

其中，ψ（x，z）表示位移勢。對于線源問題，利用傅里葉變換對：

得到如下形式的深度分離波動方程：

帶入邊界條件，求解該方程，得到深度格林函數(shù)Ψ（kx，z），利用逆傅里葉變換求得位移勢ψ（x，z），進而求得總聲場。

對于深度分離的波動方程，我們只需要考慮介質(zhì)中的上行波及下行波。圖3為雙層介質(zhì)中的上行波與下行波示意圖，在每層介質(zhì)中，都分別包含上行波與下行波。其中，上層介質(zhì)中的下行波幅值記為，上行波幅值記為；下層介質(zhì)中的下行波幅值記為，上行波幅值記為。

圖3 有限深度雙層介質(zhì)中上行波與下行波示意圖Fig.3 Schematic diagram of up-going and down-going waves in the two-layer medium with both limited depths

2.2 無限深度的雙層介質(zhì)的標準解

對于無限深度的雙層介質(zhì)，即 Pckcris波導(dǎo)中的聲場，同樣也是首先求解式（4）中的深度分離的波動方程。

對于深度分離的波動方程，只需要考慮上行波及下行波。圖4表示雙層介質(zhì)中的上行波與下行波示意圖，在上層介質(zhì)中，既包含上行波，又包含下行波；而在下層介質(zhì)中，只存在下行波。其中，上層介質(zhì)中的下行波幅值記為，上行波幅值記為；下層介質(zhì)中下行波幅值記為。

圖4 無限深度雙層介質(zhì)中上行波與下行波模擬示意圖Fig.4 Schematic diagram of up-going and down-going waves in the two-layer medium with an unlimited depth

3 使用有限元方法計算聲場

使用有限元方法計算聲場主要分為以下步驟：（1） 將物理域離散為許多小的計算單元；（2） 根據(jù)單元的形狀和結(jié)點數(shù)構(gòu)造形函數(shù)；（3） 寫出近似解表達式；（4） 推導(dǎo)出離散的有限元方程；（5） 建立單元剛度矩陣和單元外載荷向量；合成總體剛度矩陣和載荷向量；（6） 帶入邊界條件求解聲場；（7） 后處理和誤差分析。

3.1 有限元方程的推導(dǎo)

3.2 計算區(qū)域離散后的有限元方程

由于計算區(qū)域形狀規(guī)則，所以采用四節(jié)點四邊形單元對計算區(qū)域進行離散。同時，四節(jié)點四邊形單元比三節(jié)點三角形單元具有更高的精度，擁有更高階的連續(xù)性。四節(jié)點四邊形單元如圖5所示，單元中點坐標為 （x0，z0） ，單元的長和寬分別為2a，2b。

圖5 四節(jié)點四邊形單元的局部編號、中點坐標、長度與寬度示意圖Fig.5 Schematic illustration of local number, midpoint coordinates, length and width of the four-node quadrilateral element

用函數(shù)Nn（x，z） 表示形函數(shù)，其中n=1 ，2，3，4分

3.3 出射邊界的處理

在水下聲傳播計算方面，有限元方法面對的最大的挑戰(zhàn)之一是如何模擬沒有邊界的無限大海洋環(huán)境，即滿足無限遠處的Sommcrfcld輻射條件[13]：

其中：R表示距離；k表示波數(shù)；p表示聲壓。下面，分別使用三種常用的出射邊界計算方法來求解聲場。

3.3.1 傳統(tǒng)輻射條件

傳統(tǒng)的輻射條件是指：待求量的徑向?qū)?shù)等于待求量在聲源一定距離處的值與常數(shù)的乘積。對于本文討論的雙層介質(zhì)中的聲傳播問題，可以將右側(cè)出射邊界用傳統(tǒng)輻射條件表示為

代入式（30）的第一項，可以得到：

將式（33）中的聲壓p和式（34）中的試驗函數(shù)q代入式（38）即可求得最右側(cè)單元對應(yīng)的單元剛度矩陣。

最后，利用式（35）計算得到的單元剛度矩陣，在集成總體剛度矩陣時，將每個單元的局部坐標與總體坐標一一對應(yīng)，即可得到總體剛度矩陣。

3.3.2 DtN非局部算子

傳統(tǒng)輻射條件并不能完全吸收出射波，回波會干擾聲場，產(chǎn)生誤差。1978年，F(xiàn)ix和 Marin[14]對軸對稱情況下的分層介質(zhì)中的聲場進行了研究，并且第一次提出了非局部算子，同時指出了非局部算子相對于傳統(tǒng)輻射條件的優(yōu)越性。

3.3.3 完美匹配層

傳統(tǒng)輻射條件可能會引入較大誤差，非局部算子會改變?nèi)謩偠染仃嚨南∈栊?。下面引入完美匹配層，在保持剛度矩陣稀疏性的同時，盡量減小誤差。

完美匹配層是Bcrcngcr[15-16]在1994年提出的一種概念，希望用它來模擬出射邊界，使得只有少量的能量甚至沒有任何能量反射回我們計算的物理場中。

如圖6所示，對于深度有限的雙層介質(zhì)，我們設(shè)需要計算的物理域的水平距離為R，在物理域的右側(cè)設(shè)置一個厚度為δx的完美匹配層。如圖6中右側(cè)PML1和PML2所示，完美匹配層中的聲速和密度與相鄰的左側(cè)介質(zhì)相同。

圖6 添加完美匹配層的有限深度雙層介質(zhì)示意圖Fig.6 Schematic diagram of two-layered medium with both limited depths plus perfectly matched layers

而對于無限深度的雙層介質(zhì)，設(shè)需要計算的物理域的水平距離為R，深度為D2，在物理域的右側(cè)設(shè)置一個厚度為δx的完美匹配層，在物理域的下側(cè)設(shè)置一個厚度為δz的完美匹配層，如圖7所示。

圖7 添加完美匹配層的無限深度雙層介質(zhì)示意圖Fig.7 Schematic diagram of two-layer medium with an unlimited depth plus perfectly matched layers

3.4 其他邊界條件的處理

對于左側(cè)邊界，即x=0處，滿足對稱條件，可視作位移為零的絕對硬邊界；對于有限深度的雙層介質(zhì)z=D2處的下邊界也是絕對硬邊界；z=0 處的上邊界是絕對軟邊界，聲壓p=0 。

對于絕對硬邊界，法向位移為0，即?p/?n=0，對剛度矩陣貢獻為0，不用進行特殊的處理和計算。對絕對軟邊界，我們采用“對角元素置 1”法進行計算，步驟如下：

（1） 設(shè)總體剛度矩陣為K，其元素可以表示為Kmn，外載荷向量為f，聲壓向量為p，它們滿足線性方程組Kp=f，設(shè)絕對軟邊界的節(jié)點總體編號為l。

（2） 將K的第l列與第l行置零，第l個對角元素置1。

（3） 將f的第l個元素設(shè)置為 0。

（4） 得到更新后的線性方程組Kp=f，求解即可得到每個節(jié)點處的聲壓。

4 計算結(jié)果

4.1 有限深度雙層介質(zhì)中的聲場

取如圖8所示的環(huán)境參數(shù)，上層介質(zhì)的厚度為 60 m，聲速c1=1500 m· s-1，密度ρ1=1 000 kg·m-3，下層介質(zhì)厚度為 40 m，聲速c2=1800 m· s-1，密度ρ2=1 800 kg·m-3，聲源深度為20 m。

圖8 雙層介質(zhì)環(huán)境參數(shù)示意圖Fig.8 Specific environmental parameters of a two-layer medium

圖9表示使用不同方法計算的25 Hz聲源產(chǎn)生的聲場。其中，圖9（a）為使用波數(shù)積分方法計算得到的解，可視為標準解。圖9（b）是使用傳統(tǒng)輻射條件處理出射邊界的有限元解。可以看出有限元解與標準解的圖案基本一致，但是在水平方向存在明顯的干涉條紋。這是由于傳統(tǒng)輻射條件沒有完全吸收出射波，向右傳播的出射波與向左傳播的反射波產(chǎn)生干涉現(xiàn)象，從而引入了誤差。圖9（c）是使用DtN非局部算子處理出射邊界的有限元解?？梢钥吹?，相比傳統(tǒng)輻射條件而言，干涉現(xiàn)象有了輕微的削弱，但是依然存在，且不可忽略。同時，DtN非局部算子只能用于處理單向出射邊界，即難以處理可穿透海底的問題。圖9（d）是使用完美匹配層處理出射邊界的有限元解，可以看到，該計算結(jié)果與標準解吻合得較好。

圖9 聲源頻率為25 Hz時使用不同方法計算的雙層介質(zhì)中的聲場Fig.9 Sound field diagrams in the two-layer medium calculated by different methods at 25 Hz

圖10表示深度為48 m處的傳播損失。該深度下，有限元解和標準解的平均誤差為 0.09 dB，兩組計算結(jié)果吻合得較好。

圖10 聲源頻率為25 Hz時接收深度48 m處的傳播損失Fig.10 Transmission loss at 25 Hz and at the receiving depth of 48 m

圖11 聲源頻率為25 Hz時接收深度80 m處的傳播損失Fig.11 Transmission loss at 25 Hz and at the receiving depth of 80 m

從以上結(jié)果可以看出，本文所提的有限元模型計算得到的結(jié)果與標準解吻合得較好。下面我們需要驗證，對于更高的頻率，該有限元模型是否依然具有較高的精度。

圖12表示聲源頻率為100 Hz時分別使用波數(shù)積分方法和有限元方法計算的聲場，其中圖12（a）表示使用波數(shù)積分方法計算得到的聲場標準解，圖12（b）是使用完美匹配層處理出射邊界的有限元解。從計算結(jié)果可以看出，本文所提的有限元解與波數(shù)積分方法計算的標準解吻合得較好。

圖12 聲源頻率為100 Hz時使用不同方法計算的雙層介質(zhì)中的聲場Fig.12 Sound field diagrams in the two-layer medium calculated by different methods at 100 Hz

圖13為深度60 m處的傳播損失。該深度下，因為存在一個傳播損失波谷，所以有限元解和標準解的平均誤差稍大，近似為0.19 dB。兩組計算結(jié)果總體上吻合得較好。

圖13 聲源頻率為100 Hz時接收深度60 m處的傳播損失Fig.13 Transmission loss at 100 Hz and at the receiving depth of 60 m

4.2 無限深度雙層介質(zhì)中的聲場

對于無限深度的雙層介質(zhì)，即 Pckcris波導(dǎo)，我們?nèi)∪鐖D14所示的環(huán)境參數(shù)，上層介質(zhì)的厚度為40 m，聲速c1=1500 m· s-1，密度ρ1=1 000 kg·m-3；下層介質(zhì)無限大，聲速c2=1650 m· s-1，密度ρ2=1 500 kg·m-3。聲源深度為20 m，接收深度為40 m。

圖14 無限深度雙層介質(zhì)環(huán)境參數(shù)示意圖Fig.14 Specific environmental parameters of a two-layer medium with an unlimited depth

首先計算聲源頻率為105～125 Hz的傳播損失。如圖15所示為不同聲源頻率在接收深度40 m處的傳播損失，其中圖15（a）為常用的簡正波模型KRAKEN計算的傳播損失，圖15（b）為使用波數(shù)積分方法推導(dǎo)出的標準解，圖15（c）為使用有限元方法計算的傳播損失。從圖15可以看出，有限元解與波數(shù)積分方法計算的標準解吻合得較好，而KRAKEN計算的結(jié)果在某些頻率上存在顯著誤差，這是由于對聲場有主要貢獻的某號本征值位于割線附近時，KRAKEN無法準確計算該號本征值導(dǎo)致的[18]?？梢钥闯鲈诼曉搭l率為 112、118 Hz時，KRAKEN計算的傳播損失與另外兩種方法計算的傳播損失有明顯的差別。下面分別分析這兩個頻率處的本征值與傳播損失的關(guān)系。

圖15 聲源頻率為105～125 Hz時不同方法計算的接收深度為40 m處的傳播損失Fig.15 Transmission loss diagrams calculated by different methods at the receiving depth of 40 m and at the frequencies from 105 to 125 Hz

圖16為不同聲源頻率時使用波數(shù)積分方法計算的深度格林函數(shù)。圖16（a）中聲源頻率為112 Hz，對應(yīng)的下層介質(zhì)波數(shù)約為0.426 5 m-1。在這個頻率下，KRAKEN得到的前五號簡正波的水平波數(shù)（即本征值）分別為 0.464 455、0.449 802、0.381 011+i0.009 161、0.309 880+i0.017 855、0.188 504+i0.039 225，對于當前問題，由于第3號簡正波的水平波數(shù)與海底波數(shù)非常接近，KRAKEN未能計算出該號簡正波，從而導(dǎo)致KRAKEN模型的傳播損失結(jié)果出現(xiàn)了較大的誤差。圖16（b）中聲源頻率為 118 Hz，對應(yīng)的下層介質(zhì)波數(shù)約為0.449 3 m-1。在此頻率下，KRAKEN得到的前五號簡正波的水平波數(shù)分別為 0.489 759、0.475 685、0.451 877、0.411 541+i0.007 906、0.346 650+i0.015 463。比較簡正波的水平波數(shù)與圖16（b）可以發(fā)現(xiàn)，聲源頻率為118 Hz時，第3號簡正波的水平波數(shù)與海底波數(shù)的差異足夠大，KRAKEN可以準確求出這號簡正波，因此在聲源頻率為118 Hz時，KRAKEN計算結(jié)果與標準解的一致性非常好。

圖16 不同聲源頻率對應(yīng)的深度為40 m處的深度格林函數(shù)，F(xiàn)ig.16 Magnitudes of the depth-dependent Green’s functions at the depth 40 m and at the frequencies of 112 Hz and 118 Hz

然后分別對聲源頻率為112、118 Hz時，三種不同方法計算的傳播損失進行比較。圖17（a）、17（b）分別為聲源頻率為112、118 Hz時接收深度40 m處的傳播損失。可以看到，聲源頻率為 112 Hz時，KRAKEN遺漏了有主要貢獻的第三號簡正波，導(dǎo)致計算誤差較大，而有限元方法不需要計算簡正波，與標準解吻合得非常好。聲源頻率為 118 Hz時，KRAKEN包含了前三號有主要貢獻的簡正波，所以與有限元解和標準解吻合得非常好。

圖17 不同方法計算的在40 m深度不同聲源頻率傳播損失結(jié)果Fig.17 Transmission losses calculated by different methods at the depth of 40 m and at the frequencies of 112 Hz and 118 Hz

通過對比發(fā)現(xiàn)，簡正波模型KRAKEN需要求解簡正波，而在一些特定的情況下，可能會遺漏對聲場有主要貢獻的簡正波，從而產(chǎn)生較大誤差。而有限元模型不需要求解簡正波，所以對幾乎所有頻率都具有較高的精度，即有限元方法的普適性優(yōu)于簡正波方法。

5 結(jié) 論

本文研究了雙層介質(zhì)中的聲傳播問題，分別討論了有限深度的雙層介質(zhì)和無限深度的雙層介質(zhì)的情況。

由于有限深度雙層介質(zhì)問題不存在解析解，因此首先使用波數(shù)積分方法推導(dǎo)出標準解。然后提出一個有限元模型，并研究不同出射邊界的處理方法對有限元解的影響，發(fā)現(xiàn)傳統(tǒng)輻射條件會引入較大的誤差，非局部算子會破壞總體剛度矩陣的稀疏性，而且不適用于沿不同方向出射的聲波，而完美匹配層擁有較高的精度，且不會破壞總體剛度矩陣的稀疏性，使用也十分靈活。本文所提的有限元模型選取精度最高、使用方法最靈活的完美匹配層方法來處理出射聲場，并進行后續(xù)計算。數(shù)值結(jié)果表明本文所提的有限元解和標準解吻合較好。

對于無限深度的雙層介質(zhì)，同樣先使用波數(shù)積分方法推導(dǎo)出標準解，然后分別使用本文所開發(fā)的結(jié)合了完美匹配層的有限元模型和常用的簡正波模型KRAKEN，計算不同頻率下的聲場。研究發(fā)現(xiàn)在一些特定的頻率，由于KRAKEN在計算時遺漏了一些對聲場有主要貢獻的簡正波，會產(chǎn)生明顯的誤差。而有限元方法不需要計算簡正波，因此在各頻率都與標準解吻合得較好。

本文基于傳統(tǒng)的Galcrkin方法推導(dǎo)出水下聲傳播有限元方程，然后使用完美匹配層處理出射聲場，設(shè)計并實現(xiàn)了一個水下聲場計算有限元模型。通過該模型計算得到的雙層介質(zhì)中的聲場，與標準解吻合得較好。同時，由于有限元方法不需要求解簡正波，因此對于某些簡正波方法不能準確處理的問題，有限元方法仍能給出準確的聲場結(jié)果，表明有限元方法的普適性優(yōu)于簡正波方法。

值得一提的是，有限元方法的優(yōu)缺點都非常突出，優(yōu)點在于普適性強且精度高，缺點在于計算量大且物理意義不清晰。因此，本文選擇研究有限元方法并非是希望建立一個通用的聲場計算模型來解決實際應(yīng)用中的聲傳播問題，而是希望使用有限元方法提供參考解，來優(yōu)化其他聲場計算模型。例如，在計算復(fù)雜海洋環(huán)境中的聲場時，如復(fù)雜地形或內(nèi)波等問題，很難驗證現(xiàn)有模型的計算結(jié)果的正確性，這時就可以使用有限元方法提供標準解，用于校驗現(xiàn)有模型的結(jié)果，從而對其進行改進，使其適用于解決更加復(fù)雜的環(huán)境中的聲場問題。

本文的內(nèi)容是有限元方法在水下聲場計算中的一些初步工作，接下來將對本文所提的有限元聲場計算模型進一步改進，提升模型精度，并使之可以計算更加復(fù)雜環(huán)境中的聲場。