插值算法對(duì)透射邊界數(shù)值模擬結(jié)果的影響分析

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

[關(guān)鍵詞] 透射邊界; 插值算法; 自由場; 散射場

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

基金項(xiàng)目： 國家自然科學(xué)基金項(xiàng)目（U2039208）資助。

0" 引言

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

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

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

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

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

1" 理論推導(dǎo)

1.1" 透射邊界理論

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

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

式中，C_j^N為二次多項(xiàng)式系數(shù)，離散位移u_j^p表達(dá)式如下：

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）

由于在一個(gè)時(shí)間步內(nèi)外行波以人工波速走過的距離c_a Δt與網(wǎng)格尺寸Δx往往不一致，因此在應(yīng)用式（3）時(shí)，一般通過時(shí)間插值或空間插值的方法將透射邊界的計(jì)算結(jié)點(diǎn)用有限元結(jié)點(diǎn)表示（圖1）。用有限元三節(jié)點(diǎn)插值來表示式（3）中等式右端，計(jì)算表達(dá)式如下：

（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是插值系數(shù)。本文將推導(dǎo)不同插值算法的一階透射公式，著重探討不同插值算法對(duì)數(shù)值計(jì)算精度的影響。

1.2" 插值系數(shù)推導(dǎo)

1.2.1" 插值基本概念

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

1.2.2" 拉格朗日插值

Lagrange插值公式為：

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

據(jù)此可以得到插值多項(xiàng)式：

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）中，可以得到由拉格朗日插值法計(jì)算得到的一次透射公式的插值系數(shù)：

（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）

根據(jù)一階、二階差商公式，可以得到N階差商的表達(dá)式：

（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）

當(dāng)x_i-x_j=Δx為定值時(shí)，

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

由此可得插值多項(xiàng)式：

（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），可以得到由牛頓插值法計(jì)算得到的一次透射公式的插值系數(shù)，如下：

（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），可以發(fā)現(xiàn)，拉格朗日與牛頓插值算法一階透射公式的插值系數(shù)相同，這是由于兩種插值算法雖表達(dá)形式各異，但在本質(zhì)上具有相同的階次和系數(shù)多項(xiàng)式，因此對(duì)于高階透射公式，同樣會(huì)得出一致的插值系數(shù)。

1.2.4" 三次埃爾米特插值

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

（1） G（x）在每個(gè)小區(qū)間上的多項(xiàng)式次數(shù)為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個(gè)節(jié)點(diǎn)xi上的分段三次埃爾米特插值多項(xiàng)式。

根據(jù)埃爾米特插值要求，

（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，進(jìn)而計(jì)算得到插值系數(shù)公式：

（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" 三次樣條插值

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

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

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

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

采用三點(diǎn)插值，可以將整個(gè)區(qū)間劃分為兩段，為[0，Δx]與[Δx，2Δx]。在每個(gè)區(qū)間內(nèi)求解一個(gè)三次方程：

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

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

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

上式也提供了2個(gè)方程。至此可以求解每個(gè)區(qū)間內(nèi)的三次方程，具體方程如下：

（（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）便可得到每個(gè)區(qū)間內(nèi)的三次方程。

考慮到caΔt ≤ Δx，因此透射邊界的計(jì)算結(jié)點(diǎn)落在第一個(gè)區(qū)間內(nèi)，滿足第一個(gè)區(qū)間內(nèi)的三次方程，可得：

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）

聯(lián)合式（19）與式（4），便可得到插值系數(shù)：

（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" 數(shù)值計(jì)算

2.1" 入射波

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

{（（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）為赫維賽德函數(shù)，t是時(shí)間，A是脈沖的峰值，T是脈沖的作用時(shí)間。入射波位移脈沖時(shí)程與傅里葉頻譜如圖2所示。需要說明的是：對(duì)于自由場問題，本文取T=1；對(duì)于散射場問題，本文取T=2。

2.2" 有限元模型

本文數(shù)值模擬所取的材料參數(shù)分別為：泊松比v=0.25，介質(zhì)密度ρ=2000 kg/m3，彈性模型E= 5×109 Pa。根據(jù)式（22）可以計(jì)算得出，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）

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

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

要求一個(gè)時(shí)間步內(nèi)，波傳播的距離不大于一個(gè)網(wǎng)格尺寸，進(jìn)而可以確定時(shí)間步長，計(jì)算公式如下：

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

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

2.3" 彈性半空間

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

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

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

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

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

2.4" 半圓河谷地形

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

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

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

3" 結(jié)論

為了將拉格朗日插值、埃爾米特插值和三次樣條插值算法應(yīng)用到透射邊界中，本文首先導(dǎo)出了不同插值算法的一階透射邊界插值系數(shù)，進(jìn)而溝通了透射節(jié)點(diǎn)和有限元節(jié)點(diǎn)的關(guān)系，為透射理論與有限元理論結(jié)合提供了基礎(chǔ)。隨后，通過比較3種插值算法計(jì)算得出的自由場觀測點(diǎn)位移時(shí)程與理論解析解之間的差異，系統(tǒng)地評(píng)估了各插值方法在計(jì)算精度方面的表現(xiàn)。此外，還深入探討了這3種插值算法在散射場觀測點(diǎn)位移時(shí)程計(jì)算中的差異性，以此來詳盡剖析它們?cè)谕干溥吔鐢?shù)值模擬結(jié)果上的影響程度。得出了如下結(jié)論：

（1） 插值算法對(duì)透射邊界數(shù)值模擬結(jié)果有影響，特別是對(duì)于散射問題。

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

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

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

致謝

本文由國家自然科學(xué)基金項(xiàng)目（U2039208）資助。感謝審稿人的建設(shè)性意見，這些意見極大地改進(jìn)了本文。

參考文獻(xiàn)

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

廖振鵬. 近場波動(dòng)的數(shù)值模擬[J]. 力學(xué)進(jìn)展，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

李小信，何超，周順華，等. 具有不規(guī)則界面的層狀地基三維動(dòng)力響應(yīng)的薄層法[J]. 巖土力學(xué)，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

劉晶波，王振宇，杜修力，等. 波動(dòng)問題中的三維時(shí)域粘彈性人工邊界[J]. 工程力學(xué)，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

杜修力，趙密. 基于黏彈性邊界的拱壩地震反應(yīng)分析方法[J]. 水利學(xué)報(bào)，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

吳翔，丁海平. 基于粘彈性邊界的三維土體地震響應(yīng)分析[J]. 蘇州科技大學(xué)學(xué)報(bào)（工程技術(shù)版），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

秦得順，郭永剛，蘇立彬，等. 基于粘彈性邊界在重力壩動(dòng)力分析中的應(yīng)用[J]. 云南水力發(fā)電，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

劉喜珠，江琦，張健，等. 考慮管-土耦合的穿管堤防地震動(dòng)力響應(yīng)分析[J]. 水電能源科學(xué)，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

陳志偉，史振華，鐘紅，等. 陡崖地形上重力壩地震響應(yīng)分析研究[J]. 水利規(guī)劃與設(shè)計(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

廖振鵬，劉晶波. 波動(dòng)有限元模擬的基本問題[J]. 中國科學(xué)：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

周正華，廖振鵬. 消除多次透射公式飄移失穩(wěn)的措施[J]. 力學(xué)學(xué)報(bào)，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

廖振鵬，劉晶波. 波動(dòng)的有限元模擬：基本問題和基本研究方法[J]. 地震工程與工程振動(dòng)，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

李小軍，劉愛文. 動(dòng)力方程求解的顯式積分格式及其穩(wěn)定性與適用性[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

周正華，周扣華. 有阻尼振動(dòng)方程常用顯式積分格式穩(wěn)定性分析[J]. 地震工程與工程振動(dòng)，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

李小軍. 非線性場地地震反應(yīng)分析方法的研究[D]. 哈爾濱：中國地震局工程力學(xué)研究所，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]. 計(jì)算力學(xué)學(xué)報(bào)，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

唐暉，李小軍. 求解散射問題中消除多次透射邊界飄移失穩(wěn)措施的效果分析[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

章旭斌，廖振鵬，謝志南. 透射邊界高頻失穩(wěn)機(jī)理及穩(wěn)定實(shí)現(xiàn)：P-SV波動(dòng)[J]. 地球物理學(xué)報(bào)，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

邢浩潔，劉愛文，李小軍，等. 多人工波速優(yōu)化透射邊界在譜元法地震波動(dòng)模擬中的應(yīng)用[J]. 地震學(xué)報(bào)，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

邢浩潔，李小軍，劉愛文，等. 波動(dòng)數(shù)值模擬中的外推型人工邊界條件[J]. 力學(xué)學(xué)報(bào)，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

邢浩潔，李鴻晶，李小軍. 一維波動(dòng)有限元模擬中透射邊界的時(shí)域穩(wěn)定條件[J]. 應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，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