波動(dòng)數(shù)值模擬中的外推型人工邊界條件1)

邢浩潔 李小軍 劉愛(ài)文 李鴻晶 周正華 陳 蘇

?(中國(guó)地震局地球物理研究所,北京 100081)

?(北京工業(yè)大學(xué)建筑工程學(xué)院,北京 100124)

??(南京工業(yè)大學(xué)土木工程學(xué)院,南京 211816)

??(南京工業(yè)大學(xué)交通運(yùn)輸工程學(xué)院,南京 211816)

引言

波動(dòng)數(shù)值模擬是當(dāng)前地球物理、地震學(xué)、地震工程、海洋聲學(xué)、電磁波等多個(gè)學(xué)科領(lǐng)域廣泛應(yīng)用和研究的前沿?zé)狳c(diǎn)技術(shù).人工邊界條件為該技術(shù)的基礎(chǔ)問(wèn)題之一,它表示有限計(jì)算模型邊界上的輻射條件,其精度和穩(wěn)定性對(duì)于模擬結(jié)果能否很好地反映實(shí)際無(wú)限介質(zhì)中的波動(dòng)過(guò)程至關(guān)重要.

自20 世紀(jì)60 年代以來(lái),已有大量人工邊界條件(artificial boundary condition,ABC)被學(xué)者和工程師們所提出,文獻(xiàn)[1-2] 對(duì)它們作了較為系統(tǒng)的介紹.對(duì)于數(shù)量眾多的ABC,如何有效掌握并在不同問(wèn)題中進(jìn)行選用已成為一個(gè)突出的實(shí)際問(wèn)題.一直以來(lái)關(guān)于各個(gè)特定ABC 的應(yīng)用及探討較多,而從整體上針對(duì)ABC 理論與公式體系的系統(tǒng)性研究[3]則相對(duì)較少.這不僅導(dǎo)致較為復(fù)雜的新近研究成果不易被理解及吸收,甚至對(duì)于一些簡(jiǎn)單且有著廣泛應(yīng)用的經(jīng)典人工邊界條件,也仍然存在認(rèn)識(shí)方面的不足.

廖氏透射邊界[4-6]是中國(guó)學(xué)者在人工邊界領(lǐng)域做出的重要成果.它采用等距離散點(diǎn)的外推公式(Newton-Gregory 公式) 來(lái)推算任意人工邊界節(jié)點(diǎn)在各個(gè)時(shí)刻的運(yùn)動(dòng).廖振鵬等將該邊界條件與“誤差波多次透射”的物理機(jī)制聯(lián)系起來(lái),據(jù)此命名為“多次透射公式”(multi-transmitting formula,MTF),也可稱為廖氏透射邊界.MTF 以簡(jiǎn)潔的離散表達(dá)式來(lái)描述任意外行波的一般傳播過(guò)程,完美地將嚴(yán)格的數(shù)學(xué)公式和清楚的物理機(jī)制融合在一起,其基本思想和表達(dá)形式均具有普適性.在廖振鵬等的工作基礎(chǔ)上,作者經(jīng)過(guò)理論分析和波動(dòng)模擬實(shí)踐發(fā)現(xiàn),MTF 邊界的上述特征實(shí)際上奠定了一大類人工邊界條件的思想和公式基礎(chǔ).本文將這類邊界稱為外推型人工邊界條件,其特點(diǎn)是人工邊界節(jié)點(diǎn)在每一時(shí)刻的運(yùn)動(dòng),直接由一組臨近離散網(wǎng)格點(diǎn)在前若干時(shí)刻的運(yùn)動(dòng)進(jìn)行時(shí)空外推得到.具有時(shí)空外推特點(diǎn)的人工邊界條件除了MTF 之外,主要包括經(jīng)典的旁軸近似邊界[7-8]、Higdon 邊界[9-10]、Givoli-Neta 等[11-12]的輔助變量高階邊界、Hagstrom 等[13-14]的輔助變量高階邊界、Guddati 等[15-16]的任意大角度波動(dòng)方程(arbitrarily wide-angle wave equations,AWWE)、Lindman 數(shù)值優(yōu)化邊界[17]、Peng-Toks¨oz 數(shù)值優(yōu)化邊界[18],以及Randall 勢(shì)分解方法[19]、Liu-Sen 混合邊界方法[20-21]等.這些邊界使用相同的控制參數(shù)(邊界階次和一組計(jì)算波速),通過(guò)某種等價(jià)的中間形式相聯(lián)系,并且在數(shù)值模擬中表現(xiàn)出相似的精度和穩(wěn)定性.

現(xiàn)有研究雖然早已指出 MTF 邊界和經(jīng)典的Clayton-Engquist 旁軸近似邊界以及Higdon 邊界存在密切聯(lián)系[9,22],但是在討論和應(yīng)用中仍然將它們視作不同的人工邊界條件.由于上述邊界至今尚未被歸并到一個(gè)理論框架當(dāng)中,導(dǎo)致對(duì)于它們的理論聯(lián)系以及與精度和穩(wěn)定性相關(guān)的諸多應(yīng)用問(wèn)題還缺乏系統(tǒng)性認(rèn)識(shí),如:表達(dá)形式與波動(dòng)類型、內(nèi)域離散格式、幾何構(gòu)型的無(wú)關(guān)性,邊界階次和計(jì)算波速參數(shù)對(duì)邊界精度的控制原理,獨(dú)立的邊界離散格式所導(dǎo)致的穩(wěn)定性問(wèn)題等.這不利于該類人工邊界條件的研究成果的相互借鑒和進(jìn)一步發(fā)展與應(yīng)用.另外值得指出的是,真正與上述邊界在表達(dá)形式、數(shù)值離散格式、精度控制原理以及穩(wěn)定性各個(gè)方面都存在顯著不同的是另外兩大類人工邊界條件,分別為以黏性邊界、黏彈性邊界等為代表的應(yīng)力型邊界[2,23-24]和以函數(shù)衰減層、完美匹配層等為代表的衰減層型邊界[1,25-26].

本文工作以新提出的兩個(gè)基本邊界公式作為出發(fā)點(diǎn),通過(guò)公式推導(dǎo)證明它們的等價(jià)性并論證其與具有時(shí)空外推特點(diǎn)的上述各種ABC 的理論聯(lián)系,建立外推型人工邊界條件理論.深入研究外推型ABC的精度控制原理,分析并強(qiáng)調(diào)可調(diào)的計(jì)算波速參數(shù)對(duì)于處理復(fù)雜波動(dòng)問(wèn)題的價(jià)值.最后給出基本邊界公式的數(shù)值離散格式,并通過(guò)算例驗(yàn)證外推型ABC理論的合理性.

1 外推型人工邊界條件的基本公式

1.1 兩個(gè)基本邊界公式

本文提出兩個(gè)基本公式作為應(yīng)用和討論外推型人工邊界條件的出發(fā)點(diǎn),即

這兩個(gè)公式本身是最簡(jiǎn)單有效的外推型ABC,它們奠定了以時(shí)間-空間外推(簡(jiǎn)稱時(shí)空外推)方式計(jì)算任意時(shí)刻單個(gè)人工邊界節(jié)點(diǎn)運(yùn)動(dòng)的理論和方法基礎(chǔ).其他外推型ABC 大多可以從這兩個(gè)基本邊界公式轉(zhuǎn)化得到,或者通過(guò)某種等價(jià)的中間形式與之相關(guān)聯(lián).式(1)為MTF 引入多個(gè)人工波速caj(j=1,2,···,N,這里N為邊界階次)得到的優(yōu)化形式,本文記為優(yōu)化的MTF 或caj-MTF 邊界.它為離散表達(dá)式,人工邊界節(jié)點(diǎn)的波場(chǎng)值由一系列參考點(diǎn)的波場(chǎng)值外推得到.式(2) 為Higdon 邊界在一個(gè)統(tǒng)一局部坐標(biāo)s0t上定義并使用人工波速caj作為邊界參數(shù)的一般形式,記為統(tǒng)一的Higdon 或caj-Higdon 邊界.它為連續(xù)表達(dá)式,是由多個(gè)一階單向波動(dòng)算子相乘形成的波動(dòng)微分方程.

式(1) 和式(2) 定義在如圖1 所示的以各個(gè)待求解的人工邊界節(jié)點(diǎn)為原點(diǎn)的統(tǒng)一局部坐標(biāo)s0t中,坐標(biāo)空間軸s位于經(jīng)過(guò)人工邊界節(jié)點(diǎn)的一條指向內(nèi)域的離散網(wǎng)格線上(正方向向內(nèi)),時(shí)間軸t與整體模型一致.式中u=u(s,t)為局部坐標(biāo)s0t中的波場(chǎng)函數(shù);N為邊界階次,在式(1) 中對(duì)應(yīng)于離散時(shí)空外推算子(I-B(caj))的個(gè)數(shù),在式(2)中為集成的一階單向波動(dòng)微分算子(?/?t-caj?/?s)的個(gè)數(shù);caj(j=1,2,···,N) 為一組計(jì)算波速參數(shù)稱之為人工波速,其取值具有任意性,但較優(yōu)的選擇通常為等于或稍大于介質(zhì)物理波速的值(后文將具體探討);Δt為時(shí)間步長(zhǎng),取值同樣具有任意性,不過(guò)為便于計(jì)算,通常取為與內(nèi)域格式一致;I為單位算子,在乘積中可以略去.在局部坐標(biāo)s0t中,離散算子B(caj) 定義為:B(caj)u(s,t)=u(s+cajΔt,t-Δt),即空間上向計(jì)算區(qū)域內(nèi)部、時(shí)間上向以前時(shí)刻移動(dòng)(cajΔt,-Δt),它用于定出caj-MTF 邊界的各個(gè)參考點(diǎn)的位置(相對(duì)于人工邊界節(jié)點(diǎn)).離散算子的求和、乘積運(yùn)算法則為:[B(ca1)+B(ca2)]u(s,t)=u(s+ca1Δt,t-Δt)+u(s+ca2Δt,t-Δt),B(ca1)B(ca2)u(s,t)=u(s+ca1Δt+ca2Δt,t-2Δt),依次類推.邊界表達(dá)式(1)和式(2)可直接用于處理單分量波動(dòng)問(wèn)題,如聲波、SH 波動(dòng),也可以分別應(yīng)用于多分量波動(dòng)問(wèn)題的各個(gè)分量,如彈性波、電磁波、固液兩項(xiàng)介質(zhì)波動(dòng)等.

圖1 用于定義基本邊界公式(1)和(2)的統(tǒng)一局部坐標(biāo)示意圖Fig.1 Sketch map of the unified local coordinate used to define the basic boundary formulas(1)and(2)

圖2 給出了二階、三階廖氏透射邊界(MTF)和本文外推型邊界公式(1) (caj-MTF) 所涉及的離散參考點(diǎn)示意圖.基于單一人工波速的廖氏透射邊界是邊界(1) 的一個(gè)特例.利用上述離散算子廖氏透射邊界可以表示為(I-B(ca))Nu=0.如圖2(a)和圖2(c) 所示,MTF 的參考點(diǎn)為等距離分布且在每個(gè)t-jΔt時(shí)刻只有一個(gè)點(diǎn),各個(gè)參考點(diǎn)的算子表示依次為:u(s+caΔt,t-Δt)=B(ca)u(s,t),u(s+2caΔt,t-2Δt)=B2(ca)u(s,t),u(s+3caΔt,t-3Δt)=B3(ca)u(s,t),其余類推.二階MTF 采用前兩個(gè)參考點(diǎn),其系數(shù)分別為2 和-1;三階MTF 采用前3 個(gè)參考點(diǎn),其系數(shù)分別為3,-3 和1.而從圖2(b) 可以看到,二階caj-MTF 與二階MTF 相比,2B(ca)u(s,t)被優(yōu)化成[B(ca1) +B(ca2)]u(s,t),-B2(ca)u(s,t) 被優(yōu)化成-B(ca1)B(ca2)u(s,t).類似地,根據(jù)圖2(d),三階caj-MTF 與三階MTF 相比,3B(ca)u(s,t) 被優(yōu)化成[B(ca1)+B(ca2)+B(ca3)]u(s,t),-3B2(ca)u(s,t) 被優(yōu)化成-[B(ca1)B(ca2)+B(ca2)B(ca3)+B(ca3)B(ca1)]u(s,t),B3(ca)u(s,t) 被優(yōu)化成B(ca1)B(ca2)B(ca3)u(s,t).因此,caj-MTF 邊界在各個(gè)時(shí)刻可能有一個(gè)或多個(gè)參考點(diǎn),它對(duì)傳統(tǒng)MTF 的優(yōu)化是通過(guò)在各個(gè)離散時(shí)刻上由多個(gè)不同參考點(diǎn)的波場(chǎng)值代替原來(lái)單個(gè)參考點(diǎn)的波場(chǎng)值而實(shí)現(xiàn)的.如果從MTF“誤差波多次透射”的物理機(jī)制角度來(lái)分析caj-MTF,那么后者可以被看作是保留了“多次透射”過(guò)程,但優(yōu)化了對(duì)各階“誤差波”的描述.

圖2 二階、三階MTF 和caj-MTF 邊界的參考點(diǎn)示意圖Fig.2 Sketch map of the reference points in 2nd-and 3rd-order MTF and caj-MTF boundaries

Higdon 邊界對(duì)于聲波和彈性波分別具有不同形式,它們基于整體坐標(biāo)x,y或z給出,且根據(jù)外行波沿坐標(biāo)的正、負(fù)方向傳播而有所不同.式(2)可以視為這些現(xiàn)有Higdon 邊界公式的統(tǒng)一形式.理由如下:考慮到邊界參數(shù)取值的任意性,現(xiàn)有Higdon 邊界的聲波和彈性波形式實(shí)際上是一致的[9-10];用統(tǒng)一局部坐標(biāo)的空間軸s代替整體坐標(biāo)x,y或z,在推算人工邊界節(jié)點(diǎn)運(yùn)動(dòng)時(shí)結(jié)果并無(wú)差異;同樣,統(tǒng)一考慮s軸指向內(nèi)域,與考慮s軸指向外域亦具有相同效果,Higdon[27]對(duì)此已有結(jié)論.用統(tǒng)一局部坐標(biāo)的空間軸s代替整體坐標(biāo)x,y或z還有一個(gè)優(yōu)勢(shì),就是s軸可以不與整體坐標(biāo)軸平行.式(2) 使得現(xiàn)有Higdon 邊界針對(duì)不同方向邊界分別給出表達(dá)式再進(jìn)行數(shù)值離散的做法得到大大簡(jiǎn)化,它略去了不必要的冗余,回歸到時(shí)空外推的本質(zhì).

式(1) 給出的多人工波速離散邊界表達(dá)式無(wú)法單獨(dú)從傳統(tǒng)廖氏透射邊界理論導(dǎo)出,其借鑒了Higdon 邊界的計(jì)算波速設(shè)置.式(2)給出的簡(jiǎn)潔而統(tǒng)一的連續(xù)邊界表達(dá)式同樣無(wú)法單獨(dú)從現(xiàn)有Higdon 邊界導(dǎo)出,其得益于廖氏透射邊界思想和參數(shù)的引入.Liao[28]曾指出過(guò)廖氏透射邊界與Higdon 邊界的這種互補(bǔ)性.

1.2 添加小量修正的邊界公式

在基本邊界公式(1)和式(2)中添加一個(gè)小量對(duì)其進(jìn)行修正所得到的形式,在消除高階外推型ABC的計(jì)算漂移失穩(wěn)[29]和處理強(qiáng)衰減外行波方面具有重要價(jià)值.添加小量修正的基本邊界公式為

式(3) 中的小量γ 為一個(gè)小正數(shù),如0.001～0.05 之間.式(4)中的小量εj(j=1,2,...,N)為非負(fù)數(shù),其中至少一個(gè)不為零.

周正華和廖振鵬[30]指出添加小量在物理上考慮了介質(zhì)的幾何擴(kuò)散特性或引入了介質(zhì)的阻尼特性.可認(rèn)為這兩種物理特性對(duì)應(yīng)于同一種數(shù)學(xué)機(jī)制,即考慮了波場(chǎng)的衰減性.由此不難得出,陳少林和廖振鵬[31]所提出的衰減波場(chǎng)MTF、以及Bayliss 和Turkel[32]針對(duì)柱面波或球面波擴(kuò)散問(wèn)題(即外行波含有顯著的幾何衰減)而提出的邊界表達(dá)式,分別與式(3)、式(4)在形式上具有一致性.

2 caj-MTF 與caj-Higdon 邊界的等價(jià)性

本節(jié)證明外推型ABC 理論與公式體系中的一個(gè)關(guān)鍵結(jié)論:離散的caj-MTF 邊界與連續(xù)的caj-Higdon邊界是等價(jià)的.有兩種證明方法:一是利用二元函數(shù)泰勒展開(kāi)法則對(duì)caj-MTF 中各個(gè)參考點(diǎn)的波場(chǎng)函數(shù)在點(diǎn)(s,t) 處進(jìn)行展開(kāi)并保留至所需精度階,化簡(jiǎn)可得caj-Higdon;二是將caj-MTF 中的各個(gè)外推算子改寫(xiě)成某種與之等效的差分算子,使之形成類似于caj-Higdon 的結(jié)構(gòu),對(duì)差分算子取極限即得到后者.

證明 1:對(duì)于一階caj-MTF,其具體形式為u(s,t)=u(s+ca1Δt,t-Δt).將等號(hào)右端項(xiàng)進(jìn)行二元泰勒展開(kāi)并保留至一階項(xiàng),得到

化簡(jiǎn)式(5)即得一階caj-Higdon 表達(dá)式.對(duì)于二階caj-MTF,先寫(xiě)出其具體形式,再將等號(hào)右端各項(xiàng)進(jìn)行二元泰勒展開(kāi)并保留至二階項(xiàng),得到

化簡(jiǎn)式(6)就能夠得到二階caj-Higdon 表達(dá)式.對(duì)于更高階次邊界,可用歸納法依次類推.

證明2:考慮caj-MTF 的單個(gè)算子形式(I-B(caj))u=0.分別定義空間移動(dòng)算子K和時(shí)間移動(dòng)算子Z-1:Ku(s,t)=u(s+cajΔt,t),Z-1u(s,t)=u(s,t-Δt),則式(1)的單個(gè)算子形式可以表示為(I-Z-1K)u=0.該式可改寫(xiě)為

將式(7) 中第一個(gè)Δs用cajΔt替換然后對(duì)每一項(xiàng)乘以caj,得到

對(duì)式(8) 中的差分算子取極限,即可得到(?/?t-?/?s)u=0.這表明caj-MTF 中每個(gè)離散形式的時(shí)空外推算子的極限形式為caj-Higdon 的各個(gè)單向波動(dòng)微分算子,所以這兩個(gè)邊界表達(dá)式等價(jià).

由于證明1 和證明2 都建立在基本的數(shù)學(xué)原理之上,證明過(guò)程中并未涉及任何波動(dòng)方程,所以caj-MTF 與caj-Higdon 邊界的等價(jià)性是普遍存在的.這遠(yuǎn)遠(yuǎn)強(qiáng)于已有研究所指出的Higdon 邊界或MTF 邊界與Clayton-Engquist 聲波邊界之間的相似性[9,22].這種等價(jià)性將離散caj-MTF 邊界與連續(xù)caj-Higdon 邊界緊密聯(lián)系在一起,在討論和應(yīng)用中二者通?？梢曰ハ嗵鎿Q.

3 外推型ABC 的理論聯(lián)系

這一部分將通過(guò)公式推導(dǎo)來(lái)論證其他外推型ABC 大多可以從基本邊界公式(1) 和(2) 轉(zhuǎn)化得到,或通過(guò)某種等價(jià)的中間形式與之相關(guān)聯(lián).

3.1 Reynolds 邊界、Keys 邊界

Reynolds[33]和Keys[34]各自獨(dú)立提出了一種ABC.以外法向?yàn)閤軸負(fù)向的邊界為例,其表達(dá)式分別為

將式(9)、式(10) 與式(2) 相比較可知,它們都是式(2)的特例.因此,caj-Higdon 邊界已將這兩種邊界概括在內(nèi).

3.2 旁軸近似邊界

旁軸近似邊界主要為經(jīng)典的 Clayton-Engquist聲波和彈性波邊界(CE 邊界)[7]、Halpern-Trefethen 大角度單向波動(dòng)方程[8]、以及Fuyuki-Matsumoto(FM)[35]、Emerman-Stephen (ES)[36]、Stacey[37]對(duì)CE彈性波邊界的改進(jìn).

對(duì)于聲波旁軸近似邊界,Higdon[9]已證明它們可以由Higdon 邊界結(jié)合聲波方程來(lái)導(dǎo)出.以外法向?yàn)閦軸負(fù)向的邊界為例,考慮二階情形,CE 和Halpern-Trefethen 所提出的各種二階聲波邊界表達(dá)形式為

其中c為聲波波速,p0和p2為不同邊界的控制參數(shù),具體見(jiàn)文獻(xiàn)[8].容易驗(yàn)證,在式(2)的二階形式中引入聲波方程c2(?2u/?x2+?2u/?z2)=?2u/?t2,將關(guān)于z的二階偏導(dǎo)用關(guān)于x和t的二階偏導(dǎo)所替換,就導(dǎo)出了式(11).式(2) 中參數(shù)ca1,ca2與式(11) 中參數(shù)p0,p2存在如下對(duì)應(yīng)關(guān)系:p0=(1+α1α2)/(α1+α2),p2=α1α2/(α1+α2),其中α1=ca1/c,α2=ca2/c.

對(duì)于彈性波旁軸近似邊界,作者分析發(fā)現(xiàn)它們?nèi)狈?yán)格的理論基礎(chǔ),難以達(dá)到較好的精度.由于彈性波方程不能由其頻散關(guān)系唯一地確定,旁軸近似技術(shù)實(shí)際上無(wú)法用于推導(dǎo)彈性波邊界.Clayton 和Engquist[7]僅通過(guò)粗略地模仿聲波邊界給出了如下彈性波邊界(z軸負(fù)向)

然而,式(12)既無(wú)法從旁軸近似角度進(jìn)行解釋,也無(wú)法由caj-Higdon 邊界結(jié)合彈性波方程導(dǎo)出.Fuyuki 和Matsumoto[35]、Emerman 和Stephen[36]、Stacey[37]在式(12) 基礎(chǔ)上所做的改進(jìn)并未解決這一問(wèn)題.數(shù)值試驗(yàn)表明這些邊界精度較差.為了從理論上揭示所謂的旁軸近似彈性波邊界的不合理性,在波數(shù)平面上繪出了CE,FM,ES 和Stacey 彈性波邊界的頻散曲線,如圖3 所示.作為對(duì)比,圖中還給出了聲波邊界(11)的頻散曲線,A1,A2,A3 分別表示一階、二階、三階聲波邊界.

圖3 中各曲線與上半圓的接近程度反映了人工邊界條件的精度.各個(gè)彈性波邊界頻散曲線與上半圓的吻合度很差,這證實(shí)了旁軸近似彈性波邊界的不合理性.與之相比,旁軸近似聲波邊界精度較好,且隨著階次上升迅速提高,與理論結(jié)果相符.

圖3 彈性波、聲波旁軸近似邊界的頻散曲線Fig.3 Dispersion curves of acoustic-and elastic-wave paraxial-approximation boundaries

3.3 輔助變量高階邊界

針對(duì)高階Higdon 邊界數(shù)值離散比較困難的問(wèn)題,Givoli-Neta 等[11-12]、Hagstrom-Warburton 等[13-14]分別提出一種可順利離散至任意階次的輔助變量高階邊界(G-N 邊界和H-W 邊界).這兩種邊界都是先給出一個(gè)與Higdon 邊界等價(jià)的輔助變量遞推關(guān)系,再引入聲波方程將遞推關(guān)系中的邊界法向?qū)?shù)轉(zhuǎn)換為切向?qū)?shù),得到只含有邊界切向?qū)?shù)和時(shí)間導(dǎo)數(shù)的輔助變量遞推關(guān)系,形成實(shí)用的高階邊界.由于聲波方程的引入,這兩種為聲波邊界.

G-N 邊界表達(dá)式見(jiàn)文獻(xiàn)[11-12],它所使用的輔助變量遞推關(guān)系為

逐步消去輔助變量,可得式(13)的合并形式為

將式(14) 與式(2) 進(jìn)行比較可知它們是一致的.式(14)中的計(jì)算波速參數(shù)c1,c2,···,cN對(duì)應(yīng)于式(2)中的人工波速參數(shù)ca1,ca2,···,caN.所以在聲波問(wèn)題中,相同參數(shù)的G-N 邊界與caj-Higdon 邊界具有相同的精度.

H-W 邊界表達(dá)式見(jiàn)文獻(xiàn)[13-14],它基于的輔助變量遞推關(guān)系為

B′ecache 等[38]經(jīng)過(guò)推導(dǎo),證明式(15) 的合并形式為

式 (16) 中計(jì)算波速由聲波波速c和任意系數(shù)a0,a1,a2,···,aN組合而成,為c/a0,c/a1,c/a2,···,c/aN.比較式(16) 和式(2) 可知,在聲波問(wèn)題中,1,2,3,··· 階H-W 邊界分別對(duì)應(yīng)于3,5,7,··· 階caj-Higdon 邊界,即存在2N+1 倍關(guān)系.

3.4 任意大角度波動(dòng)方程

Guddati 等[15-16]提出一種矩陣形式的輔助變量高階邊界AWWE,他們通過(guò)兩種方式導(dǎo)出這種邊界:一是將聲波方程的旁軸近似遞推關(guān)系按偶數(shù)階提取出來(lái),將其寫(xiě)成矩陣形式,再返回到時(shí)域成為人工邊界條件;二是將半無(wú)限介質(zhì)分層,利用各層動(dòng)力剛度的有限單元近似,在頻域內(nèi)導(dǎo)出矩陣關(guān)系式,進(jìn)而返回時(shí)域得到人工邊界條件.根據(jù)本文第3.2 節(jié)所指出的聲波旁軸近似邊界可以由caj-Higdon 邊界與聲波方程相結(jié)合得到,可以推斷出1,2,3,··· 階AWWE邊界分別對(duì)應(yīng)于2,4,6,··· 階caj-Higdon 邊界,即存在2N倍關(guān)系.AWWE 邊界與caj-Higdon 邊界的等價(jià)中間形式為

3.5 其他相關(guān)邊界方法

Lindman[17]提出一種數(shù)值優(yōu)化邊界,表達(dá)式見(jiàn)文獻(xiàn)[17].Kausel[39]證明了該邊界的連續(xù)形式為一種高階旁軸近似邊界,與Clayton-Engquist 聲波邊界屬于一個(gè)類型.

Peng-Toks¨oz[18]提出另一種數(shù)值優(yōu)化邊界,表達(dá)式為

式中un(i,j,k)=u(iΔx,jΔy,kΔz,nΔt)表示離散網(wǎng)格中波場(chǎng)的節(jié)點(diǎn)值.編號(hào)0 為邊界節(jié)點(diǎn),1 和2 為臨近的內(nèi)域節(jié)點(diǎn).a01,a02,···,a22為控制參數(shù),由一系列限定條件建立的方程組來(lái)確定.顯然,式(18)借鑒了二階Higdon 邊界數(shù)值離散格式所采用的3×3 節(jié)點(diǎn)組合,只是確定計(jì)算系數(shù)的方式有所不同.

Randall[19]針對(duì)彈性波邊界精度不高的問(wèn)題,提出一種將彈性波場(chǎng)進(jìn)行勢(shì)分解,對(duì)分解出來(lái)的波勢(shì)函數(shù)分別使用高精度的Lindman 聲波邊界,再將結(jié)果合成為彈性波場(chǎng)的方法.彈性波場(chǎng)向勢(shì)函數(shù)的分解與合成分別表示為

式中,u為彈性波位移場(chǎng)向量,Φ和Ψ分別為P 波和S 波勢(shì)函數(shù).?為向量微分算子,cp和cs分別為P 波和S 波波速.Randall 勢(shì)分解方法是提高外推型ABC 模擬彈性波的精度的一個(gè)有益嘗試,但由于勢(shì)分解與合成過(guò)程以及Lindman 邊界所帶來(lái)的雙重復(fù)雜度,作者未見(jiàn)其他學(xué)者應(yīng)用或探討該邊界方法.

Liu 和Sen[20-21]將外推型ABC 與過(guò)渡區(qū)技術(shù)相結(jié)合,發(fā)展出一套混合邊界方法,如圖4 所示.該方法在計(jì)算模型的外邊界附近設(shè)置一定寬度的過(guò)渡區(qū),過(guò)渡區(qū)波場(chǎng)由內(nèi)域離散格式和邊界數(shù)值格式分別計(jì)算的波場(chǎng)進(jìn)行加權(quán)和得到.二者權(quán)系數(shù)之和始終為1,其分配比例沿著過(guò)渡區(qū)寬度的變化如圖中三角形所示.

圖4 Liu-Sen 混合邊界示意圖Fig.4 Schematic of Liu-Sen hybrid ABC

大量波動(dòng)模擬實(shí)踐表明Liu-Sen 混合邊界能夠顯著提高人工邊界的精度,并有效降低外推型ABC可能面臨的數(shù)值失穩(wěn)風(fēng)險(xiǎn).過(guò)渡區(qū)的引入在相當(dāng)程度上阻止了反射誤差向內(nèi)域傳播,并延緩了不穩(wěn)定結(jié)果的積累.在數(shù)值試驗(yàn)中觀察到過(guò)渡區(qū)內(nèi)波場(chǎng)表現(xiàn)出類似完美匹配層內(nèi)波場(chǎng)的快速衰減特征.我們認(rèn)為本文外推型ABC 與Liu-Sen 混合邊界方法相結(jié)合,應(yīng)是解決波動(dòng)數(shù)值模擬中人工邊界問(wèn)題一個(gè)值得重視的方向.

4 外推型ABC 的精度控制原理

4.1 離散邊界公式(1)的時(shí)空外推原理

離散邊界公式(1)的單人工波速形式(即廖氏透射邊界)存在一種明確的幾何解釋,如圖5 所示.外行波沿人工邊界節(jié)點(diǎn)所在的局部坐標(biāo)空間軸s的傳播過(guò)程,可以在sot時(shí)空坐標(biāo)平面上展開(kāi)成一個(gè)二維“靜態(tài)”波場(chǎng).在這個(gè)二維靜態(tài)波場(chǎng)中,沿著由單一人工波速和離散時(shí)間步距所確定的時(shí)空外推方向,可以截取出一維波場(chǎng)曲線g.人工邊界節(jié)點(diǎn)運(yùn)動(dòng)的實(shí)際值由連續(xù)曲線g決定,但在數(shù)值模擬中只能利用曲線g上若干個(gè)參考點(diǎn)的值來(lái)推算其端部人工邊界節(jié)點(diǎn)的值.這個(gè)推算過(guò)程的精度即為MTF 邊界的精度.

從圖5 中可以看出,邊界階次N(即參考點(diǎn)數(shù)目)代表基于N-1 階多項(xiàng)式的外推精度.曲線g的彎曲程度則由外行波復(fù)雜程度和所選取的時(shí)空外推方向(它與人工波速有關(guān)) 共同決定.于是這種幾何解釋可以用一句話概括為:N階MTF 是以基于N-1 階多項(xiàng)式的外推精度對(duì)由人工波速所確定的時(shí)空外推方向上的外行波曲線的一種近似.顯然,提高邊界階次N,或者選取與實(shí)際視波速相接近的人工波速以獲得較為平緩的曲線g,都是提高M(jìn)TF 邊界精度的有效措施.

圖5 MTF 邊界的時(shí)空外推原理Fig.5 Time-space extrapolation of MTF boundary

MTF 邊界這種明確的幾何解釋,加之其簡(jiǎn)潔的數(shù)學(xué)形式和“多次透射”的物理機(jī)制,奠定了外推型人工邊界理論的思想基礎(chǔ).

4.2 連續(xù)邊界公式(2)的反射系數(shù)乘積原理

連續(xù)邊界公式(2) 是Higdon 邊界的統(tǒng)一形式.Higdon 邊界所采用的多個(gè)單向波動(dòng)微分算子乘積形式是從滿足邊界要求的一系列離散表達(dá)式中總結(jié)而來(lái)的[9].這表明連續(xù)邊界公式(2)的理論源頭可以追溯到離散邊界公式(1),即上述MTF 邊界的時(shí)空外推原理.因此,外推型人工邊界條件具有統(tǒng)一的思想基礎(chǔ),它來(lái)源于離散邊界.各種微分方程形式的連續(xù)邊界主要為外推型ABC 提供了豐富的變化形式.

連續(xù)邊界公式(2) 的精度可以利用反射系數(shù)加以分析.該邊界對(duì)于入射角為θ 的外行平面波的反射系數(shù)為

式中c為介質(zhì)物理波速,caj為邊界的人工波速參數(shù).式(21)為多個(gè)因子的乘積,這對(duì)應(yīng)了式(2)的單向波動(dòng)微分算子乘積.圖6 繪出了幾種不同N和caj取值的反射系數(shù)R(θ).

在圖6 中,根據(jù)一階邊界反射系數(shù)R1 可知單個(gè)反射因子的幅值總是小于或等于1,所以高階邊界所具有的多個(gè)反射因子的乘積能夠迅速提高邊界精度.另外,多個(gè)反射因子如果具有不同的零反射角度(由不同人工波速取值決定),能夠在更大透射角范圍內(nèi)提高邊界精度.

圖6 連續(xù)邊界公式的反射系數(shù)乘積原理Fig.6 Product of reflection coefficients of multiple first-order wave-propagation operators

圖7 繪出了caj-Higdon 邊界(Hig)、Hagstrom-Warburton 邊界(HW)、AWWE 邊界的反射系數(shù).根據(jù)第3.3 節(jié)和第3.4 節(jié)內(nèi)容可知,Givoli-Neta 邊界(GN)反射系數(shù)與caj-Higdon 邊界相同,故未在圖中列出;N階HW 邊界、AWWE 邊界分別對(duì)應(yīng)于2N+1 階、2N階caj-Higdon 邊界.

圖7 曲線表明相同人工波速取值下,1,2,3 階HW 反射系數(shù)與3,5,7 階Higdon 相同,1,2,3 階AWWE 反射系數(shù)與2,4,6 階Higdon 相同.這證實(shí)了前文理論分析所給出的結(jié)論.由此可見(jiàn),GN、HW和AWWE 3 種輔助變量高階邊界使用與caj-Higdon相同數(shù)量的一階單向波動(dòng)微分算子(?/?t-caj?/?s)實(shí)現(xiàn)了與之相同的精度,即它們提高邊界精度的原理與caj-Higdon 并無(wú)二致.

圖7 Higdon,G-N,H-W,AWWE 邊界的反射系數(shù)Fig.7 Reflection coefficients of Higdon,G-N,H-W,AWWE boundaries

4.3 基于視波速的多人工波速分析

基本邊界公式(1)和(2)所采用的多人工波速參數(shù)caj(j=1,2,···,N)并不僅僅是形式上的優(yōu)化,它對(duì)于具有多種物理波速的復(fù)雜波動(dòng)問(wèn)題具有非常重要的實(shí)用價(jià)值.這些問(wèn)題包括多層介質(zhì)問(wèn)題、彈性波、兩相介質(zhì)波動(dòng)、頻散波動(dòng)等.上述基于透射角度θ 的反射系數(shù)式(21)未考慮物理波速的變化.為同時(shí)考慮介質(zhì)物理波速和透射角度的變化,這里采用外行波的視波速cn=c/cos θ 作為基本變量,得出caj-Higdon邊界反射系數(shù)的另一種表達(dá)形式

圖8 所示為存在兩種物理波速時(shí)caj-Higdon 邊界的反射系數(shù)幅值與外行波視波速cn的關(guān)系.這里兩種物理波速為c1=v和c2=3v,v表示某個(gè)可約去的波速,它的取值不影響分析結(jié)果.3 倍的物理波速差異將對(duì)人工邊界性能提出挑戰(zhàn).

在圖8 中,主要波動(dòng)能量的視波速范圍有兩段,其中雖然有透射角度的影響,但差異較大的兩種物理波速起了主要作用.理想的人工邊界反射系數(shù)應(yīng)當(dāng)在所關(guān)心的兩個(gè)區(qū)段中都處于低值.可以看出,采用單一人工波速的3 種邊界R(v,v),R(2v,2v) 和R(3v,3v)都不夠理想,只有采用多人工波速的二階邊界R(v,3v)和三階邊界R(v,3v,3v)較好地兼顧了對(duì)兩段波動(dòng)能量的吸收.因此,在具有多種差異較大的物理波速的波動(dòng)問(wèn)題中,使用以多人工波速為參數(shù)的外推型ABC 很有必要也非常實(shí)用.

圖8 多人工波速分析Fig.8 Analysis of various artificial wave velocities

5 數(shù)值算例

5.1 聲波大角度透射問(wèn)題

以一個(gè)均勻介質(zhì)聲波問(wèn)題算例來(lái)檢驗(yàn)上述外推型ABC 的精度.模型尺度為200 m×500 m,介質(zhì)波速為c=500 m/s.頂面為自由邊界,左、右、底邊界為人工邊界.在頂部中點(diǎn)施加中心頻率為10 Hz 的Ricker 子波激勵(lì).計(jì)算模型的空間、時(shí)間離散均采用標(biāo)準(zhǔn)的二階中心差分格式[9].網(wǎng)格尺寸為2 m×2 m,時(shí)間步長(zhǎng)滿足穩(wěn)定條件cΔt/Δx≤0.71.人工邊界分別采用本文提出的外推型ABC 基本公式caj-MTF 和caj-Higdon,以及現(xiàn)有的Clayton-Engquist (CE) 旁軸近似邊界[7]、輔助變量高階邊界Givoli-Neta(GN)[11]、Hagstrom-Warburton(HW)[13]和AWWE[16].本文ABC 的數(shù)值離散格式見(jiàn)附錄.

圖9 給出了0.9 s 時(shí)的波場(chǎng)快照.此時(shí)外行波在側(cè)邊界上的透射角度約75°,為大角度透射,這種情形對(duì)ABC 性能要求較高.圖中括號(hào)內(nèi)參數(shù)如(c,1.5c,2.5c)表示所采用的人工波速參數(shù)caj(j=1,2,···,N).

圖9 采用不同外推型人工邊界的聲波模擬結(jié)果Fig.9 Modeling results of acoustic wave propagation using different ABCs

圖9 中不同ABC 的反射誤差對(duì)比證實(shí)了前文理論分析所給出的結(jié)論,主要包括:(1)caj-MTF 和caj-Higdon 等價(jià)(見(jiàn)第3 部分),這在本算例中表現(xiàn)為盡管用于實(shí)現(xiàn)caj-MTF 和caj-Higdon 的數(shù)值計(jì)算格式有所不同,但只要當(dāng)二者所采用的邊界參數(shù)N和caj(j=1,2,···,N)相同時(shí),其模擬結(jié)果就會(huì)保持一致;(2)隨著階次N升高,邊界精度提高,且多人工波速方案能夠進(jìn)一步改善精度(見(jiàn)圖6),圖9(a)～圖9(e)中不同邊界階次結(jié)果對(duì)比和相同階次下多人工波速方案與單一人工波速方案結(jié)果對(duì)比證實(shí)了這個(gè)結(jié)論;(3)CE 邊界和輔助變量高階邊界GN,HW 和AWWE 具有與基本邊界公式caj-Higdon 相同的精度控制原理(見(jiàn)第3.2～3.4 節(jié)和圖7),這預(yù)示著相同邊界參數(shù)下它們的模擬結(jié)果一致,如圖9 中二階CE 與二階caj-Higdon在0.9 s 時(shí)的結(jié)果非常接近,二階HW(相當(dāng)于五階caj-Higdon)和三階AWWE(相當(dāng)于六階caj-Higdon)模擬結(jié)果與三階caj-Higdon 結(jié)果相似.不同外推型ABC除了呈現(xiàn)出上述共同特征之外,它們各自還有一些不同的表現(xiàn),如圖中CE 邊界存在角點(diǎn)反射問(wèn)題(見(jiàn)0.5 s 快照?qǐng)D),GN 邊界出現(xiàn)了失穩(wěn)(該失穩(wěn)首先出現(xiàn)在邊界上距離震源較近的區(qū)域,為典型的人工邊界條件失穩(wěn)特征),HW 和AWWE 邊界由于數(shù)值離散誤差影響顯著,其三階以上的高階精度難以實(shí)現(xiàn).這些問(wèn)題從側(cè)面說(shuō)明本文所提出的基本邊界公式(1) 和式(2)最為簡(jiǎn)單實(shí)用.

5.2 均勻介質(zhì)彈性波問(wèn)題

以一個(gè)均勻介質(zhì)彈性波問(wèn)題算例來(lái)進(jìn)一步檢驗(yàn)上述外推型ABC 的性能.模型尺度為2500 m ×2500 m,介質(zhì)縱、橫波速為cp=2000 m/s,cs=1000 m/s.四邊均為人工邊界.在模型中心施加中心頻率為10 Hz 的Ricker 子波激勵(lì).計(jì)算模型的空間、時(shí)間離散均采用標(biāo)準(zhǔn)的二階中心差分格式[35].網(wǎng)格尺寸為5 m×5 m,時(shí)間步長(zhǎng)滿足穩(wěn)定條件Δt/Δx≤1.人工邊界分別采用caj-MTF 和caj-Higdon,以及第3.2 節(jié)提及的幾種彈性波旁軸近似邊界.此時(shí)附錄A 中的邊界離散格式分別用于人工邊界節(jié)點(diǎn)的水平和豎向波場(chǎng)分量.

圖10 給出了0.9 s 和1.5 s 時(shí)水平分量的波場(chǎng)快照.此時(shí)存在兩種差異較大的物理波速(cp等于2倍cs),這種情形能夠凸顯采用多人工波速方案的必要性及優(yōu)勢(shì).caj-Higdon 與caj-MTF 結(jié)果相同,未在圖中列出.考慮了3 種單一人工波速和一種多人工波速的二階caj-MTF 邊界,并以幾種旁軸近似邊界作為對(duì)比.

圖10 中caj-MTF 邊界模擬結(jié)果證明了對(duì)于多物理波速的復(fù)雜波動(dòng)問(wèn)題采用多人工波速外推型ABC的必要性及優(yōu)勢(shì)(第4.3 節(jié)).此時(shí)3 種單一人工波速的二階caj-MTF 均不能很好地實(shí)現(xiàn)對(duì)P 波和S 波的同步透射,而多人工波速(cs,cp) 方案幾乎完美地實(shí)現(xiàn)了這一目標(biāo).與基本邊界caj-MTF 相比,所謂的彈性波旁軸近似邊界反射量很大,這對(duì)應(yīng)了圖3 關(guān)于這幾個(gè)彈性波邊界理論上具有不合理性的分析結(jié)果.Clayton-Engquist 和Fuyuki-Matsumoto(它與前者高度相似) 邊界出現(xiàn)失穩(wěn),Emerman-Stephen 和Stacey 邊界(文獻(xiàn)[36-37] 指出它們改善了CE 彈性波邊界的穩(wěn)定性)以及本文提出的caj-MTF 邊界在20 s 長(zhǎng)持時(shí)模擬結(jié)果中保持了穩(wěn)定.因此,綜合考慮精度與穩(wěn)定性,模擬彈性波的外推型ABC 建議采用基本邊界公式caj-MTF 或caj-Higdon.

圖10 采用不同外推型人工邊界的彈性波模擬結(jié)果Fig.10 Modeling results of elastic wave propagation by using different ABCs

5.3 縱、橫向不均勻介質(zhì)彈性波模擬

設(shè)計(jì)如圖11 所示的縱、橫向不均勻介質(zhì)場(chǎng)地模型,檢驗(yàn)本文邊界條件的吸收效果.該模型為右下方整塊基巖與左側(cè)、上部成層土體形成的復(fù)雜場(chǎng)地,巖土體參數(shù)如圖所示.在(x,z)=(300 m,100 m)處節(jié)點(diǎn)的豎向分量上施加中心頻率為10 Hz 的Ricker 子波激勵(lì).上部為自由地表,其余3 邊分別采用自由邊界(作為對(duì)照) 或本文提出的人工邊界,進(jìn)行兩組模擬.采用集中質(zhì)量有限元與中心差分相結(jié)合的數(shù)值模擬方法[29-30].單元尺寸為2.5 m×2.5 m,時(shí)間步長(zhǎng)滿足穩(wěn)定條件

圖11 縱、橫向不均勻介質(zhì)模型Fig.11 Computational model of a longitudinal and transverse inhomogeneous media

這個(gè)模擬選用所提出的caj-MTF (或caj-Higdon,與前者結(jié)果一致)邊界,邊界參數(shù)為N=2,caj=cs,cp.不考慮其他外推型ABC,因?yàn)榍拔囊阎赋?Givoli-Neta,Hagstrom-Warburton,AWWE 等高階輔助變量邊界和Lindman 數(shù)值優(yōu)化邊界只適用于聲波;算例5.2 表明幾種旁軸近似彈性波邊界的精度和穩(wěn)定性遠(yuǎn)不如caj-MTF 或caj-Higdon;引言中提及的Peng-Toks¨oz 數(shù)值優(yōu)化邊界和Randall 勢(shì)分解邊界十分復(fù)雜,少有應(yīng)用和討論,其是否有精度或穩(wěn)定性方面的優(yōu)勢(shì)尚無(wú)定論;Liu-Sen 混合邊界是外推型ABC 與過(guò)渡區(qū)技術(shù)相結(jié)合的綜合性解決方案,已發(fā)表大量文獻(xiàn),而本文關(guān)注的僅是外推型ABC 本身的性能.模擬結(jié)果如圖12 所示.

在圖12 中,由于自由邊界會(huì)完全反射波動(dòng)能量,該模型像一個(gè)封閉波動(dòng)能量的“盒子”,其模擬結(jié)果與實(shí)際問(wèn)題中波可以無(wú)阻擋地穿過(guò)(因?yàn)樵瓎?wèn)題中此處無(wú)界面)人工邊界再傳向遠(yuǎn)處的過(guò)程相去甚遠(yuǎn).caj-MTF(或caj-Higdon)邊界則幾乎看不到虛假反射波,說(shuō)明它們所推算的邊界節(jié)點(diǎn)運(yùn)動(dòng)與原來(lái)無(wú)限介質(zhì)中該處運(yùn)動(dòng)非常接近(或基本一致),因此這個(gè)有限計(jì)算模型給出了與實(shí)際問(wèn)題相符的模擬結(jié)果.

圖12 不均勻介質(zhì)彈性波模擬結(jié)果Fig.12 Simulation results of elastic waves in the inhomogeneous model

5.4 關(guān)于復(fù)雜介質(zhì)模擬效果和穩(wěn)定性的討論

外推型ABC 模擬復(fù)雜介質(zhì)的效果已被現(xiàn)有文獻(xiàn)給出的各種波動(dòng)問(wèn)題模擬結(jié)果所證實(shí),如文獻(xiàn)[4]混凝土壩開(kāi)裂模擬中的邊界處理,文獻(xiàn)[11] 中頻散介質(zhì)波動(dòng)模擬,文獻(xiàn)[16] 對(duì)力學(xué)性質(zhì)漸變介質(zhì)的反演,文獻(xiàn)[21] 對(duì)三維VTI 介質(zhì)的波動(dòng)過(guò)程正演,文獻(xiàn)[31]中適用于兩相介質(zhì)波動(dòng)問(wèn)題的邊界,文獻(xiàn)[35]對(duì)Rayleigh 波邊界的探討等.在本文建立的外推型人工邊界條件理論與公式體系中,這些復(fù)雜波動(dòng)過(guò)程的絕大多數(shù)特征已得到充分考慮,具體為:(1)多種人工波速的參數(shù)配置能夠很好地考慮由多種物理波速(包括彈性波、兩相介質(zhì)波動(dòng)、頻散波動(dòng)、介質(zhì)交界面等)以及不同透射角度所決定的沿外推型ABC計(jì)算方向視傳播速度的復(fù)雜性;(2)由幾何擴(kuò)散、阻尼效應(yīng)或其他原因引起的波場(chǎng)的衰減性,可以由添加小量修正的邊界表達(dá)式(3) 和式(4) 很好地模擬,或者只是單純地將邊界(1)和邊界(2)的階次提得足夠高,也能夠達(dá)到滿意精度(其原理參考圖5);(3)對(duì)于角點(diǎn)、介質(zhì)交界面、曲線邊界等不同幾何構(gòu)型的適應(yīng)性,caj-MTF 和caj-Higdon 可以無(wú)差別地應(yīng)用,因?yàn)樗鼈冎簧婕暗揭粭l由邊界節(jié)點(diǎn)指向內(nèi)域的離散網(wǎng)格線上的信息,與邊界形狀無(wú)關(guān),而其他外推型ABC若涉及到邊界本身的信息,則會(huì)受到邊界幾何構(gòu)型的制約.除了上述特征之外,在不同復(fù)雜波動(dòng)問(wèn)題中是否還存在其他影響邊界精度的重要因素,如多種波動(dòng)成分之間的耦合效應(yīng)等,將在后續(xù)工作中具體探討.

數(shù)值試驗(yàn)表明,本文基本邊界公式caj-MTF 和caj-Higdon 的穩(wěn)定性高度一致,且與現(xiàn)有MTF 和Higdon 邊界基本沒(méi)有區(qū)別.外推型ABC 的穩(wěn)定性可以從所引文獻(xiàn)及相關(guān)學(xué)者的研究成果中總結(jié)得到,這個(gè)研究課題目前仍然活躍.本文和相關(guān)文獻(xiàn)數(shù)值試驗(yàn)表明了外推型ABC 的實(shí)用性,或者說(shuō)它們一般能夠在所關(guān)心的模擬時(shí)間內(nèi)給出穩(wěn)定結(jié)果,這對(duì)于解決大多數(shù)持時(shí)不太長(zhǎng)的波動(dòng)問(wèn)題已經(jīng)足夠.對(duì)外推型ABC 穩(wěn)定性問(wèn)題完整而詳細(xì)的探討顯然不是本文所能完成的任務(wù),不過(guò),關(guān)于穩(wěn)定性、邊界失穩(wěn)特征及其解決方法,作者從已有研究中總結(jié)出以下幾點(diǎn)主要認(rèn)識(shí):(1)不同于應(yīng)力型ABC(如黏彈性邊界)和衰減層型ABC(如完美匹配層邊界)需要附著在內(nèi)域離散模型之上,外推型ABC 有著獨(dú)立于內(nèi)域離散模型的數(shù)值計(jì)算格式,后者兩種不同離散格式的耦合常常不如前兩者的整體計(jì)算格式穩(wěn)定.(2)當(dāng)邊界出現(xiàn)某種失穩(wěn)時(shí),其不穩(wěn)定值的積累速度由邊界條件的反射幅值(高階邊界大于低階邊界)和波在整個(gè)模型與所有失穩(wěn)節(jié)點(diǎn)之間的反復(fù)反射放大共同決定,因此三維模型(邊界節(jié)點(diǎn)多) 失穩(wěn)放大最快,二維模型次之,一維模型最慢;同樣維度下,小尺寸模型(來(lái)回反射距離短)比大尺寸模型失穩(wěn)積累得更快;高波速問(wèn)題(來(lái)回反射所需時(shí)間少)比低波速問(wèn)題失穩(wěn)積累得要快;復(fù)雜介質(zhì)問(wèn)題(多個(gè)界面之間反射次數(shù)多)比均勻介質(zhì)問(wèn)題失穩(wěn)積累得快.(3)若邊界出現(xiàn)漂移失穩(wěn)(指邊界節(jié)點(diǎn)運(yùn)動(dòng)向正值或負(fù)值單方向的不正常累積),通常采用修正的邊界公式(3)或(4)來(lái)解決,亦可考慮將高階邊界與低價(jià)邊界組合使用,或者在邊界上附加彈簧和阻尼元件等.(4)若邊界出現(xiàn)振蕩失穩(wěn),可以考慮在整個(gè)計(jì)算模型中增加少量阻尼來(lái)解決,或者在與一些特定的內(nèi)域離散格式結(jié)合時(shí)能夠避免這一問(wèn)題,此外,增大計(jì)算模型,或使用具有耗能特性的時(shí)間積分格式,或采用Liu-Sen 混合邊界方法等都是值得考慮的選擇.

6 結(jié)論

本文新提出離散形式和連續(xù)形式 (即微分方程形式) 兩個(gè)基本邊界公式,通過(guò)證明二者的等價(jià)性并闡明其與經(jīng)典的廖氏透射邊界、旁軸近似邊界、Higdon 邊界以及Givoli-Neta、Hagstrom-Warburton、AWWE 輔助變量邊界等之間的理論聯(lián)系,據(jù)此建立了一種外推型人工邊界條件的理論與公式體系.主要研究結(jié)論如下:

(1) 上述經(jīng)典人工邊界條件可以統(tǒng)一地從外推型ABC 角度來(lái)討論和應(yīng)用.因?yàn)樗鼈冇?jì)算人工邊界節(jié)點(diǎn)運(yùn)動(dòng)的方式都是利用邊界附近一組節(jié)點(diǎn)的運(yùn)動(dòng)(不同邊界使用的節(jié)點(diǎn)組合可能有一定差異)進(jìn)行時(shí)空外推,使用相同的控制參數(shù)(邊界階次和一組計(jì)算波速),并且在數(shù)值模擬中表現(xiàn)出相似的精度和穩(wěn)定性.

(2)本文提出的基本邊界公式在精度和穩(wěn)定性方面均更有優(yōu)勢(shì),是該類邊界最簡(jiǎn)單實(shí)用的形式.

(3) 真正與外推型ABC 不同類型的是以黏性邊界、黏彈性邊界等為代表的應(yīng)力型邊界[40-41]和以函數(shù)衰減層、完美匹配層等為代表的衰減層型邊界[25-26].在后續(xù)研究和應(yīng)用中,外推型ABC[42-43]的穩(wěn)定性、高階應(yīng)力型邊界或完美匹配層邊界應(yīng)用于不同波動(dòng)問(wèn)題的實(shí)現(xiàn)方式與計(jì)算效率、融合不同方法的人工邊界問(wèn)題綜合性解決方案等,都是值得繼續(xù)探討并更好地解決的重要問(wèn)題.

附錄: caj-MTF 邊界的數(shù)值離散格式

圖A1 caj-MTF 邊界數(shù)值計(jì)算格式的局部節(jié)點(diǎn)編號(hào)Fig.A1 A1 Local index numbers of the numerical scheme of caj-MTF boundaries