二維楔形海域之聲波傳播數值計算

行政院國家科學委員會專題研究計畫成果報告

二維楔形海域之聲波傳播數值計算

Computation on the two-dimensional sound propagation

in a wedge-shaped ocean

計畫編號：NSC 89-2611-E-002-012

執行期限：88 年 8 月 1 日至 89 年 7 月 31 日

主持人：謝傳璋

國立台灣大學造船及海洋工程學系

一、中文摘要 本文以有限元素方法，求解一線聲源 在二維楔形淺海域中聲場變化，分析楔頂點 附近陰影區與聲源頻率、楔頂角大小，聲源 深度，聲源位址等參數之關係。 關鍵詞：有限元素法、楔形海域、聲波傳播 Abstract

In order to investigate the sound propagation properties of the shadow zone, the finite element method will be used to compute the sound pressure distribution of a wedge-shape (of rigid sea bottom) sea due to a line source. The relationships of the size of the zone among those parameters, say, the location and frequency of the source, the apex angle… are analyzed.

Keywords: finite element method,

wedge-shape sea , sound propagation 背景介紹 在 range-dependent 海洋聲傳問題中， 由於聲傳模型日趨完備而漸漸受到重視，不 論是以傳統數學解析方法或數值計算均能得 到相當一致的結果。楔形聲傳問題中，陰影 區的存在無庸置疑，但於文獻中吾人只能找 到以不同方式求解聲場分佈特性，對於陰影 區受那些參數影響，而各項參數影響程度如 何則未見討論，本文將以數值方法對楔形聲 場中陰影區特性做初步探討。 數值解法大致可分成三類，分為別 PE （Parabolic Equation）、正態解（normal mode）、有限元素/有限差分法，當中 PE 法為 三種方法中通常須要較長久的計算時間，特 別是在聲場中存在有陰影區（shadow zone） 時，另外，PE 法在使用上也有較多限制。正 態解則對於計算時所需模態數目的多寡，必 須經過多次計算比較後，方能得到適當模 型，猶其在愈接近聲源時計算量將隨之增 加。有限元素及有限差分法具有使用方便、 較少量計算時間等優點。 研究動機 在一個二維楔形淺海域內，如下圖所示，若 其楔頂角為θ ，當一點聲源或線聲源置於其 海域時，聲源所發出上坡方向傳播射線，在 聲線在楔形海域行程意識圖 經歷一次自由液面及海床反射後，此入射聲 線之入射掠角將增為2θ，經過多次反射後， 此射線將開始往下坡方向傳遞，亦就是說， 在接近楔頂角附近聲波將無法到達，此區域 稱為陰影區。另外，因聲波每次碰撞海床後， 總會損失其強度，特別在接近楔頂角時，傳 輸損失將更加明顯，是故將更形強化陰影區 產生。

本文將以有限元素方法，求解一線聲源在二 維楔形淺海域中聲場變化，分析楔頂點附近 陰影區與聲源頻率、楔頂角大小，聲源深度， 聲源位址等參數之關係。 文獻回顧 楔形聲場之研究及計算，最早可能是在西元 １８９６年由 Sommerfeld1 提出，該文以理 想反射邊界組成楔形區域內，由點聲源或稱 為線聲源所引發擾動，以分離變數法求解非 齊性 Helmholtz 方程式，進而得到正態解 （normal mode）。 於１９８８年，Buckingham 與 Tolstoy2 在一 均質介質內，以聲壓釋放面為邊界，求解非 齊性 Helmholtz 方程式，經 Fourier 及 Hankel 轉換，而求得非耦合（uncouple）正 態速度位勢（velocity potential）解

(

)

( )

( )

∑

∞ = = Ψ 1 ' ' 0 sin sin 2 m v r r v v I θ θ θ ， ， （1.1） 其中 0 θ π m v= ， （1.2）

( )

_∫

∞

( )

( )

− = 0 ' 2 2 ' dp pr J pr J k p p r r I_v ， _{v v} ， （1.3） 上式中 k 為聲源波數（wavenumber），J 為 v_v 階第一類 Bessel 方程，θ 為楔頂角，r 與_o ' r 分別為收聲器及發聲器至楔頂角矩離，θ 與 ' θ 為收聲器及發聲器至楔頂角連線與水面夾 角，（1.1）至（1.3）式中各符號幾何意義如 圖 1.1 所示。將 1.3 式代入 1.1 式得到最終 速度位勢

(

)

_∑

∞

( )

( )

_{( ) ( )}

( )

= < > = Ψ 1 ' 1 0 ' ' sin sin , , , m v v kr H kr v v J i r r θ θ θ π θ θ ， （1.4） 其中

( )

' , min r r r_< = ，

( )

' , maxr r r_> = ，

( )

12 1 − = i 及 ( )1

( )

v H 為 v 階第一類 Hankel 方程式。1.4 式適用於任何楔頂角，亦即不限定楔頂角為 π 的約數。 發聲器信號產生由下式表示

( )

i _H( )

_{( )}

_kR R 1 0 4     = Φ ， （1.5） 當中 R 為聲場中某一點與聲源的距離， ( )1 0 H 為零諧第一類 Hankel 方程式。為比較上方 便，將聲場由下式正常化

(

)

( )

1 ' ' Φ Ψ = Λ r，r，θ，θ ， （1.6） 並定義傳輸損耗為

( )

Λ =20log10 TL ，（1.7） 此一結果將在第三章中與程式執行結果做比 對。 Glegg 與 Yoon3 對 Buckingham 理論建立一水 槽加以驗證，同時對於在楔頂角附近陰影區 （shadow zone）與聲源深度、收聲器與發聲 器橫向距離、聲源波長做了比較。 在使用數值法求解楔形聲場時，多數文獻利 用拋物線方程式（Parabolic Equation, PE） 法，此種方法首先由 Tapper4_{所推得。考慮} 一聲場，在以圓柱座標為系統，可將聲波方 程表示為

(

r z

)

t

(

r z t

)

F

(

r z t

)

c , , , , , s , , , 1 2 2 2 2 φ θ θ θ  =    ∂ ∂ − ∇ （1.8）

其中φ 為 scalar wavefield potential，

(

r zt

)

F_s ,θ, , 為聲源。若（1.8）式中聲場為軸 對稱，即φ 與θ 項無關，可將（1.8）式簡化 θ α+2 α θ

為

( )

r z t

(

r z t

)

F

(

r z t

)

c , , , s , , 1 2 2 2 2 =     ∂ ∂ − ∇ φ （1.9） 利用 Fourier 轉換對

( )

{ }

( )

( )

( )

( )

( )

      ≡       ℑ ≡ ≡ ℑ ≡

∫

∫

∞ ∞ − − ∞ ∞ − − ^ ^ 1 ^ 2 1 ω ω ω π ω ω ω d e f f t f dt e t f t f f t j t j （1.10） 可將（1.9）表成

{ }

{ } { }

tt Fs c ℑ =ℑ − ℑ ∇2 φ 1₂ φ （1.11） 在對上式分別運算後，可整理得

( )

(

,

)

φ

(

, ,ω

)

(

, ,ω

)

^ ^ 2 2 z r F z r z r k_m = s + ∇ （1.12） 上式中，km =ω c

( )

r,z 介質波數（media wavenumber）。假若聲場環境不依水平距離而 變，即c

( ) ( )

r,z =c z ，同時將聲源置於r =0 處，而聲源為線聲源，那麼（1.12）式可表 示為

( )

(

)

( )

( ) ( )

r r z f S z r z k s m _π δ φ ω 2 , ^ 2 2+ = ∇ （1.13） 再經 Hankel 轉換，便可將上式由 r-domain 轉換至 k-domain，將（1.13）式做一次零階 Hankel 轉換，便得

( )

(

)

φ

( )

ω _π

( )

2 , 2 2 2 2 z f S z k z k k dz d s m  =      − − − （1.14） 上式為一特解，再加上齊性方程式通解

( )

(

2 2

)

( )

, 0 2 2 =       _{− −} − z k z k k dz d m φ （1.15） 即

( )

( )

( ) ( )

( ) ( )

− − − + − − − − + + = k z A k k z A k k z z k, φ* , φ , φ , φ （1.16） （1.16）式中，A−

( )

k 及A+

( )

k 為邊界條件求 得之係數。最後將上式依次做 Hankel 反轉換 及 Fourier 反轉換，可得到時域反應（time response）φ

(

r,z,t

)

解。 Jensen 與 Ferla5 使用了三種不同運算 法求解楔形聲場，分別為（一）coupled-mode、（二）有限差分 PE 及（三）slipt-step PE，對不同底層性質與不同聲剖面做測試， 在多次測試各項參數及證實程式收歛後，得 知在這三種方法中以雙向（two-way） couple-mode 最為準確，但所須執行時間也 最長。 數 學 模 型 由於本文著重於數值計算，並使用既有 程式為工具，故本章將只簡略陳述統御方程 式及其邊界條件，另外也會大略說明使用 Matlab PDE toolbox 方法及本文程式執行步 驟，在本章最後一節將簡略描述本方程式所 將進行的各項參數，以期對整個測試項目有 概觀認識。 統御方程式 圖 2.1 為問題幾何外形，線聲源平行於楔頂 點，介質為均勻 流體，聲源產生一簡諧信號

(

_{x y}

)

_e i t f =δ ， −ω ，其中δ 為 delta 方程式， 則聲波U

(

x,y,t

)

之波動方程式可寫成

(

)

i t e y x f U c t U _{− ∇ = =δ} −ω ∂ ∂ ， 2 2 2 2 ， （2.1） 令

(

) (

)

i t e y x p t y x U ， ， = ， −ω ，i= −1，則 （2.1）式可表為 Helmholtz 方程式

(

x y

)

p k p− =δ ， ∇ − 2 2 _{， （2.2）} 表中 k 為波數，其相對頻率為 f ，波長為λ ， k 與 f 之間關係為 λ π π ω ₌ 2 ₌ 2 = c f c k ， （2.3） 圖 2.1 中，邊界條件為 在 上，自由液面邊界條件：p=0

（2.4） 在 上，剛性海床邊界條件：→n⋅∇p=0 （2.5） 在 上：n→⋅∇p=ikp （2.6） 邊界條件 3 稱為 Sommerfeld 邊界條件，或 稱為幅射邊界條件。其中→n 為各邊界上之單 位法線向量。 有限元素方法 在不失一般狀況下可將（2.2）式寫成

( )

c∇u +au = f ⋅ ∇ − in Ω （2.7） 其中Ω為在平面的定義域， c 、 a 、 f 及待 定解 u 為在平面區域Ω之複數方程式。邊界 條件係在邊界上，由 u 與 u 的導數混合而 成，可分成三類： Dirichlet：hu=r on boundary ∂Ω。 Generalized Neumann：n→⋅

( )

c∇u +qu= gon Ω ∂ 。

混合型：Dirichlet 及 generalized Neumann 的混合。 以有限元素法求解橢圓型偏微分方程近 似解大致可分成三個步驟： 確定定義域Ω及邊界條件。 網格產生。 偏微分方程式及其邊界條件離散化，而得一 聯立代數方程組Ku=F，其中 u 為在網格點 上的近似解。 在進行理論推導前，為避免邊界應用時 的困擾，在 Matlab PDE users guide 手冊中， 將 Dirichlet 邊界條件擴充成一般化 Neumann 型邊界條件，因此，所有邊界條件 只有一種型式，這對程式寫作將較為方便。 假設在（2.7）式中，u 為該式的解，在 兩端乘以一測試方程式 v ，並在定義域Ω內 積分，即

( )

(

)

[

]

_∫

[ ]

∫

Ω − ∇⋅ c∇u v+auvdx= Ω fvdx （2.8） 對上式做一次部份積分

( )

[

]

_∫

( )

_∫

∫

∂Ω Ω → Ω c∇u ⋅∇v+auvdx− n⋅ c∇u vds= fvdx （2.9） （2.9）式中，第二項邊界積分可由邊界條件 代換，而得

( )

[

∇ ⋅∇ +

]

−

∫

−

∫

(

− +

)

=0

∫

Ω c u v auvdx Ω fvdx ∂Ω qu g vds v∀ （2.10） 上式稱為 variational 或 weak 式型的微分方 程式。很顯然地，任何微分方程式的解亦為 變分方程式的解，反過來說亦成立，只要定 義域與方程式係數在適當限制之下。（2.10） 式的解 u 及測試函數 v 屬於某一函數空間 V 。有限元素法下一步則是選擇一個函數集 p N ，使得次空間V V p N ⊂ ，而將（2.10）式 變為

( )

[

∇ ⋅∇ + −

]

−

∫

(

− +

)

=0

∫

Ω c u φi auφi fφi dx ∂Ω qu g φids ， i=1Κ N （2.11）_p 將 u 以φ 為基礎展開即_p

( )

_∑

( )

= = Np j j j x U x u 1 φ （2.12） j U 方程式為待係數，故（2.11）

(

)

[

]

{

}

∑ ∫

₌ Ω ∇ ⋅∇ + +

∫

∂Ω p N j j i j i j i j a dx q dsU c 1 φ φ φ φ φ φ

∫

∫

Ω + ∂Ω = fφ_idx gφ_ids （2.13） 若定義

(

)

∫

Ω ∇ ⋅∇ = c dx K_i,_j φ_j φ_i （2.14）

∫

Ω = a dx M_i,_j φ_jφ_i （2.15）

∫

∂Ω = q ds Q_i,_j φ_jφ_i （2.16）

∫

Ω = f dx F_i φ_i （2.17）

ds f G_i

∫

_i Ω ∂ = φ （2.18） 則（2.13）式可表為較簡捷方式

(

K+M +Q

)

U =F +G （2.19） 上式中 K 、 M 、Q 為Np×Np矩陣， F 與 G 為N 向量，通常將 K 、 M 、 Q 合並為 K ，_p 而 F 、G 簡單表為 F ，即橢圓形偏微分方程 式可簡單寫成KU = F。

Matlab PDE toolbox

商用軟體 Matlab 之 PDE（Partial Differential Equation）工具組提供了在處 理二維偏微分方程式上方便有利的解決之 道，PDE toolbox 包函三大部份 定義偏微分方程式，即定義問題區域、邊界 條件及偏微分方程式係數。 解微分問題，即網格產生、離散化問題與求 得近似數值解。 輸出結果，包括各結點值、網格資料，繪出 二維或三維結果。 同時 PDE toolbox 也提供使用者圖形介面， 使得在整個問題定義及處理上更加方便。當 然此一工具組也遵循開放架構理念，也就是 說，使用者可依情況自定函式或呼叫已有的 程式，以達到可擴充、易整合之目的。 此工具組可適用於 time-dependence、 特徵值及非線性等各類偏微分方程式問題。 在使用上，除了有完整手冊、線上說明之外， 亦可在網路論壇或 Matlab 公司網站求得必 要幫助。 程式概述

在使用 Matlab PDE toolbox 時，可分成兩種 模式執行，一為視窗模式，使用者依計算目 的點選所須項目，在適當時機再輸入必要參 數，例如方程式係數、網格產生參數、重要 參考點等，最後點選執行，便可得到所須結 果，但因視窗模式下，並非所有內定函式均 能使用，同時在呼叫自有函數上的不便，加 上本文必須執行諸多算例，故在第三章所有 算例均在第二種模式下執行，即命令列模 式。 本文所使用程式執行依程序大約可分成五個 步驟： 問題幾何外形定義 邊界條件 方程式定義 求解網格點壓力 繪出與儲存結果 在附錄中列了 main.m、maing.m、mainb.m、 mainf.m 及 mainpick.m 五個程式，當中 main.m 為主程式，亦即在 Matlab 命令視窗 中鍵入 main 便可執行，該程式呼叫了 maing.m、mainb.m 及 mainf.m，定義統御方 程式係數，產生網格並求解，最後繪出結果 及儲存結果。maing.m 為問題幾何定義檔， 當中定義了楔形三個頂點，網格產生即利用 此三個頂點做線性內插而得。mainf.m 定義 點聲源位址。mainfpick.m 定義二維點聲源 產生方式，以點聲源所在元素（element）面 積倒數為其大小，亦即必須盡量使點聲源所 屬元素面積愈小，就能使聲源更接近似點聲 源。mainb.m 說明邊界條件，限定了程式是 否要以線性或非線性方式解題。另外，在各 個程式中附加必要說明，以使程式更清楚易 解。 測試項目摘要 本文將在第三章對楔形聲傳特性做適當討 論，藉著對參數改變以了解各參數對聲場影 響程度。計算內容包括： 理想楔形聲場：驗証程式正確性，經由比較 解析解與程式結果，以期驗証程式是否正 確，並觀查兩者之間差異點，做為深入探討 的基石。 幅射邊界條件檢驗：邊界條件建置合理性， 對幅射邊界條件唯一能知道的原則為邊界必 須距聲源足夠遠，但基於程式計算時間及效 益考量，必先對此一邊界進行測試，以期定 出合理而適當的計算區域。 聲源深度的影響：聲源對陰影之影響，將聲

源放置於不同深度，於聲場中合理量測聲壓 值，以陰影區變化情形探討聲源對聲場的探 討。 頻率變化的影響：將問題固定只改變頻率條 件下，對低頻聲源與陰影區間之關係，做一 特性上的探討。 楔頂角大小的影響：同樣地，在固定所有參 數的情形下，變　　動楔頂角大小，觀察楔 頂角與陰影區間的關係。 結論與展望 本文經由對平面海洋楔形聲場影響陰影 區大小的各個參數做測試，雖然在海床底層 設為剛性邊界不是完全合乎實際海洋狀況， 但做為三維楔形海洋聲場傳播特性研究上， 應可做為基本的參考標的。 在參數測驗上共分為四大類： 1. 幅射邊界條件：此一類參數主要目的在了 解楔形聲場邊界條件設立上是否恰當，經 由結果顯示，不論聲源是位在邊界上或楔 形內，只要聲源不要太接近楔形點或位在 楔形內且聲源深度不為波長整數倍，聲源 不一定要離幅射邊界”足夠”遠的定義。 2. 聲源與水面距離：當聲源距水面距離逐 漸由淺而深時，陰影區將逐漸變小，且聲 源在距水面某一距離內，由於某些模態會 被截斷，而使得陰影區會突然增大。 3. 頻率變化：經聲源變化陰影區大小大致 可依頻率高低分成兩大族群，在頻率低於 150Hz 時，陰影區將隨頻率提高而成反向 變化，而在 150Hz 之後，陰影區與頻率成 正向變動，當頻率高於 300Hz 後，聲波即 不再向外傳播。 4. 楔頂角大小：此一測試項目所得結果最 易被了解，即當楔頂角變大時，陰影區將 逐漸變小。 由於本文所使用工具只適合計算平面問題， 在往後研究上勢必得使用其它方法，方能以 數值方法模擬更複雜海洋環境，同時也必須 對海洋環境做一些實際量測，諸如聲速剖 面、海床特質與海床外觀等，這些工作勢必 花費極大心血方能完成，另外，在程式結果 驗証上至今在實際海洋量測也存在著諸多困 難等待解決。 參考文獻

[1] A. Somerfeld, “Mathematische Theorie der difraktion,” Math. Ann., 47, 1896 [2] Michael J. Buckingham and Alexandra Tolstoy, “An analytical solutionfor benchmark problem 1: The “ideal” wedge,” J. Acoust. Soc. Am., vol.87, pp1511-1513,1990.

[3] Stewart A. L. Glegg and Jong R. Yoon , “Experimental measurements of

three-dimensional propagation in a wedge-shaped with pressure-release boundary conditions,” J. Acost. Soc. Am., vol.87(1), pp101-105

[4] F. D. Tappert, “The Parabolic Approximation Method,” in Wave Propagation and Underwater Acoustics, edited by J. B. Keller and J. S. Papadakis, Lecture Notes in Physics 70 (Springer-Verlag, New York, 1977) [5] Finn B. Jensen and Carlo M. Ferla,

“Numerical solutions of range-dependent benchmark problems in ocean acoustics,” J. Acost. Soc. Am. vol.87(4), pp1499-1510

[6] Matlab PDE toolbox User’s Guide, The Mathworks, Inc., 1995

成果自評

目標、研究成果兼具學術及應用價值、適合 在學術期刊發表。