插值算法對透射邊界數值模擬結果的影響分析

[摘要]" " 一個時間步內外行波以人工波速走過的距離與網格尺寸往往不一致，因此在應用透射邊界公式時，一般通過時間插值或空間插值的方法將透射邊界的計算結點用有限元結點表示。本文導出了拉格朗日插值算法、埃爾米特插值算法以及三次樣條插值算法的透射邊界一次透射插值公式，探討不同插值算法對透射邊界數值模擬結果的影響。首先將彈性半空間中入射P/SV波自由場問題的數值解與解析解對比，討論3種插值算法對計算精度的影響，然后以半圓河谷地形入射P/SV波散射問題為例，對比不同插值算法計算得到的結果。數值結果表明：不同插值算法對透射邊界數值模擬結果有顯著影響，在3種插值算法中，三次樣條插值算法計算精度最高，其次是拉格朗插值算法，相比之下，埃爾米特插值算法計算得到的精度較低。相比于自由場，散射問題的計算結果受3種插值算法的影響較大，拉格朗日插值算法與三次樣條插值算法計算結果相近，而埃爾米特插值算法的計算結果差異較大。本研究成果為在應用透射邊界模擬場地動力反應時，選擇合適的插值算法提供了科學依據。

[關鍵詞] 透射邊界; 插值算法; 自由場; 散射場

[DOI] 10.19987/j.dzkxjz.2024-046

基金項目： 國家自然科學基金項目（U2039208）資助。

0" 引言

近場波動數值模擬的關鍵性問題在于建立近場區域的邊界條件，由于這一邊界是人為設置的，故稱之為人工邊界條件[1]。人工邊界要以保證有限近場區域與無限外部區域之間波能量的交換為設置原則，其又根據是否能夠精確滿足外部無限域內的所有場方程和物理邊界條件分為全局人工邊界和局部人工邊界[2]。

全局人工邊界主要包括邊界元法[3-4]、波函數展開法[5]、薄層法[6]等。由于全局人工邊界不僅將所有邊界節點的運動耦聯起來，而且還將全部過去時刻的運動與當前時刻的運動耦聯起來，這就導致在應用全局人工邊界進行時空域內的數值模擬時存在計算量巨大的問題，難以滿足實際需要[2]。相比之下，局部人工邊界的主要特征是時空解耦，即一個邊界節點在某一時刻的運動僅與其臨近節點臨近時刻的運動有關[1]，這極大地簡化了數值計算，因此逐漸成為研究人員關注的焦點。

局部人工邊界是模擬外行波穿過人工邊界向無窮遠傳播的性質來實現的[2]，主要包括黏彈性邊界[7-8]、Higdon邊界[9]、完美匹配層邊界[10]和透射邊界[11]等。黏彈性邊界是針對特殊波場建立的，一般情況下其僅有零階精度，但是實現簡單，對于工程所關心的、遠離人工邊界的結構反應而言具有實用價值[12-15]。Higdon邊界是能夠完全吸收不同入射角入射的若干平面波的高階人工邊界，但存在高次時空導數，難以直接采用位移有限元對其進行離散，同時構建差分格式也較為困難。完美匹配層邊界是可以實現外行波無反射地進入邊界吸收層的人工邊界方案，但需要引入復變量，因此在數學上，實現起來較為復雜。相較于上述人工邊界，透射邊界具有以下優勢：①透射邊界是從一般無限域模型出發，可直接用于各種場方程，適用性廣；②透射邊界是由單向波導出，已計入無限域中邊界的影響，可適用人工邊界上的各節點，包括角點和分層界面點；③透射邊界是一維表達式，易于數值實現，且計算量小；④透射邊界的精度階數含義明確且可控。基于以上優勢，很多學者將透射邊界應用到局部地形場地效應研究上。

盡管透射邊界具有很多優勢，但其自身存在高頻振蕩[16]和低頻飄移[17]等穩定性問題，這些問題制約了透射邊界的工程應用。導致透射邊界失穩的主要原因可歸結為：源自數值計算過程中產生的附加內行波能量從人工邊界持續滲透至計算區域內，進而導致有限計算區內波動能量的累計增大。一種是由人工邊界反射引起的附加內行波，當反射系數大于1時，將導致透射邊界的高頻振蕩失穩[16]；另一種是數值積分過程中形成的附加內行波，由于舍入誤差等，將導致能量不斷輸入有限計算區域內，從而導致低頻飄移失穩[17]。針對高頻振蕩失穩，廖振鵬和劉晶波[18]給出了平滑因子抑制高頻失穩的措施；李小軍和劉愛文[19]、周正華和周扣華[20]從顯示積分格式的能耗特性考慮，通過增加與應變速率成正比的阻尼來抑制高頻失穩。對于低頻飄移失穩，李小軍[21]建議采用降階的方式來抑制飄移失穩；周正華等[22]提出了在多次透射公式中引入修正因子的方式來消除飄移失穩的措施。唐暉和李小軍[23]系統地分析了求解散射問題中消除多次透射邊界漂移失穩措施的效果。章旭斌等[24]等探究了透射邊界高頻失穩的機理，并提供了一種穩定實現的方法。另外，透射邊界與普單元法結合，邢浩潔等[25-27]提出了一種高效的人工邊界條件實現方法。

一個時間步內外行波以人工波速走過的距離與網格尺寸往往不一致，因此在應用透射邊界公式時，一般通過時間插值或空間插值的方法將透射邊界的計算結點用有限元結點表示。插值技術作為連接透射邊界理論與有限元分析的重要環節，在以往研究中普遍采用了拉格朗日插值方法對透射邊界節點進行數值插值。然而，關于插值算法選取對透射邊界數值模擬精度及其最終結果的影響探究尚顯不足。基于此，本文采用不同插值算法來計算透射節點，通過對比不同插值算法計算的結果，來分析插值算法對透射邊界數值模擬結果的影響，研究成果為地震工程研究者在應用透射邊界理論進行場地動力反應數值模擬過程時，選擇合適的插值算法提供科學依據。

1" 理論推導

1.1" 透射邊界理論

透射邊界是通過模擬外行波透過邊界的過程來實現邊界截斷的。這一理論的要點包括兩方面：第一，給出穿過人工邊界上一點并沿邊界外法線方向傳播的單向波動的一般表達式，即將單向波動表示成一系列外行平面波的疊加。第二，利用多次透射技術模擬上述一般表達式。首先假定所有單向波動以同一人工波速沿法線方向從邊界透射出去，由此得到一次透射公式，然后證明一次透射的誤差也是具有相同外行性質的單向波動，由此建立了二次進而多次透射公式（MTF）。MTF可以寫作：

u_0^（p+1）=∑_（j=1）^n?〖〖（-1）〗^（j+1） C_j^N u_j^（p-j+1） 〗 （1）

式中，C_j^N為二次多項式系數，離散位移u_j^p表達式如下：

u_j^p=u（jc_a Δt，pΔt） （2）

式中，c_a為人工波速。考慮一次透射邊界，式（1）可以寫為：

u[0\"，\"（p+1）Δt]=u（c_a Δt\"，\" pΔt） （3）

由于在一個時間步內外行波以人工波速走過的距離c_a Δt與網格尺寸Δx往往不一致，因此在應用式（3）時，一般通過時間插值或空間插值的方法將透射邊界的計算結點用有限元結點表示（圖1）。用有限元三節點插值來表示式（3）中等式右端，計算表達式如下：

（u（c_a Δt，pΔt）=amp;t_11 u（0，pΔt）+t_12 u（Δx，pΔt）+@amp;t_13 u（2Δx，pΔt）） （4）

式中，t11，t12和t13是插值系數。本文將推導不同插值算法的一階透射公式，著重探討不同插值算法對數值計算精度的影響。

1.2" 插值系數推導

1.2.1" 插值基本概念

設函數y=f（x）有n個已知數據點[xi，yi]，i = 0，1，…，n，若存在一個簡單函數P（x），使得P（x_i）= yi，i = 0，1，…，n成立。則P（x）為f（x）的插值函數。

1.2.2" 拉格朗日插值

Lagrange插值公式為：

l_i （x）=∏_（j=0@j≠i）^n" （x-x_j）/（x_i-x_j ）（i=0，1，…，n） （5）

據此可以得到插值多項式：

L_n （x）=∑_（i=0）^n y_i l_i （x）=∑_（i=0）^n y_i （∏_（j=0@j≠i）^n" （x-x_j）/（x_i-x_j ）） （6）

然后，將式（6）帶入式（4）中，可以得到由拉格朗日插值法計算得到的一次透射公式的插值系數：

（amp;t_11=1/2（-S+1）（-S+2）@amp;t_12=S（-S+2）@amp;t_13=1/2 S（S-1）） （7）

式中，S=（c_a Δt）/Δx。

1.2.3" 牛頓插值

牛頓插值的一階、二階差商公式分別為：

f[x_i，x_j]=（f（x_i）-f（x_j））/（x_i-x_j ） （8）

f[x_i，x_j，x_k]=（f[x_i，x_j]-f[x_j，x_k]）/（x_i-x_k ） （9）

根據一階、二階差商公式，可以得到N階差商的表達式：

（f[x_0，x_1，…，x_n]=amp;（f[x_0，x_1，…，x_（n-1）]）/（x_0-x_n ）-@amp;（f[x_1，…，x_n]）/（x_0-x_n ）） （10）

當x_i-x_j=Δx為定值時，

f[x_0，x_1，…，x_n]=（-1）/（（n-1）！〖（Δx）〗^（n-1） ） C_（n-1）^i f（iΔx） （11）

由此可得插值多項式：

（amp;N_n （x）=f（x_0）+（x-x_0）f[x_0，x_1]+…+@amp;（x-x_0）（x-x_1）…（x-x_（n-1））f[x_0，x_1，…，x_n]） （12）

將式（12）帶入式（4），可以得到由牛頓插值法計算得到的一次透射公式的插值系數，如下：

（amp;t_11=1/2（-S+1）（-S+2）@amp;t_12=S（2-S）@amp;t_13=S/2（S-1）） （13）

觀察式（7）和式（13），可以發現，拉格朗日與牛頓插值算法一階透射公式的插值系數相同，這是由于兩種插值算法雖表達形式各異，但在本質上具有相同的階次和系數多項式，因此對于高階透射公式，同樣會得出一致的插值系數。

1.2.4" 三次埃爾米特插值

假定已知函數f（x）在插值區間[p，q]上的n+1個互不相同節點xi（i = 0，1，…，n）處滿足f（x_i）=f_i以及f^' （x_i）=〖f_i〗^'（i = 0，1，…，n），如果函數G（x）的存在滿足下列條件：

（1） G（x）在每個小區間上的多項式次數為3；

（2） G（x）∈ C[a，b]；

（3）G（x_i）=f（x_i）以及G^' （x_i）=f^' （x_i）（i = 0，1，…，n）。

則稱G（x）是f （x）在n+1個節點xi上的分段三次埃爾米特插值多項式。

根據埃爾米特插值要求，

（G（x）=amp;（1-3t^2+2t^3）y_0+t（1-2t+t^2）m_0+@amp;t^2 （3-2t）y_1+t^2 （t-1）m_1 ） （14）

其中，t=（x-x_0）/（x_1-x_0），m_0=（u（Δx，t）-u（0，t））/Δx，m_1=（u（2Δx，t）-u（Δx，t））/Δx，進而計算得到插值系數公式：

（amp;t_11=9/5 S^3-14/5 S^2-1/5 S+1@amp;t_12=-2S^3+14/5 S^2+1/5 S@amp;t_13=1/5 S^2 （S-1）） （15）

1.2.5" 三次樣條插值

將區間[a，b]劃分成n個區間，[（x0，x1），（x1，x2），…，（xn-1，xn）]，共有n+1個點，其中兩個端點x0=a，xn=b。三次樣條確保每個小區間的曲線都是一個三次方程，三次樣條方程滿足以下條件：

（1） 在每個分段區間[xi，xi+1]上，滿足S（x） = Si（x）都是一個三次方程；

（2） 滿足插值條件，Si（x） = yi；

（3） 曲線光滑，即S（x），S'（x），S''（x）連續。

采用三點插值，可以將整個區間劃分為兩段，為[0，Δx]與[Δx，2Δx]。在每個區間內求解一個三次方程：

y=a_i+b_i x+c_i x^2+d_i x^3 （i=1，2） （16）

根據式（16）可知，待求未知系數共有8個。為了求解兩個區間內的三次函數，需要得到8個線性無關的方程。滿足插值條件可以提供4個方程，滿足S（x），S'（x）連續條件可以提供2個方程，另外考慮自然邊界條件，即

S^″ （\"0\"）=0=S^″ （\"2\" Δx） （17）

上式也提供了2個方程。至此可以求解每個區間內的三次方程，具體方程如下：

（（1amp;0amp;0amp;0amp;0amp;0amp;0amp;0@1amp;Δxamp;〖（Δx）〗^2amp;〖（Δx）〗^3amp;0amp;0amp;0amp;0@0amp;0amp;0amp;0amp;1amp;Δxamp;〖（Δx）〗^2amp;〖（Δx）〗^3@0amp;0amp;0amp;0amp;1amp;2Δxamp;4〖（Δx）〗^2amp;8〖（Δx）〗^3@0amp;1amp;2Δxamp;3〖（Δx）〗^2amp;0amp;-1amp;-2Δxamp;-3〖（Δx）〗^2@0amp;0amp;2amp;6Δxamp;0amp;0amp;-2amp;-6Δx@0amp;0amp;2amp;0amp;0amp;0amp;0amp;0@0amp;0amp;0amp;0amp;0amp;0amp;2amp;12Δx））（a_1@b_1@c_1@d_1@a_2@b_2@c_2@d_2 ）}={（u_0@u_1@u_1@u_2@0@0@0@0）} （18）

式中，u0 = u（0，pΔt），u1 = u（Δx，pΔt），u2 = u（2Δx，pΔt）。求解式（18）便可得到每個區間內的三次方程。

考慮到caΔt ≤ Δx，因此透射邊界的計算結點落在第一個區間內，滿足第一個區間內的三次方程，可得：

u（c_a Δt，t）=a_1+b_1 c_a Δt+c_1 〖（c_a Δt）〗^2+d_1 〖（c_a Δt）〗^3 （19）

聯合式（19）與式（4），便可得到插值系數：

（amp;t_11=1-5/4 S+1/4 S^3@amp;t_12=3/2 S-1/2 S^3@amp;t_13=-1/4 S+1/4 S^3 ） （20）

2" 數值計算

2.1" 入射波

通過分析自由場來驗證不同插值算法的計算精度，并以半圓河谷地形為例，分析不同插值算法對散射場計算結果的影響。入射波位移脈沖為：

{（（u（t）=amp;16A[Z（t/T）-4Z（t/T-1/4）+6Z（t/T-1/2）-@amp;4Z（t/T-3/4）+Z（t/T-1）]）@Z（a）=a^3 H（a）） （21）

式中，H（a）為赫維賽德函數，t是時間，A是脈沖的峰值，T是脈沖的作用時間。入射波位移脈沖時程與傅里葉頻譜如圖2所示。需要說明的是：對于自由場問題，本文取T=1；對于散射場問題，本文取T=2。

2.2" 有限元模型

本文數值模擬所取的材料參數分別為：泊松比v=0.25，介質密度ρ=2000 kg/m3，彈性模型E= 5×109 Pa。根據式（22）可以計算得出，P/SV波的波速分別為vP=1732 m/s，vS=1000 m/s。

v_S=√（E/（2ρ（1+v））），v_P=√（（E（1-v））/（ρ（1+v）（1-2v））） （22）

假定一個波長包含10個網格，由此可以計算網格的最大尺寸。首先確定入射波的截止頻率，然后根據下式可得：

Δx?λ/10=c_min/10f=1000/（18×10）=5.56\"m\" （23）

要求一個時間步內，波傳播的距離不大于一個網格尺寸，進而可以確定時間步長，計算公式如下：

Δt?Δx/c_max =5/1732≈0.00289\"s\" （24）

本文為了方便計算以及排除網格尺寸與時間步長對計算結果的影響，統一取Δx = 5 m，Δt = 0.001 s。根據選定的時間步長和網格尺寸并結合式（7）、（15）和（20），計算得到了3種插值算法的插值系數（表1）。

2.3" 彈性半空間

本文研究所使用的自由場模型如圖3所示，模型為矩形，長寬分別為400 m和200 m。矩形的左上頂點坐標為（0，0），右下頂點為（400，?200）。入射波從左下方入射到計算區域內，入射點為（0，?200）。A（200，0），B（200，?100）和C（200，?200）為觀測點。

圖4展示了P波垂直入射時，不同插值算法計算得到的場區觀測點的位移時程以及其與解析解對比后的誤差。P波垂直入射時，3種插值算法計算得到的場區觀測點的水平分量位移與解析解完美匹配，誤差幾乎為零。相比于水平分量，豎向分量位移在不同插值算法下計算得到的結果與精確解析解相比，出現了明顯誤差，但絕對誤差最大值均小于0.005。另外，3種插值算法計算得到的誤差曲線走勢基本相同，但不同觀測點最大誤差出現的時間點不盡相同，基本滿足隨著觀測點y坐標值越小誤差最大值出現的時間越靠后的規律。

圖5展示了P波60°斜入射時，不同插值算法計算得到的場區觀測點的位移時程以及其與解析解對比后的誤差。相比于垂直入射，斜入射誤差值明顯增大。拉格朗日插值與三次樣條插值算法計算得到的位移時程誤差曲線基本相近，而埃爾米特插值算法計算得到的位移時程誤差曲線與其他兩種插值算法的結果明顯不同，且埃爾米特插值算法計算得到的誤差更大。3種插值算法計算得到的水平分量位移時程誤差比豎向分量大。另外，隨觀測點y坐標值減小，誤差最大值越大。

圖6展示了SV波垂直入射時，不同插值算法計算得到的場區觀測點的位移時程以及其與解析解對比后的誤差。與P波垂直入射結果相反，3種插值算法計算得到的場區觀測點的豎向分量位移與解析解完美匹配，誤差幾乎為零，而水平分量位移誤差較P波入射時情況顯著增大，絕對誤差最大值達到0.02。3種插值算法對水平分量位移時程影響顯著，對于觀測點C，三次樣條插值算法計算得到水平分量位移時程誤差結果與另外兩種插值算法計算得到的結果明顯不同，其誤差最大值最小，基本小于0.001，而另外兩種插值算法計算的誤差均大于0.01，但上述規律對于觀測點A和B不滿足。整體上，對比于拉格朗日插值和三次樣條插值算法，埃爾米特插值算法計算得到的誤差曲線不同，且誤差更大。

圖7展示了SV波30°斜入射時的計算結果，其規律與P波60°斜入射時相似，最大的不同在于SV斜入射誤差最大值比垂直入射小，且絕對誤差最大值均小于0.01。

2.4" 半圓河谷地形

本文研究所使用的散射場模型如圖8所示，場地的長寬尺寸分別為400 m和200 m，區域的左上頂點坐標為（0，0），右下頂點為（400，?200）。入射波從左下方入射到計算區域內，入射點為（0，?200），半圓形河谷的半徑為50 m，河谷圓心坐標為（200，0）。根據網格劃分的尺寸，沿地表水平方向每隔5 m取一個觀測點，并對觀測點進行編號Pi（i =1，2，…，81）。

垂直入射時，觀測點位移時程關于地表中點中心對稱，因此選取左半區域中的4個典型觀測點來分析其位移時程受插值算法的影響。4個典型觀測點分別為P1，P16，P31和P41，其位置如圖8所示。圖9展示了4個觀測點在不同插值算法計算下P、SV波垂直入射時位移時程的對比結果。P波垂直入射時，觀測點P41水平分量位移時程始終為0，而對于SV波，觀測點P41豎向分量位移時程為0，且不同插值算法計算結果基本一致。不同插值算法對位移時程的計算結果有顯著影響，特別是靠近河谷區域的觀測點。圖9所示表明拉格朗日插值和三次樣條插值算法計算得到的位移時程結果基本相近，而埃爾米特插值算法計算結果卻顯著不同。

為了進一步分析斜入射時不同插值算法對位移時程計算結果的影響，圖10a展示了P波60°斜入射時的位移時程對比結果。考慮SV入射時存在臨界角（在本文選取的參數下為35.26°）問題，因此計算了SV波30°斜入射時的位移時程，如圖10b所示。本節主要研究在散射問題中插值算法對透射邊界計算結果的影響，不同于垂直入射，地表點的位移時程具有對稱性，斜入射時，局部地形兩側的觀測點的位移時程顯著不同，特別是靠近局部地形的觀測點，需要重點分析，因此本節所考察觀測點為P1，P31，P41和P51。觀察圖10，不難發現相較于靠近河谷的觀測點P31，P41和P51，3種插值算法計算得到的觀測點P1的水平和豎向位移時程結果基本一致。采用拉格朗日插值和三次樣條插值算法計算得到的位移時程結果基本相近，而采用埃爾米特插值算法計算得到的位移時程表現出來顯著的差異，特別是觀測點P31和P51。盡管3種插值算法計算得到的位移時程曲線具有差異，但這種差異多出現在位移時程的后半段，這也導致3種插值算法計算得到的位移時程的峰值點基本相同。

3" 結論

為了將拉格朗日插值、埃爾米特插值和三次樣條插值算法應用到透射邊界中，本文首先導出了不同插值算法的一階透射邊界插值系數，進而溝通了透射節點和有限元節點的關系，為透射理論與有限元理論結合提供了基礎。隨后，通過比較3種插值算法計算得出的自由場觀測點位移時程與理論解析解之間的差異，系統地評估了各插值方法在計算精度方面的表現。此外，還深入探討了這3種插值算法在散射場觀測點位移時程計算中的差異性，以此來詳盡剖析它們在透射邊界數值模擬結果上的影響程度。得出了如下結論：

（1） 插值算法對透射邊界數值模擬結果有影響，特別是對于散射問題。

（2） 3種插值算法中，三次樣條插值算法計算精度最高，其次是拉格朗日插值算法，最后是埃爾米特插值算法。

（3） 散射問題中，拉格朗日插值算法與三次樣條插值算法計算結果相近，而埃爾米特插值算法的計算結果差異較大。

本文采用的一階透射公式數值模擬精度很高，誤差均不大于1%，特別是三次樣條插值算法，精度最高。本文提供了三次樣條插值算法的插值系數計算公式（式（20）），研究人員可根據數值模擬所取網格尺寸和時間步長進行插值系數計算。值得注意的是：在本研究中，所獲得的所有插值系數均基于一階透射公式構建，而對于更高階數的透射公式對應的插值系數，其精確推導有待深入研究。此外，通過大量數值模擬案例的系統性分析，揭示了插值算法對透射邊界的計算結果所產生的影響規律，但欠缺針對這一現象背后數學機制的理論解析。

致謝

本文由國家自然科學基金項目（U2039208）資助。感謝審稿人的建設性意見，這些意見極大地改進了本文。

參考文獻

廖振鵬. 工程波動理論導論[M]. 2版. 北京：科學出版社，2002" " Liao Z P. Introduction to wave motion theories in engineering[M]. 2nd ed. Beijing：Science Press，2022

廖振鵬. 近場波動的數值模擬[J]. 力學進展，1997，27（2）：193-216" " Liao Z P. Numerical simulation of near-field wave motion[J]. Advances in Mechanics，1997，27（2）：193-216

Gao Z Y，Li Z L，Liu Y J. A time-domain boundary element method using a kernel-function library for 3D acoustic problems[J]. Engineering Analysis with Boundary Elements，2024，161：103-112

Liu Z X，Ai T C，Huang L，et al. The scattering of seismic waves from saturated river valley with water layer：Modelled by indirect boundary element method[J]. Engineering Analysis with Boundary Elements，2023，149：282-297

Lee V W，Liu W Y. Two-dimensional scattering and diffraction of P-and SV-waves around a semi-circular canyon in an elastic half-space：An analytic solution via a stress-free wave function[J]. Soil Dynamics and Earthquake Engineering，2014，63：110-119

李小信，何超，周順華，等. 具有不規則界面的層狀地基三維動力響應的薄層法[J]. 巖土力學，2023，44（增刊1）：655-668" " Li X X，He C，Zhou S H，et al. Thin layer method for three-dimensional dynamic response of layered foundation with irregular interfaces[J]. Rock and Soil Mechanics，2023，44（S1）：655-668

劉晶波，王振宇，杜修力，等. 波動問題中的三維時域粘彈性人工邊界[J]. 工程力學，2005，22（6）：46-51" " Liu J B，Wang Z Y，Du X L，et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J]. Engineering Mechanics，2005，22（6）：46-51

杜修力，趙密. 基于黏彈性邊界的拱壩地震反應分析方法[J]. 水利學報，2006，37（9）：1063-1069" " Du X L，Zhao M. Analysis method for seismic response of arch dams in time domain based on viscous-spring artificial boundary condition[J]. Journal of Hydraulic Engineering，2006，37（9）：1063-1069

Higdon R L. Absorbing boundary conditions for acoustic and elastic waves in stratified media[J]. Journal of Computational Physics，1992，101（2）：386-418

Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics，1994，114（2）：185-200

Liao Z P，Wong H L. A transmitting boundary for the numerical simulation of elastic wave propagation[J]. Soil Dynamics and Earthquake Engineering，1984，3（4）：174-183

吳翔，丁海平. 基于粘彈性邊界的三維土體地震響應分析[J]. 蘇州科技大學學報（工程技術版），2023，36（4）：1-7" " Wu X，Ding H P. Seismic response analysis of three-dimensional soil based on viscoelastic boundary[J]. Journal of Suzhou University of Science and Technology （Engineering and Technology Edition），2023，36（4）：1-7

秦得順，郭永剛，蘇立彬，等. 基于粘彈性邊界在重力壩動力分析中的應用[J]. 云南水力發電，2024，40（3）：74-77" " Qin D S，Guo Y G，Su L B，et al. Application of viscoelastic boundary in dynamic analysis of gravity dam[J]. Yunnan Water Power，2024，40（3）：74-77

劉喜珠，江琦，張健，等. 考慮管-土耦合的穿管堤防地震動力響應分析[J]. 水電能源科學，2023，41（12）：129-133" " Liu X Z，Jiang Q，Zhang J，et al. Seismic dynamic response analysis of pipe-piercing dike considering pipe-soil coupling[J]. Water Resources and Power，2023，41（12）：129-133

陳志偉，史振華，鐘紅，等. 陡崖地形上重力壩地震響應分析研究[J]. 水利規劃與設計，2023（6）：117-122" " Chen Z W，Shi Z H，Zhong H，et al. Seismic response analysis of gravity dams on steep cliff terrain[J]. Water Resources Planning and Design，2023（6）：117-122

廖振鵬，劉晶波. 波動有限元模擬的基本問題[J]. 中國科學：B輯，1992（8）：874-882" " Liao Z P，Liu J B. Basic problems of wave finite element simulation[J]. Science in China （Series B），1992（8）：874-882

周正華，廖振鵬. 消除多次透射公式飄移失穩的措施[J]. 力學學報，2001，33（4）：550-554" " Zhou Z H，Liao Z P. A measure for eliminating drift instability of the multi-transmitting formula[J]. Acta Mechanica Sinica，2001，33（4）：550-554

廖振鵬，劉晶波. 波動的有限元模擬：基本問題和基本研究方法[J]. 地震工程與工程振動，1989，9（4）：1-14" " Liao Z P，Liu J B. Finite element simulation of wave motion：Basic problems and conceptual aspects[J]. Earthquake Engineering and Engineering Vibration，1989，9（4）：1-14

李小軍，劉愛文. 動力方程求解的顯式積分格式及其穩定性與適用性[J]. 世界地震工程，2000，16（2）：8-12" " Li X J，Liu A W. Explicit step-by-step integration formulas for dynamic differential equation and their stability and applicability[J]. World Information on Earthquake Engineering，2000，16（2）：8-12

周正華，周扣華. 有阻尼振動方程常用顯式積分格式穩定性分析[J]. 地震工程與工程振動，2001，21（3）：22-28" " Zhou Z H，Zhou K H. Stability analysis of explicit integral methods for damped vibration equation[J]. Earthquake Engineering and Engineering Vibration，2001，21（3）：22-28

李小軍. 非線性場地地震反應分析方法的研究[D]. 哈爾濱：中國地震局工程力學研究所，1993" " Li X J. Research on nonlinear site seismic response analysis method[D]. Harbin：Institute of Engineering Mechanics，China Earthquake Administration，1993

周正華，魏景芝，王玉石，等. 修正算子" γB_0^0的物理含義及精度分析[J]. 計算力學學報，2011，28（1）：20-24" " Zhou Z H，Wei J Z，Wang Y S，et al. The physical implication and the precision analysis of modified operator" "γB_0^0[J]. Chinese Journal of Computational Mechanics，2011，28（1）：20-24

唐暉，李小軍. 求解散射問題中消除多次透射邊界飄移失穩措施的效果分析[J]. 地震研究，2019，42（4）：493-502" " Tang H，Li X J. Analysis on measures of eliminating drift instability of multi-transmission formula in solving scattering problems[J]. Journal of Seismological Research，2019，42（4）：493-502

章旭斌，廖振鵬，謝志南. 透射邊界高頻失穩機理及穩定實現：P-SV波動[J]. 地球物理學報，2021，64（10）：3646-3656" " Zhang X B，Liao Z P，Xie Z N. Mechanism of high frequency instability and stable implementation for transmitting boundary：P-SV wave motion[J]. Chinese Journal of Geophysics，2021，64（10）：3646-3656

邢浩潔，劉愛文，李小軍，等. 多人工波速優化透射邊界在譜元法地震波動模擬中的應用[J]. 地震學報，2022，44（1）：26-39" " Xing H J，Liu A W，Li X J，et al. Application of an optimized transmitting boundary with multiple artificial wave velocities in spectral-element simulation of seismic wave propagation[J]. Acta Seismologica Sinica，2022，44（1）：26-39

邢浩潔，李小軍，劉愛文，等. 波動數值模擬中的外推型人工邊界條件[J]. 力學學報，2021，53（5）：1480-1495" " Xing H J，Li X J，Liu A W，et al. Extrapolation-type artificial boundary conditions in the numerical simulation of wave motion[J]. Chinese Journal of Theoretical and Applied Mechanics，2021，53（5）：1480-1495

邢浩潔，李鴻晶，李小軍. 一維波動有限元模擬中透射邊界的時域穩定條件[J]. 應用基礎與工程科學學報，2021，29（3）：617-632" " Xing H J，Li H J，Li X J. Time-domain stability conditions of multi-transmitting formula in one-dimensional finite-element simulation of wave motion[J]. Journal of Basic Science and Engineering，2021，29（3）：617-632

Analysis of the influence of interpolation algorithm on the numerical simulation results of transmission boundary

Duan Xueliang， Zhou Zhenghua*， Bian Zhu， Han Yi， Zhao Ling， Li Zheng， Liao Chengliang， He Jiacong， Liu Wei

College of Transportation Engineering， Nanjing Tech University， Nanjing 210009， China

[Abstract]" " "The distance traveled by contour wave with artificial wave velocity in a step-time is often inconsistent with the mesh size. Therefore， when the transmission boundary formula is applied， the calculation nodes of transmission boundary are usually represented by finite element nodes through the method of time interpolation or spatial interpolation. In this paper， the transmission boundary interpolation formulas of Lagrange interpolation algorithm， Hermite interpolation algorithm and cubic spline interpolation algorithm are derived， and the influence of different interpolation algorithms on the numerical simulation results of transmission boundary is discussed. Firstly， the numerical and analytical solutions of the incident P/SV wave free field problem in elastic half space are compared， and the influence of three interpolation algorithms on the calculation accuracy is discussed. Then， taking the scattering problem of incident P/SV wave in semi-circular valley terrain as an example， the results obtained by different interpolation algorithms are compared. The numerical results show that： Different interpolation algorithms have significant influence on the numerical simulation results of transmission boundary. Among the three interpolation algorithms， the cubic spline interpolation algorithm has the highest accuracy， followed by Lagrange interpolation algorithm， and the Hermite interpolation algorithm has lower accuracy. Compared with the free field， the results of scattering problem are more affected by three interpolation algorithms， and the results of Lagrange interpolation algorithm and cubic spline interpolation algorithm are similar， while the results of Hermite interpolation algorithm are different. The findings of this study can serve as a scientific foundation for the selection of an appropriate interpolation algorithm when utilizing transmission boundaries to simulate the dynamic reaction process of the site.

[Keywords] transmission boundary; interpolation algorithm; free field; scattering field