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

2021-12-21 08:02:04孔曦駿邢浩潔李鴻晶周正華

力學學報 2021年11期

關鍵詞：波場反射系數(shù)波速

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

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

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

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

引言

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

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

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

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

1 MTF 的小量修正形式

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

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

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

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

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

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

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

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

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

如果定義

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

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

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

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

3 反射系數(shù)

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

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

若定義

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

事實上，式(10)給出的就是周-廖方案的邊界反射系數(shù)表達式[2].顯然，若取N=1，本文方法給出的MTF 反射系數(shù)(見式(11)) 退化為式(9)，即一階MTF 反射系數(shù)，這與本文方法在N=1 條件下即為一階MTF 的結論一致.而形如式(10)的周-廖方案的MTF 反射系數(shù)則不同，為了比較，圖1 示出了周-廖方案和本文方案在人工邊界處的反射系數(shù)，從中可以看出消飄因子 γ 對不同階次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ù).在入射角為0 的情形下本文方法的反射系數(shù)保持為0，而周-廖方案無法達到零反射系數(shù).從圖1 中還可知，在本文方法中，消飄因子 γ 取值增大時主要對中等和大角度透射波動的模擬精度造成影響，而對波動能量比較集中的法向和小角度透射范圍影響很小.周-廖方案和本文方法的邊界反射系數(shù)均隨著消飄因子 γ 取值單調(diào)增大，在 γ 取值較小時兩種方案的反射系數(shù)差別很小，但當 γ 取值較大時本文方法的反射系數(shù)明顯要小得多.因此，本文方法能夠適應更大的消飄因子取值范圍(如1.0 以內(nèi)的正數(shù))

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

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

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

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

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

4 數(shù)值算例

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

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

4.1 波源問題

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

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

以初始波場輸入的波源問題作為研究對象，設定以(0.5，0)為圓心，半徑為0.45 的圓形擴散波場作為初始波場.其表達式如下

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

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

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

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

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

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

圖4 不同方法情況下測點A 的位移Fig.4 Displacement at point A by using different methods

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

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

圖5 人工波速對波場精度的影響Fig.5 The effect of artificial wave velocity on the accuracy of the wave field

4.2 散射問題

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

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

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

圖7 為采用不同MTF 邊界方案計算得到的波場快照圖.圖中采用了兩種不同尺度的顏色標尺.對0.2 s，0.4 s 和0.8 s 的波場快照采用尺度較大的顏色標尺，便于觀察波動傳播過程.由于波到達邊界處的反射比較小，對1.6 s 和1.9 s 的波場快照采用尺度較小的顏色標尺，以便比較不同方案的效果.

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

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

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

圖8 測點A 和測點B 的位移時程Fig.8 Displacement time histories at point A and point B

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

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

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

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

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

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

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

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

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

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

如果定義記號

則當m=0 時由式(17)立即得到

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

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

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

6 結論

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

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

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

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

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

附錄A