多次透射公式飄移問題的控制方法1)

孔曦駿 邢浩潔 李鴻晶，2) 周正華

* (南京工業(yè)大學(xué)工程力學(xué)研究所，南京 211816)

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

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

引言

多次透射公式(multi-transmitting formula，MTF)是由廖振鵬等[1-2]提出的一種基于離散參考點(diǎn)表示的局部人工邊界條件[3-4]，具有普適性、易實(shí)施和精度可控等優(yōu)點(diǎn).近年來，MTF 被廣泛地應(yīng)用于復(fù)雜地形場(chǎng)地[5-7]、多種介質(zhì)場(chǎng)地[8-10]和大型結(jié)構(gòu)[11]的地震反應(yīng)分析.還有研究表明[12-15]，MTF與Clayton-Engquist 邊界[16]和Higdon 邊界[17]具有一定的包含關(guān)系，并被應(yīng)用于具有高精度特性的譜元法[18-20]中.高頻振蕩失穩(wěn)[21-23]和飄移失穩(wěn)[24-35]是MTF 可能會(huì)出現(xiàn)的兩種失穩(wěn)現(xiàn)象，本文將研究后者.

對(duì)飄移失穩(wěn)的機(jī)理解釋主要有兩種:一種觀點(diǎn)[24]認(rèn)為將MTF 與有限元法結(jié)合后，用離散模型擬合連續(xù)介質(zhì)產(chǎn)生的波場(chǎng)誤差引發(fā)了失穩(wěn)，且失穩(wěn)導(dǎo)致數(shù)值解將以pN-1的形式增長[24-25]，這里p為計(jì)算步數(shù)，N為MTF 階次;另一種觀點(diǎn)[26-28]認(rèn)為零頻和零波數(shù)的MTF 情形違背了用于判斷雙曲型微分方程初邊值穩(wěn)定性的GKS 準(zhǔn)則，零頻和零波數(shù)成分不斷地從邊界處進(jìn)入波場(chǎng)從而產(chǎn)生了誤差的累積.

基于第1 種解釋，李小軍等發(fā)展了實(shí)時(shí)降階法[24]和數(shù)值輸入法[29]等消飄方法.實(shí)時(shí)降階法認(rèn)為一階MTF不會(huì)發(fā)生失穩(wěn)，高階MTF 則可能發(fā)生失穩(wěn).在應(yīng)用高階MTF 的過程中，通過判斷失穩(wěn)現(xiàn)象的發(fā)生將高階MTF 降為一階MTF，然后再適時(shí)調(diào)回高階.因此，實(shí)時(shí)降階法需要不斷地進(jìn)行失穩(wěn)判斷.事實(shí)上，杜修力等[30]的研究表明一階MTF 也是可能發(fā)生失穩(wěn)的，只是在大多數(shù)情況下不易發(fā)生[31].數(shù)值輸入法則是建立邊界計(jì)算區(qū)，從而獲得入射波離散數(shù)值解，代替連續(xù)介質(zhì)解析解作為波動(dòng)的輸入，減少誤差波場(chǎng)的產(chǎn)生.基于第2 種解釋，周正華和廖振鵬[26]提出的主要消飄方法是對(duì)MTF 添加消飄因子 γ 得到一種修正形式，形式上表現(xiàn)為對(duì)MTF 的各個(gè)時(shí)空外推節(jié)點(diǎn)項(xiàng)進(jìn)行衰減.該方法在物理概念上被理解為考慮了介質(zhì)的幾何擴(kuò)散特性或增加了阻尼特征[27].李小軍和楊宇[32]在MTF 的基礎(chǔ)上對(duì)邊界附近的節(jié)點(diǎn)添加阻尼器或彈簧，結(jié)合應(yīng)力型人工邊界來消除失穩(wěn)，本質(zhì)上與添加消飄因子的方法類似.Zhang 和Yu[33]將一階MTF 與高階MTF 進(jìn)行加權(quán)求和，用來解決電磁波差分模擬中的高階MTF 飄移失穩(wěn)問題.在上述幾種消飄方法中，實(shí)時(shí)降階法、數(shù)值輸入法、附加黏彈性元件法和MTF 加權(quán)求和法都增加了編程的復(fù)雜性，并一定程度上降低了計(jì)算效率.添加消飄因子的方法從形式上并沒有過多地改變MTF 的計(jì)算步驟，非常容易實(shí)現(xiàn).然而添加消飄因子的方法對(duì)消飄因子 γ 的取值有較嚴(yán)苛的限制，若γ取值太小，消飄效果有限，而當(dāng) γ 取值太大時(shí)，則計(jì)算精度受到影響[34].同時(shí)，MTF 階次越高，γ 的取值也要求越大[35].因此，在對(duì)不同模型進(jìn)行分析時(shí)，傳統(tǒng)的消飄因子添加方案需要根據(jù)經(jīng)驗(yàn)不斷地調(diào)試γ值，難于掌握和使用.值得注意的是，消飄因子 γ 的存在必然會(huì)降低MTF 的精度.以一維波動(dòng)問題為例，當(dāng)人工波速(MTF 邊界所使用的計(jì)算波速參數(shù))取值等于視波速時(shí)，除了離散舍入誤差外，一階MTF 本身不產(chǎn)生任何精度損失.然而，添加的消飄因子則先天性的對(duì)MTF 進(jìn)行了衰減，使計(jì)算結(jié)果的幅值低于實(shí)際值.由此導(dǎo)致的精度損失與 γ 取值的大小有關(guān).

本文將采用消飄因子法的技術(shù)路徑，提出一種新的消飄因子修正MTF 方案，用于控制多次透射公式的飄移問題.該方法僅針對(duì)二階及其以上各階MTF 的透射誤差項(xiàng)進(jìn)行修正，而保持一次透射項(xiàng)不變，從而在保有消飄因子法的消飄能力和高階MTF適應(yīng)不同人工波速能力的前提下提供了至少一階MTF 的精度保證.在此基礎(chǔ)上，進(jìn)一步將該方案推廣至任意m階MTF 情形，分別從反射系數(shù)和波動(dòng)有限元數(shù)值算例的角度，對(duì)基于不同 γ 取值的消飄方案的精度和消飄能力進(jìn)行對(duì)比分析，從而驗(yàn)證本文方法的可行性與適用性.

1 MTF 的小量修正形式

MTF 是基于多次透射的概念建立起來的，即假定所有外行波及其透射誤差都是具有相同性質(zhì)的外行單向波動(dòng)，它們都以相同的人工波速沿著法線方向從人工邊界處透射出去.MTF 的一般表達(dá)式寫為

式中，N為透射階次，為二項(xiàng)式系數(shù)，表示參考點(diǎn)j在t=pΔt時(shí)刻的波場(chǎng)值，由下式定義

其中Δt為時(shí)間步距，ca為人工波速.這里人工波速指的是MTF 邊界所使用的計(jì)算波速參數(shù)，理論上它可以設(shè)定為任意值，在實(shí)際計(jì)算中通常取為等于或略大于介質(zhì)物理波速的某個(gè)值.

為了抑制MTF 飄移失穩(wěn)，周正華和廖振鵬[26，27]提出了一種在MTF 中直接添加帶有幾何擴(kuò)散特性或阻尼特征的小量 γ (本文稱之為消飄因子)的措施，即將MTF 表達(dá)式(1)修正為如下形式

這里消飄因子 γ 被規(guī)定為取值遠(yuǎn)小于1 的正數(shù)，因而其對(duì)有限元或有限差分模擬的波動(dòng)分量只產(chǎn)生微小影響，可以忽略不計(jì)[2].本文中簡稱該措施為周-廖方案(Zhou-Liao’s method)，附錄A 給出了該方案在前4 階MTF 情況下的具體表達(dá)式.

根據(jù)GKS 準(zhǔn)則，添加消飄因子 γ 后的MTF 防止了零頻率和零波數(shù)的內(nèi)行波進(jìn)入計(jì)算區(qū)，從而避免了邊界節(jié)點(diǎn)的飄移失穩(wěn).然而，消飄因子 γ 的取值要求比較嚴(yán)格，既不能過小亦不可太大:若 γ 取值過小，則抑制零頻和零波數(shù)內(nèi)行波穿過人工邊界而進(jìn)入計(jì)算域的效果不明顯，邊界節(jié)點(diǎn)處位移依然可能會(huì)發(fā)生飄移失穩(wěn);若 γ 取值較大，由式(3)可知，人工邊界節(jié)點(diǎn)處的位移值將被大大壓低，這會(huì)嚴(yán)重地?fù)p失MTF 模擬精度.故 γ 值的確定需要平衡失穩(wěn)抑制與精度損失這對(duì)矛盾，具體實(shí)施時(shí)需要經(jīng)過繁瑣的試錯(cuò)和調(diào)試尋求比較合理的 γ 取值.

2 MTF 飄移失穩(wěn)控制方案

一階MTF 不易發(fā)生飄移失穩(wěn)現(xiàn)象，但精度有限.當(dāng)人工波速偏離視波速越多，其擬合波動(dòng)的精度就越低.此情形下可通過高階MTF 的方式提高模擬精度，而高階MTF 通常又會(huì)帶來失穩(wěn)問題.通過在MTF 上添加小量的消飄因子能夠有效地抑制MTF 的飄移失穩(wěn)，從而保持人工波速與視波速不一致情況下使用高階MTF 時(shí)的穩(wěn)定性.觀察式(3)可以發(fā)現(xiàn)，這種消飄因子添加方案必然會(huì)降低所有階次MTF 的模擬精度，如果控制不好將導(dǎo)致應(yīng)用MTF 模擬波動(dòng)問題時(shí)失真.

事實(shí)上，高階MTF 是對(duì)誤差波再次實(shí)施透射的結(jié)果，故僅針對(duì)誤差波自身進(jìn)行調(diào)整即可實(shí)現(xiàn)抑制失穩(wěn)的目標(biāo)，從而將MTF 精度損失盡可能控制在較小范圍內(nèi).考慮將一階MTF 表達(dá)為如下形式

如果定義

則N階MTF 可以統(tǒng)一地寫為

為實(shí)現(xiàn)抑制飄移失穩(wěn)的目的，將式(5)改寫為

將式(7)代入式(6)，即得到MTF 飄移失穩(wěn)控制方案.顯然，按照這種控制方案，經(jīng)過修正的前3 階MTF 表示為

如此修正MTF 所形成的消飄方案是基于一階MTF 不易失穩(wěn)這一前提.容易發(fā)現(xiàn)，該方案僅對(duì)透射誤差項(xiàng)進(jìn)行了衰減處理，而一階透射項(xiàng)保持不變，因而至少能夠保留一階MTF 的精度.當(dāng)N=1 時(shí)，該方法實(shí)際上即轉(zhuǎn)變?yōu)橐浑AMTF.

3 反射系數(shù)

下面通過人工邊界反射系數(shù)分析消飄因子 γ 對(duì)修正MTF 精度的影響.不失一般性，這里使用文獻(xiàn)[1]第4.2.3 節(jié)所定義的反射系數(shù).令ca=c(這里c代表物理波速)，當(dāng)輸入波為理想暫態(tài)波時(shí)，一階MTF 的反射系數(shù)寫為

式中，θ為入射平面波與人工邊界法線的夾角，T為平面諧波的周期.

若定義

則按照本文方法，MTF 反射系數(shù)可表達(dá)為

事實(shí)上，式(10)給出的就是周-廖方案的邊界反射系數(shù)表達(dá)式[2].顯然，若取N=1，本文方法給出的MTF 反射系數(shù)(見式(11)) 退化為式(9)，即一階MTF 反射系數(shù)，這與本文方法在N=1 條件下即為一階MTF 的結(jié)論一致.而形如式(10)的周-廖方案的MTF 反射系數(shù)則不同，為了比較，圖1 示出了周-廖方案和本文方案在人工邊界處的反射系數(shù)，從中可以看出消飄因子 γ 對(duì)不同階次MTF 反射系數(shù)的影響.

圖1 邊界反射系數(shù)比較(Δt/T=1/20)Fig.1 Comparison of boundary reflection coefficients with respect to different γ (Δt/T=1/20)

圖1 邊界反射系數(shù)比較(Δt/T=1/20) (續(xù))Fig.1 Comparison of boundary reflection coefficients with respect to different γ (Δt/T=1/20) (continued)

分別比較圖1(a)、圖1(b)和圖1(c)，可以發(fā)現(xiàn)，MTF 階次相同時(shí)，本文方法的邊界反射系數(shù)都明顯低于周-廖方案的反射系數(shù).在入射角為0 的情形下本文方法的反射系數(shù)保持為0，而周-廖方案無法達(dá)到零反射系數(shù).從圖1 中還可知，在本文方法中，消飄因子 γ 取值增大時(shí)主要對(duì)中等和大角度透射波動(dòng)的模擬精度造成影響，而對(duì)波動(dòng)能量比較集中的法向和小角度透射范圍影響很小.周-廖方案和本文方法的邊界反射系數(shù)均隨著消飄因子 γ 取值單調(diào)增大，在 γ 取值較小時(shí)兩種方案的反射系數(shù)差別很小，但當(dāng) γ 取值較大時(shí)本文方法的反射系數(shù)明顯要小得多.因此，本文方法能夠適應(yīng)更大的消飄因子取值范圍(如1.0 以內(nèi)的正數(shù))

由此可見，對(duì)于引入 γ 帶來的精度損失本文方法比周-廖方案要小.比較同一種方案中的不同階次MTF 情況，容易發(fā)現(xiàn)MTF 的階次越高則反射系數(shù)越低，這符合高階MTF 的精度比低階MTF 精度高的特性.

下面從反射系數(shù)公式的角度來闡釋本文方法反射系數(shù)低的原因.一階MTF 方案、周-廖方案和本文方法的邊界反射系數(shù)分別按照式(9)、式(10)和式(11) 進(jìn)行計(jì)算.若分別以R1，R2和R3表示θ=0 情形下上述3 種方案的反射系數(shù)，則有

由于消飄因子 γ 為小值正數(shù)，周-廖方案的反射系數(shù)R2必定大于0，即周-廖方案必定導(dǎo)致一定的精度損失.由式(13)不難發(fā)現(xiàn)，消飄因子 γ 取值越大，周-廖方案的反射系數(shù)亦越大.而一階MTF 方案則不受 γ 的影響，在該情形下的反射系數(shù)R1恒為0.本文方法建立在一階MTF 的基礎(chǔ)上，故該情形下的反射系數(shù)R3同樣恒為0，即不會(huì)因?yàn)橄h因子 γ 的引入而產(chǎn)生精度損失.

當(dāng)θ逐漸增大時(shí)，人工波速ca與邊界上實(shí)際波速的差值越大，MTF 的精度越低，反射系數(shù)也越大(見圖1).此時(shí)，消飄因子 γ 取值越大，周-廖方案和本文方案的反射系數(shù)整體上都有變大的趨勢(shì)，但本文方案的增長趨勢(shì)較小，而周-廖方案的增長趨勢(shì)則十分顯著.這也是引入消飄因子 γ 后對(duì)本文方案的精度影響低于傳統(tǒng)方案的原因.

值得指出，MTF 在實(shí)際波動(dòng)模擬中可以采用不同的插值方案[1，12，21，31，36]來實(shí)現(xiàn).Wang 和Tang[36]研究了MTF 一種簡潔直觀的線性內(nèi)插方案，提出利用Taylor 展開導(dǎo)出與之等效的微分方程形式，進(jìn)而導(dǎo)出理論反射系數(shù)的分析方法.他們發(fā)現(xiàn)此時(shí)人工波速ca的最優(yōu)取值與介質(zhì)物理波速c有所不同，這表明插值方案對(duì)MTF 的反射系數(shù)有重要影響.廖振鵬[1，21]以及邢浩潔等[12，31]對(duì)這個(gè)問題亦做過一定探討.顯然，本文所提出的消飄因子修正MTF 可以結(jié)合不同的插值方案對(duì)其反射系數(shù)進(jìn)行進(jìn)一步理論分析，不過，具體插值方案并不影響本文方法的適用性，限于篇幅，本文對(duì)此不做展開論述.

4 數(shù)值算例

反射系數(shù)分析從解析解的角度客觀地評(píng)估了各參數(shù)對(duì)計(jì)算精度的影響.然而，考慮到實(shí)際波場(chǎng)構(gòu)成的復(fù)雜性以及有限元與人工邊界相結(jié)合的特性，對(duì)相關(guān)數(shù)值模型進(jìn)行時(shí)域分析是必要的.

下面分別以初始波場(chǎng)輸入的二維波源模型和平面波輸入的二維散射模型為例，檢驗(yàn)本文方法的可靠性.為便于分析比較，輸入波的類型皆為SH 波.將前述MTF 人工邊界方案施加到對(duì)應(yīng)場(chǎng)地的有限元離散模型中，并通過集中質(zhì)量的顯式時(shí)域有限元法[1]計(jì)算出各種方案的波場(chǎng)值.在本文數(shù)值算例中，MTF 時(shí)空外推節(jié)點(diǎn)的位移由有限元節(jié)點(diǎn)位移用三點(diǎn)內(nèi)插方法得到[26]，詳見附錄B.所有數(shù)值計(jì)算皆通過自編程序完成.

4.1 波源問題

圖2 為二維均勻介質(zhì)模型，分為半空間模型Ω2和全空間模型 Ω3.半空間模型 Ω2的左側(cè)為前述各種方案對(duì)應(yīng)的人工邊界，其余三側(cè)為自由邊界.全空間模型 Ω3四邊皆為自由邊界.圖中，Ω1為計(jì)算分析區(qū)域，五角星位置為初始波場(chǎng)的中心位置，圓點(diǎn)A(0，0)為測(cè)點(diǎn)位置.設(shè)定介質(zhì)的物理波速為c=1，計(jì)算時(shí)間為2，這將確保在計(jì)算時(shí)間內(nèi)自由邊界的反射波不會(huì)返回到分析區(qū)域 Ω1.

圖2 二維波源問題模型Fig.2 Two-dimensional model of wave source problem

以初始波場(chǎng)輸入的波源問題作為研究對(duì)象，設(shè)定以(0.5，0)為圓心，半徑為0.45 的圓形擴(kuò)散波場(chǎng)作為初始波場(chǎng).其表達(dá)式如下

式中，r2=(x-0.5)2+y2.

由文獻(xiàn)[13]可知，圓形擴(kuò)散波場(chǎng)的最短波長成分約為0.25.為滿足網(wǎng)格尺寸選取要求和時(shí)域積分穩(wěn)定性條件，x和y方向的有限單元尺寸取為0.025，時(shí)步長度取為0.01.為便于進(jìn)行二維波場(chǎng)的誤差分析，通過計(jì)算誤差波場(chǎng)的弗羅貝尼烏斯(Frobenius)范數(shù)，引入最大誤差波場(chǎng)比值Error(U) 的定義，計(jì)算公式如下

式中，U為各個(gè)時(shí)刻波動(dòng)位移的矩陣形式，i和j表示計(jì)算域內(nèi)所有節(jié)點(diǎn)的編號(hào)和位移分量.下標(biāo) Ω1表示計(jì)算節(jié)點(diǎn)的選取區(qū)域.上標(biāo) Ω2和 Ω3分別表示計(jì)算對(duì)應(yīng)的半空間模型和全空間模型.t為計(jì)算對(duì)應(yīng)的時(shí)刻或時(shí)間范圍.誤差波場(chǎng)為半空間模型 Ω2中計(jì)算域Ω1的數(shù)值解減去全空間模型 Ω3中計(jì)算域 Ω1的精確解(即不受邊界影響的大區(qū)域數(shù)值解)，可表示為 ΔΩ1.

經(jīng)過大量的數(shù)值模擬發(fā)現(xiàn)，當(dāng)MTF 為3 階時(shí)，上述模型的邊界節(jié)點(diǎn)發(fā)生了飄移現(xiàn)象.針對(duì)3 階MTF 的飄移失穩(wěn)，分別使用周-廖方法和本文方法對(duì)MTF 人工邊界進(jìn)行修正.圖3 為不同方案的人工邊界情況下區(qū)域 Ω1的位移全波場(chǎng)快照和誤差波場(chǎng)快照.這里的人工波速取為物理波速(ca=c).其中，圖3(a)為未被修正過的3 階MTF 人工邊界，圖3(b)和圖3(c)、圖3(d) 和圖3(e) 分別為消飄因子γ=0.01 和γ=0.1 時(shí)的周-廖方法與本文方法.全波場(chǎng)(Ω1)快照的顏色標(biāo)尺范圍較大，主要反應(yīng)了波的傳播過程.誤差波場(chǎng)(ΔΩ1)快照的顏色標(biāo)尺范圍較小，便于比較不同方案的消飄效果和精度.圖4 為上述不同方法情況下，人工邊界處測(cè)點(diǎn)A的位移時(shí)程圖(測(cè)點(diǎn)A位移用u表示).

從圖3 和圖4 中不難看出，3 階MTF 情況下的波場(chǎng)位移從人工邊界處測(cè)點(diǎn)A附近發(fā)生飄移;不同的MTF 修正方法都一定程度上消除了飄移失穩(wěn)的趨勢(shì).比較圖3(b)～圖3(d)中的全波場(chǎng)快照，可以發(fā)現(xiàn)消飄因子 γ 越大，兩種MTF 修正方法的消飄效果越好.而比較圖3 中的誤差波場(chǎng)和圖4 中t=0.5～2.0 之間的位移時(shí)程，可以發(fā)現(xiàn)消飄因子 γ=0.1 時(shí)的本文方法在所有方案中最接近全空間模型的精確解.此外，當(dāng) γ=0.1 時(shí)，周-廖方法下的測(cè)點(diǎn)A 位移雖然在波峰處接近精確解，但在波谷處比其他方案都遠(yuǎn)離精確解.

圖3 不同方法情況下的全波場(chǎng)快照和誤差波場(chǎng)快照 (ca=c)Fig.3 Snapshots of full-field wave propagation and the error by using different methods (ca=c)

圖4 不同方法情況下測(cè)點(diǎn)A 的位移Fig.4 Displacement at point A by using different methods

為了具體比較消飄因子 γ 對(duì)周-廖方法和本文方法精度的影響，用式(16)中定義的最大誤差波場(chǎng)比值Error(U) 作為參考標(biāo)準(zhǔn).圖5 給出了不同人工波速情況下，消飄因子 γ 對(duì)一階MTF、3 階MTF、周-廖方法與本文方法波場(chǎng)精度的影響.

當(dāng)ca=0.5c時(shí)，雖然3 階MTF 在該情況下已經(jīng)出現(xiàn)了飄移失穩(wěn)的跡象，但飄移趨勢(shì)相對(duì)比較緩慢，在計(jì)算時(shí)間內(nèi)3 階MTF 的最大誤差波場(chǎng)比值小于一階MTF，表明3 階MTF 的精度要高于一階MTF;當(dāng)ca≥c時(shí)，3 階MTF 受到飄移失穩(wěn)的影響，其最大誤差波場(chǎng)比值大于一階MTF，即3 階MTF 的精度反而低于一階MTF.由于周-廖方法與本文方法建立在3 階MTF 的基礎(chǔ)上，當(dāng)消飄因子 γ 極小時(shí)，兩種方案的精度與3 階MTF 接近，隨著消飄因子的變大，周-廖方法和本文方法的精度都得到提升.當(dāng)ca≤ 2c時(shí)，兩種MTF 修正方法的測(cè)點(diǎn)位移值會(huì)隨時(shí)間而逐步收斂為0，解決了3 階MTF 本身飄移失穩(wěn)的問題.當(dāng)消飄因子 γ 增長到一定值時(shí)，兩種MTF 修正方法的最大誤差波場(chǎng)比值隨 γ 的增長而變大.其中，周-廖方法的最大誤差波場(chǎng)比值最終收斂至0.4 左右，而本文方法的最大誤差波場(chǎng)比值則收斂為一階MTF 對(duì)應(yīng)的值.當(dāng)ca=3c時(shí)，3 階MTF 方案、周-廖方法和本文方法都發(fā)生了較大的失穩(wěn)，只有當(dāng)消飄因子 γ 接近1.0 時(shí)，兩種MTF 修正方法才可能收斂.根據(jù)圖5，在不同的消飄因子和人工波速情況下，本文方法的計(jì)算精度整體上要比周-廖方法高，這一點(diǎn)與反射系數(shù)的分析結(jié)果一致.

圖5 人工波速對(duì)波場(chǎng)精度的影響Fig.5 The effect of artificial wave velocity on the accuracy of the wave field

4.2 散射問題

與波源問題不同，散射問題需要考慮波場(chǎng)的輸入問題.同時(shí)，散射問題需要將邊界節(jié)點(diǎn)的全波場(chǎng)反應(yīng)分解為入射波和散射波，然后對(duì)散射波應(yīng)用MTF人工邊界，使其離開不再回到計(jì)算域內(nèi).

圖6 為SH 平面波輸入的散射問題模型.模型為長300 m、寬150 m 的矩形場(chǎng)地，介質(zhì)密度為ρ=2500 kg·m-3，剪切波速為c=500 m/s.4 個(gè)邊界均為MTF 人工邊界，MTF 的階次為N=4.各個(gè)邊界上MTF 的人工波速取為對(duì)應(yīng)邊界的視波速.SH 波從模型的左側(cè)邊界和底部邊界斜入射進(jìn)入波場(chǎng)，入射角為θ=30°.入射波的波形為Ricker 波，中心頻率取為f=10 Hz[15].有限單元網(wǎng)格的尺寸取為2 m，時(shí)間步長取為0.001 s.圖中左側(cè)底角的點(diǎn)A(0，0)和模型中間的點(diǎn)B(150，74)為觀測(cè)點(diǎn).

圖6 二維散射問題模型Fig.6 Two-dimensional model of a scattering problem

圖7 為采用不同MTF 邊界方案計(jì)算得到的波場(chǎng)快照?qǐng)D.圖中采用了兩種不同尺度的顏色標(biāo)尺.對(duì)0.2 s，0.4 s 和0.8 s 的波場(chǎng)快照采用尺度較大的顏色標(biāo)尺，便于觀察波動(dòng)傳播過程.由于波到達(dá)邊界處的反射比較小，對(duì)1.6 s 和1.9 s 的波場(chǎng)快照采用尺度較小的顏色標(biāo)尺，以便比較不同方案的效果.

圖7(a)為未被修正過的4 階MTF，從其1.6 s和1.9 s 的波場(chǎng)快照看，左側(cè)人工邊界和底部人工邊界率先發(fā)生了比較明顯的飄移失穩(wěn)，并且該失穩(wěn)現(xiàn)象隨著時(shí)間的進(jìn)行而向內(nèi)域蔓延.圖7(b)和圖7(c)分別為消飄因子 γ=0.01 時(shí)的周-廖方法和本文方法，從兩圖 1.6 s 和1.9 s 的波場(chǎng)快照看，兩種MTF 修正方法都比較好的消除了4 階MTF 的飄移失穩(wěn)現(xiàn)象.圖7(d)為消飄因子 γ=1.0 時(shí)的周-廖方法，從波場(chǎng)快照看，殘余波場(chǎng)的值較大，其計(jì)算精度不高.這是因?yàn)?γ 超過一定值后，其值越大，周-廖方法精度越低，入射波到達(dá)邊界后產(chǎn)生了較大的反射.圖7(e)為消飄因子 γ=1.0 時(shí)的本文方法，從其0.8 s，1.6 s 和1.9 s的波場(chǎng)快照看，該方案既消除了飄移失穩(wěn)現(xiàn)象，又保持了較高的計(jì)算精度.這是因?yàn)楸疚姆椒ǖ木认孪逓橐浑AMTF，且人工波速為視波速，入射波到達(dá)邊界后產(chǎn)生的反射較小.

圖7 不同方法情況下的波場(chǎng)快照Fig.7 Snapshots of wave propagation by using different methods

圖8 為測(cè)點(diǎn)A和測(cè)點(diǎn)B在不同方法情況下的位移時(shí)程圖及局部放大圖(測(cè)點(diǎn)位移用u表示).圖8(b)和圖8(d)分別為圖8(a)和圖8(c)中綠色矩形框處的局部放大圖.從圖8(a) 中可以明顯看出，4 階MTF 對(duì)應(yīng)的位移時(shí)程發(fā)生了明顯的飄移失穩(wěn)現(xiàn)象.添加消飄因子的幾種方案都較好的解決了飄移失穩(wěn)，但都存在一定的殘余波場(chǎng)，尤其是 γ=1.0 時(shí)的周-廖方法，精度較差，見圖8(c).在這幾種消飄因子添加方案中，γ=1.0 時(shí)的本文方法是殘余波場(chǎng)幅值最小的，精度最高，見圖8(b)和圖8(d).

圖8 測(cè)點(diǎn)A 和測(cè)點(diǎn)B 的位移時(shí)程Fig.8 Displacement time histories at point A and point B

與前面波源問題相似，這里以測(cè)點(diǎn)A和測(cè)點(diǎn)B在t=0.5 s～2.0 s 時(shí)間范圍內(nèi)的最大位移幅值|U|max為參考標(biāo)準(zhǔn)，來更具體的比較消飄因子 γ 對(duì)周-廖方法與本文方法計(jì)算精度的影響，如圖9 所示.從圖9可以看出，當(dāng)消飄因子 γ ＜ 0.1 時(shí)，周-廖方法與本文方法的計(jì)算精度十分接近且?guī)缀醪浑S γ 取值而改變;當(dāng) γ ＞ 0.1 時(shí)，本文方法的誤差隨著 γ 的增大而逐漸降低為0，而周-廖方法的誤差整體上則隨著 γ 的增大而迅速變大.因此，如果對(duì)計(jì)算精度有最低要求，本文方法的消飄因子取值范圍將明顯大于周-廖方法.這種情況下，有別于周-廖方法在應(yīng)用過程中對(duì)消飄因子 γ 取值需要不斷地進(jìn)行試錯(cuò)，本文方法對(duì)消飄因子的選值更加方便與寬松.

圖9 消飄因子 γ 對(duì)不同方法模擬波場(chǎng)精度的影響Fig.9 Influence of the drift-elimination factor γ on the accuracy of simulated wave field by using different boundary methods

結(jié)合數(shù)值算例的結(jié)果與反射系數(shù)的分析，可以發(fā)現(xiàn)，兩種方式都能反應(yīng)出本文方法相對(duì)于傳統(tǒng)方法的精度優(yōu)勢(shì)，但精度的表現(xiàn)存在一定差異，這主要是因?yàn)?(1) 理論反射系數(shù)基于穩(wěn)態(tài)諧波導(dǎo)出，而本文數(shù)值試驗(yàn)則是短持時(shí)的暫態(tài)波;(2) 反射系數(shù)是對(duì)每個(gè)特定頻率ω(ω=2πf=2π/T)給出的，而數(shù)值試驗(yàn)中的波動(dòng)則具有多種頻率成分;(3) 反射系數(shù)使用的是理想平面波，具有特定的透射角度，而數(shù)值試驗(yàn)中的波動(dòng)以多種不同的角度射向邊界，且往往含有幾何衰減、多波疊加等效應(yīng);(4) 離散網(wǎng)格具有數(shù)值頻散效應(yīng)，會(huì)導(dǎo)致其中傳播的波形發(fā)生變化，對(duì)人工邊界的反射特征造成一定影響.

依據(jù)應(yīng)用本文方法的大量數(shù)值模擬經(jīng)驗(yàn)，對(duì)MTF的階次N和消飄因子 γ 的取值有以下歸納:(1) MTF常用的階次一般不超過3 階，本文算例中使用高階MTF (N=4)是為了研究飄移失穩(wěn)問題，探討消飄方法的效果;此外，MTF 并非只在高階情況下才發(fā)生飄移失穩(wěn)，在有些情形下[2，6，26]，低階MTF 也可能出現(xiàn)飄移失穩(wěn);(2) 關(guān)于消飄因子 γ 的取值范圍，周-廖方法一般為0.001～0.05[1，2，6，26，35]，本文方法在同樣精度情況下則為0.01～1.0，取值范圍更大;對(duì)于飄移比較嚴(yán)重的問題，本文方法對(duì) γ 的取值可以大一些，如0.5 或接近1，對(duì)于飄移比較輕微的問題，γ 的取值則可以小一些，如0.01 左右.

5 基于任意階MTF 的消飄通用形式

前面討論的MTF 飄移失穩(wěn)控制方案是立足于一階MTF 上而建立起來的，即保持一階MTF 不變，而僅對(duì)二階及其以上各階MTF 透射誤差進(jìn)行修正，這樣就在實(shí)現(xiàn)抑制飄移失穩(wěn)的同時(shí)至少保留了一階MTF 的精度.這種MTF 飄移失穩(wěn)控制方案的思路很容易推廣至任意階MTF 情形.如果認(rèn)為前m階MTF 都能夠保持穩(wěn)定而不發(fā)生失穩(wěn)，則只需針對(duì)m+1 階及其以上的各階MTF 透射誤差添加消飄因子，即可達(dá)到抑制飄移失穩(wěn)的目的.考慮到m階MTF 可以表達(dá)為

則只需將N階MTF (N＞m)修正為

此即為任意階MTF 飄移失穩(wěn)控制的通用形式.顯然，若取m=1，則式(18)與本文前面導(dǎo)出的MTF飄移失穩(wěn)控制方案完全一致.

式(18)對(duì)應(yīng)的邊界反射系數(shù)可以統(tǒng)一地寫為

不難發(fā)現(xiàn)，當(dāng)取m=0 時(shí)，由式(19)給出的反射系數(shù)實(shí)際上就是周-廖方案的反射系數(shù)表達(dá)式.下面證明周-廖方案是本文方法的一個(gè)特例.

如果定義記號(hào)

則當(dāng)m=0 時(shí)由式(17)立即得到

即邊界點(diǎn)始終保持為0，這表明零階MTF 實(shí)際上是固定邊界條件.在該情形下，根據(jù)式(17)并利用式(18)容易導(dǎo)出各階MTF 如下

將上述各階MTF 寫成通式，有

比較發(fā)現(xiàn)，式(22)與式(3)完全相同，即其實(shí)際上就是周-廖方案所給出的N階MTF.故無論從MTF 表達(dá)式(式(22))還是邊界反射系數(shù)(式(19))上看，周-廖方案本質(zhì)上只是本文方法取m=0 時(shí)的一個(gè)特例.

6 結(jié)論

多次透射公式(MTF)建立起來的人工邊界條件有時(shí)存在飄移失穩(wěn)問題.本文將周正華和廖振鵬提出的抑制MTF 飄移失穩(wěn)的方案推廣到一般情形，發(fā)展了一種多次透射公式飄移問題的控制方法.通過邊界反射系數(shù)分析和數(shù)值試驗(yàn)，對(duì)本文方法的精度和消飄效果進(jìn)行了詳細(xì)探討，得到主要研究結(jié)論如下.

(1) 本文方法能夠十分有效地消除高階MTF 與有限元模型結(jié)合后出現(xiàn)的飄移失穩(wěn)現(xiàn)象，僅通過添加消飄因子的方式實(shí)現(xiàn)控制MTF 飄移失穩(wěn)的目標(biāo)，這對(duì)于MTF 實(shí)施步驟及其編程基本不會(huì)產(chǎn)生額外負(fù)擔(dān).

(2) 本文的MTF 飄移控制方案是在m階MTF的基礎(chǔ)上構(gòu)建出來的，僅針對(duì)(m+1)階及其以上各階MTF 的透射誤差項(xiàng)進(jìn)行修正，能夠保證波動(dòng)模擬結(jié)果至少保持m階MTF 精度(當(dāng)取m=1 時(shí)即至少保有一階MTF 精度).傳統(tǒng)的周-廖方案本質(zhì)上是本文方法的一個(gè)特例，相當(dāng)于取m=0 時(shí)的結(jié)果.

(3) 在本文方法中，消飄因子 γ 取值增大時(shí)主要對(duì)中等和大角度透射波動(dòng)的模擬精度造成影響，而對(duì)波動(dòng)能量比較集中的法向和小角度透射范圍影響很小.它能夠適應(yīng)更大的消飄因子取值范圍(如1.0 以內(nèi)的正數(shù))，且較傳統(tǒng)方案能夠在更高精度水平下實(shí)現(xiàn)對(duì)MTF 飄移問題的有效控制.

(4) 當(dāng)消飄因子取值趨向于無限大時(shí)，由傳統(tǒng)方案推算的邊界運(yùn)動(dòng)趨近于零，此情形相當(dāng)于固定邊界條件，而本文方法則趨近于m階MTF(實(shí)際應(yīng)用時(shí)可只取為一階MTF).從原理上看，本文方法融合了零頻零波數(shù)情形下GKS 準(zhǔn)則所表述的消飄原理[26，27]和通常認(rèn)為低階MTF 不易發(fā)生飄移失穩(wěn)的經(jīng)驗(yàn)觀察[24，33].

附錄A

附錄B

本文MTF 外推節(jié)點(diǎn)插值格式

MTF 外推節(jié)點(diǎn)位移通過有限元節(jié)點(diǎn)位移的3 點(diǎn)內(nèi)插方式獲得.

如圖B1 所示，多次透射公式可以看成是基于時(shí)空節(jié)點(diǎn)信息的有限差分外推公式.其形式如下

圖B1 透射邊界(MTF)示意圖Fig.B1 Schematic diagram of multi-transmitting formula (MTF)

以第一個(gè)時(shí)空外推節(jié)點(diǎn)為例(圖中節(jié)點(diǎn)1a，一次透射項(xiàng))，令s=caΔt/Δx，Δx為相鄰的有限元網(wǎng)格節(jié)點(diǎn)之間的空間間距，基于與同一時(shí)刻相鄰3 個(gè)有限元網(wǎng)格節(jié)點(diǎn)位移的空間插值關(guān)系，有

同理，令sj=jcaΔt/Δx，可求第j個(gè)外推時(shí)空節(jié)點(diǎn)的位移，以矩陣形式表示如下

將式(B3)代入式(B1)，可得用有限元網(wǎng)格離散節(jié)點(diǎn)的位移表示的MTF 形式，即