SH波斜入射時(shí)有阻尼成層介質(zhì)自由場(chǎng)的一維化時(shí)域算法

高智能， 卓衛(wèi)東，2， 谷 音，2

(1. 福州大學(xué) 土木工程學(xué)院，福州 350116； 2. 福州大學(xué) 福建省土木工程多災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，福州 350116)

SH波斜入射時(shí)有阻尼成層介質(zhì)自由場(chǎng)的一維化時(shí)域算法

高智能1， 卓衛(wèi)東1，2， 谷 音1，2

(1. 福州大學(xué) 土木工程學(xué)院，福州 350116； 2. 福州大學(xué) 福建省土木工程多災(zāi)害防治重點(diǎn)實(shí)驗(yàn)室，福州 350116)

在進(jìn)行大型大跨結(jié)構(gòu)抗震設(shè)計(jì)時(shí)，需要考慮地震波斜入射的影響。在劉晶波等提出的一維化時(shí)域算法基礎(chǔ)上，采用有限差分方法，推導(dǎo)建立了SH波斜入射情形下、考慮阻尼影響的水平成層彈性介質(zhì)出平面自由場(chǎng)求解的顯式數(shù)值逐步法公式，并采用Fortan程序語(yǔ)言編制了相應(yīng)的數(shù)值計(jì)算程序。綜合采用有限元法和數(shù)值逐步法，進(jìn)行SH波斜入射下單層土和雙層土的算例分析。結(jié)果表明：土層介質(zhì)阻尼僅影響到SH波斜入射下自由場(chǎng)的位移幅值，而對(duì)位移時(shí)程的波形沒(méi)有影響；SH波斜入射下有阻尼與無(wú)阻尼情形相比位移幅值有明顯衰減，對(duì)單層土，30°斜入射下衰減幅度可達(dá)40.5%；對(duì)雙層土，最大衰減幅度可達(dá)39.7%；有阻尼土層位移幅值的衰減幅度從底部邊界向自由表面逐漸增大。

成層介質(zhì)；出平面波動(dòng)；斜入射；一維化時(shí)域算法；介質(zhì)阻尼；逐步法

已有研究表明，地震波斜入射引起的地面運(yùn)動(dòng)非一致變化對(duì)地鐵車站、長(zhǎng)大橋梁、大壩等大型大跨結(jié)構(gòu)的地震響應(yīng)有較大的影響[1-4]。因此，有必要在這類結(jié)構(gòu)的抗震設(shè)計(jì)中考慮地震波斜入射的影響。在實(shí)際工程中，由于土層分布的復(fù)雜性，對(duì)地震波斜入射下土層自由場(chǎng)的解析求解幾乎是不可能的[5-6]。因此，采用各種數(shù)值方法，研究復(fù)雜土層介質(zhì)中的地震波動(dòng)問(wèn)題，已經(jīng)成為主要的研究方向。

實(shí)際大型大跨結(jié)構(gòu)基礎(chǔ)所在的土層往往被假定為水平成層介質(zhì)；對(duì)水平成層彈性半空間，地震波斜入射下自由場(chǎng)的計(jì)算可在頻域內(nèi)實(shí)現(xiàn)[7-8]。然而，為獲得地震動(dòng)時(shí)程，頻域方法需要進(jìn)行多次快速Fourier變換，這會(huì)消耗大量的計(jì)算機(jī)存儲(chǔ)空間和計(jì)算時(shí)間。針對(duì)頻域方法的不足，李山有等[9]根據(jù)地震波斜入射情形下水平成層彈性半空間波動(dòng)傳播水平視波速不變且已知的特點(diǎn)，采用透射邊界[10]模擬遠(yuǎn)場(chǎng)對(duì)有限域的影響，建立了入射側(cè)邊界節(jié)點(diǎn)自由場(chǎng)計(jì)算的精確內(nèi)插公式，并將其與計(jì)算內(nèi)節(jié)點(diǎn)位移的顯式差分公式相結(jié)合，得到了地震波斜入射下無(wú)阻尼水平成層彈性半空間人工邊界節(jié)點(diǎn)自由場(chǎng)的時(shí)域計(jì)算方程，可用于近場(chǎng)波動(dòng)數(shù)值模擬中內(nèi)行波場(chǎng)與散射波場(chǎng)的時(shí)域同步模擬。劉晶波等[11-13]進(jìn)一步提出一種一維化的時(shí)域算法，將地震波斜入射下水平成層彈性半空間自由場(chǎng)的計(jì)算問(wèn)題簡(jiǎn)化為時(shí)域內(nèi)的一維問(wèn)題求解，且具有與二維有限元數(shù)值解同樣的精度。趙密等[14]在劉晶波等算法的基礎(chǔ)上，提出一種模擬基巖半空間輻射阻尼的人工邊界條件，采用該人工邊界條件代替黏性邊界條件，提高了大角度斜入射下平面內(nèi)自由波場(chǎng)的計(jì)算精度。為了改善顯式有限差分時(shí)域算法的局限性，范留明等[15]基于惠更斯波動(dòng)觀點(diǎn)，提出一種稱為“界面子波算法”的新的時(shí)域數(shù)值方法，可求解成層彈性半空間出平面一維自由場(chǎng)。王篤國(guó)等[16]基于一維等效線性化頻域方法，利用成層彈性介質(zhì)頻域傳遞矩陣，推導(dǎo)了等效剪應(yīng)變公式，提出一種可計(jì)算地震波斜入射下成層介質(zhì)非線性自由場(chǎng)的數(shù)值方法。

已有關(guān)于地震波斜入射情形下水平成層彈性半空間自由場(chǎng)的顯式有限差分時(shí)域算法，均未考慮土層介質(zhì)阻尼的影響。然而，土層介質(zhì)阻尼可能對(duì)地震波斜入射下土層介質(zhì)的自由場(chǎng)有重要的影響。本文在劉晶波等算法的基礎(chǔ)上，考慮土層介質(zhì)阻尼，綜合采用有限元法與有限差分法，推導(dǎo)建立SH波斜入射下水平成層彈性介質(zhì)出平面自由場(chǎng)求解的數(shù)值逐步法公式，并編制相應(yīng)的數(shù)值計(jì)算程序；通過(guò)算例分析，討論土層介質(zhì)阻尼對(duì)SH波斜入射情形下水平成層彈性介質(zhì)出平面自由場(chǎng)的影響。

1 成層半空間出平面波動(dòng)問(wèn)題的一維化

在地下土層自由場(chǎng)的計(jì)算中，通常將土層介質(zhì)假設(shè)為由下臥彈性半空間(基巖)和其上的水平成層彈性介質(zhì)(土層)組成，各層的質(zhì)量密度和剪切波速分別為ρl和csl(l=1,2,…,L)。按照劉晶波等提出的一維化算法的思路，將水平成層彈性半空間劃分為如圖1所示的有限元離散化模型；圖1中，豎向網(wǎng)格尺寸為Δy，可取為滿足精度要求的任意值；水平方向網(wǎng)格尺寸為Δx，其取值需滿足式(1)

Δx=cx·Δt

(1)

式中：cx為入射SH波的水平視波速；Δt為時(shí)間步長(zhǎng)。

圖1 水平成層半空間有限元模型Fig.1 Finite element model of the layered half-space

(2)

式中，Mi,j、Ci,j和Ki,j分別為與節(jié)點(diǎn)(i,j)相關(guān)的集中質(zhì)量系數(shù)、阻尼系數(shù)和剛度系數(shù)。

由式(2)可見(jiàn)，節(jié)點(diǎn)(m,n)的運(yùn)動(dòng)只與該節(jié)點(diǎn)直接相鄰的節(jié)點(diǎn)運(yùn)動(dòng)有關(guān)。在彈性介質(zhì)中，沿水平方向的波動(dòng)為行波，波速為視波速cx，因此，彈性介質(zhì)內(nèi)的波場(chǎng)可以表示為

(3)

若取時(shí)間步長(zhǎng)為Δt，則由式(3)可以得到式(4)

u(x+jcxΔt,y,t)=u(x,y,t-jΔt) j∈(-∞,∞)

(4)

對(duì)圖1中的各離散節(jié)點(diǎn)，利用式(4)可得到

(5)

當(dāng)i=±1時(shí)，由式(5)有

(6)

同理，節(jié)點(diǎn)速度和加速度同樣存在以上關(guān)系。由式(6)可見(jiàn)，與節(jié)點(diǎn)(m,n)直接相鄰的同一層節(jié)點(diǎn)(m-1,n)和節(jié)點(diǎn)(m+1,n)在pΔt時(shí)刻的運(yùn)動(dòng)可用節(jié)點(diǎn)(m,n) 在相鄰時(shí)刻的運(yùn)動(dòng)來(lái)表示。

綜合式(2)和式(6)可見(jiàn)，節(jié)點(diǎn)(m,n)在pΔt時(shí)刻的運(yùn)動(dòng)方程式(2)中，未知量?jī)H包含第m列節(jié)點(diǎn)在pΔt及其相鄰時(shí)刻的運(yùn)動(dòng)。在求得第m列節(jié)點(diǎn)位移之后，根據(jù)式(3)即可確定整個(gè)半空間中的自由場(chǎng)。如此，SH波斜入射下水平成層彈性半空間出平面自由場(chǎng)的計(jì)算問(wèn)題就轉(zhuǎn)化為一維問(wèn)題求解。

2 考慮介質(zhì)阻尼的一維化時(shí)域算法

2.1 考慮介質(zhì)阻尼的數(shù)值逐步法的計(jì)算公式

首先，采用劉晶波等算法的思路，將SH波斜入射下水平成層彈性半空間出平面自由場(chǎng)的計(jì)算問(wèn)題轉(zhuǎn)化為一維問(wèn)題。其次，綜合采用有限元法與有限差分法，在時(shí)域內(nèi)直接求解式(2)所列的運(yùn)動(dòng)方程。在數(shù)值計(jì)算中可取m=0，即先求得y軸上各節(jié)點(diǎn)的位移，再根據(jù)式(3)，依次確定水平成層彈性半空間中的自由場(chǎng)。以下，采用有限差分方法，推導(dǎo)建立考慮土層介質(zhì)阻尼時(shí)求解y軸上各節(jié)點(diǎn)位移的計(jì)算列式。

根據(jù)式(2)，在考慮土層介質(zhì)阻尼時(shí)，y軸上任一內(nèi)節(jié)點(diǎn)(0,n)在pΔt時(shí)刻的運(yùn)動(dòng)方程為

(n=1,2,…,N-1)

(7)

式中，M0,n為節(jié)點(diǎn)(0,n)的集中質(zhì)量。

y軸與自由表面交界處的節(jié)點(diǎn)(0,0) 在pΔt時(shí)刻的運(yùn)動(dòng)方程為

(8)

式中，M0,0為節(jié)點(diǎn)(0,0)處的集中質(zhì)量。

在采用數(shù)值方法計(jì)算水平成層彈性半空間的自由場(chǎng)時(shí)，顯然必須從半無(wú)限介質(zhì)中切取有限的計(jì)算區(qū)域，并在區(qū)域邊界引入合適的人工邊界條件。目前已發(fā)展了多種人工邊界[17-26]，不失一般性，本文假定人工邊界為黏性邊界。對(duì)y軸與黏性邊界交界處的節(jié)點(diǎn)(0,N)(設(shè)為節(jié)點(diǎn)B)，其在pΔt時(shí)刻的運(yùn)動(dòng)方程為

(9)

(10)

式中：ρ為介質(zhì)質(zhì)量密度；cs為介質(zhì)剪切波速；θ為地震波入射角；w0(xB,yB,t)和τ0(xB,yB,t)分別為入射波在黏性邊界節(jié)點(diǎn)B上產(chǎn)生的位移和應(yīng)力。

式(7)～式(9)所列的運(yùn)動(dòng)方程中的加速度項(xiàng)可利用中心差分法近似計(jì)算

(11)

(12)

考慮式(6)，同時(shí)將式(11)和式(12)代入式(7)～式(9)中，并簡(jiǎn)記u0,n=un，M0,n=Mn，經(jīng)整理得到如下矩陣形式的方程

(13)

(14a)

(14b)

從式(14)可見(jiàn)，式(13)所列的矩陣方程左邊的系數(shù)矩陣是稀疏的三對(duì)角矩陣，右邊的向量?jī)H與邊界節(jié)點(diǎn)輸入的等效荷載以及y軸上各節(jié)點(diǎn)在pΔt時(shí)刻及其前一時(shí)刻(p-1)Δt的位移有關(guān)，故只要給定邊界節(jié)點(diǎn)等效荷載以及初始時(shí)刻各節(jié)點(diǎn)的位移和速度值，即可通過(guò)逐步法求解上述方程組，得到軸上各節(jié)點(diǎn)在時(shí)域內(nèi)的位移解，進(jìn)而確定全部自由場(chǎng)。

2.2 計(jì)算步驟及程序?qū)崿F(xiàn)

由式(13)所列的矩陣方程可見(jiàn)，本文基于有限差分方法建立的考慮介質(zhì)阻尼的水平成層彈性介質(zhì)自由場(chǎng)的數(shù)值逐步法的計(jì)算公式是顯式的，因此，求解容易實(shí)現(xiàn)。具體計(jì)算步驟如下：

步驟1 對(duì)選取的計(jì)算區(qū)域進(jìn)行離散化(見(jiàn)圖1)，采用有限元法計(jì)算離散化模型的質(zhì)量矩陣M和剛度矩陣K，其中質(zhì)量矩陣M采用集中質(zhì)量法計(jì)算，并假定阻尼矩陣C。

步驟2 根據(jù)黏性人工邊界條件，確定人工邊界節(jié)點(diǎn)等效荷載時(shí)程。

(15)

步驟6 由式(3)，確定計(jì)算區(qū)域內(nèi)各節(jié)點(diǎn)的自由場(chǎng)。

根據(jù)上述計(jì)算步驟，基于Fortran語(yǔ)言編制了相應(yīng)的計(jì)算程序。限于篇幅，這里沒(méi)有給出計(jì)算程序的源代碼。

2.3 算法的穩(wěn)定性和計(jì)算精度

本文所建立的SH波斜入射下考慮介質(zhì)阻尼的水平成層彈性介質(zhì)出平面自由場(chǎng)的一維化時(shí)域算法，其穩(wěn)定性條件仍與劉晶波等算法的穩(wěn)定性條件相同，即要求

(16)

文獻(xiàn)[11]分析表明，當(dāng)時(shí)間步長(zhǎng)滿足式(16)所表示的穩(wěn)定性條件時(shí)，一維化時(shí)域算法對(duì)于不同方向的入射波都有很高的計(jì)算精度。

3 算例分析

以下選擇物理性質(zhì)均勻的單一土層(下文簡(jiǎn)稱單層土)和基巖表面覆蓋軟土層(下文簡(jiǎn)稱雙層土)這兩種情形分別進(jìn)行分析。算例分析中，均假定土層和基巖介質(zhì)阻尼為Rayleigh阻尼。因而，可按下式計(jì)算其阻尼矩陣C[27]

C=αM+βK

(17)

式中，α、β為比例系數(shù)，它們由式(18)計(jì)算

(18)

式中：ωi和ωj分別為兩個(gè)特定的自振圓頻率；ξi和ξj分別為與ωi和ωj相關(guān)的模態(tài)阻尼比。

3.1 單層土算例分析

假定單層土由半無(wú)限空間的均勻彈性土層組成，其幾何參數(shù)和物理參數(shù)取值如表1所列，計(jì)算區(qū)域yb取為100 m。計(jì)算區(qū)域的有限元離散化模型如圖2所示，其中豎向網(wǎng)格尺寸Δy取為5 m，水平方向網(wǎng)格尺寸Δx根據(jù)式(1)確定。在計(jì)算區(qū)域的底部邊界引入人工黏性邊界條件，并假定在黏性人工邊界處入射SH波為持時(shí)0.5 s、峰值1 cm的Dirac函數(shù)形式的單位脈沖，其位移時(shí)程如圖3所示。在數(shù)值計(jì)算中，為滿足穩(wěn)定性條件，時(shí)間步長(zhǎng)取為0.005 s，持時(shí)取為2.5 s。

表1 單層土模型參數(shù)

圖2 單層土有限元離散化模型Fig.2 FEM model of one-layer soil

圖3 黏性邊界節(jié)點(diǎn)C1處入射的SH波位移時(shí)程Fig.3 Displacement time history of incident SH wave at boundary joint C1

利用ABAQUS軟件建立計(jì)算區(qū)域的有限元模型，計(jì)算得到單層土模型的質(zhì)量矩陣M和剛度矩陣K；通過(guò)動(dòng)力特性分析，得到該單層土模型前兩階的自振頻率：ω1=8.864 5rad/s，ω2=10.028 0rad/s；將其代入式(18)，并取單層土模態(tài)阻尼比ξ1=ξ2=ξ=0.1，求得α=0.941 2，β=0.010 58；將α、β的數(shù)值代入式(17)，得到單層土模型的阻尼矩陣C。

利用所編制的計(jì)算程序，計(jì)算了SH波以30°斜入射情形下單層土模型的自由場(chǎng)。選取y軸上自由表面節(jié)點(diǎn)A1、中部節(jié)點(diǎn)B1以及底部邊界靠?jī)?nèi)側(cè)節(jié)點(diǎn)C1作為觀測(cè)點(diǎn)(見(jiàn)圖2)，對(duì)計(jì)算結(jié)果進(jìn)行分析討論。由于無(wú)阻尼情形下一維化算法的計(jì)算結(jié)果已得到驗(yàn)證，這里主要討論考慮土層介質(zhì)阻尼情形下一維化算法的計(jì)算結(jié)果。

圖4繪出了SH波以30°斜入射情形下，采用一維化算法計(jì)算得到的無(wú)阻尼和有阻尼單層土模型在3個(gè)觀測(cè)點(diǎn)處的位移響應(yīng)。從圖4中可以發(fā)現(xiàn)，無(wú)阻尼與有阻尼均勻單層土的位移時(shí)程曲線波形一致；有阻尼與無(wú)阻尼均勻單層土相比，前者在各個(gè)時(shí)刻的位移響應(yīng)絕對(duì)值均減小了，尤其在自由表面觀測(cè)點(diǎn)A1處，位移幅值衰減幅度最大。這說(shuō)明土層介質(zhì)阻尼僅影響到位移幅值，而對(duì)位移時(shí)程的波形沒(méi)有影響。

(a) 觀測(cè)點(diǎn)A1位移時(shí)程的計(jì)算結(jié)果

(b) 觀測(cè)點(diǎn)B1位移時(shí)程的計(jì)算結(jié)果

(c) 觀測(cè)點(diǎn)C1位移時(shí)程的計(jì)算結(jié)果

圖4 SH波30°斜入射時(shí)無(wú)阻尼和有阻尼單層土的位移響應(yīng)比較

Fig.4 Comparison between the displacement response of the one-layer soil with and without damping by incident SH wave at 30°

根據(jù)計(jì)算結(jié)果，無(wú)阻尼與有阻尼均勻單層土相比，自由表面觀測(cè)點(diǎn)A1處的位移幅值從2.00 cm減小至1.19 cm，衰減幅度達(dá)40.5%；中部觀測(cè)點(diǎn)B1的位移幅值從1.08 cm減小至0.82 cm，衰減幅度達(dá)24.1%；底部觀測(cè)點(diǎn)C1的位移幅值從1.00 cm衰減至0.92 cm，衰減幅度達(dá)8.0%?？梢?jiàn)，SH波斜入射下土層介質(zhì)阻尼顯著影響到均勻單層土自由場(chǎng)的位移幅值，且其影響程度從底部邊界向自由表面逐漸增大。

3.2 雙層土算例分析

假定雙層土由下臥半無(wú)限空間的基巖和上覆50 m厚的均勻土層組成，其幾何參數(shù)和物理參數(shù)取值如表2所列，計(jì)算區(qū)域yb取為100 m。計(jì)算區(qū)域的有限元離散化模型如圖5所示，其中豎向網(wǎng)格尺寸Δy取為5 m，水平方向網(wǎng)格尺寸Δx根據(jù)式(1)確定。在計(jì)算區(qū)域的底部邊界引入人工黏性邊界條件，同樣假定在黏性人工邊界處入射SH波為持時(shí)0.5 s、峰值1 cm的Dirac函數(shù)形式的單位脈沖，其位移時(shí)程見(jiàn)圖3。為滿足穩(wěn)定性條件，時(shí)間步長(zhǎng)取為0.005 s，持時(shí)取為2.5 s。

首先，利用ABAQUS軟件建立計(jì)算區(qū)域的有限元模型，計(jì)算得到雙層土模型的質(zhì)量矩陣M和剛度矩陣K；進(jìn)而得到上覆均勻土層和下臥基巖層的質(zhì)量子矩陣MS和MR，以及剛度子矩陣KS和KR。

表2 雙層土模型參數(shù)

圖5 雙層土有限元離散化模型Fig.5 FEM model of two-layer soil

對(duì)雙層土模型，可假定均勻土層和基巖介質(zhì)阻尼為Rayleigh阻尼；然而，考慮到土層和基巖層阻尼機(jī)制不同，雙層土模型的阻尼矩陣C需要采用非比例阻尼矩陣構(gòu)造。通過(guò)動(dòng)力特性分析，得到該雙層土模型前兩階的自振頻率：ω1=13.585 0 rad/s，ω2=16.258 0 rad/s；取上覆均勻土層的模態(tài)阻尼比ξS1=ξS2=ξS=0.1，取下臥基巖層的模態(tài)阻尼比ξR1=ξR2=ξR=0.03，分別將其代入式(18)，求得均勻土層和基巖層的比例系數(shù)分別為αs=1.480 3，βs=0.006 7，以及αR=0.444 2，βR=0.002 0；將質(zhì)量子矩陣MS、MR和剛度子矩陣KS、KR以及各自的比例系數(shù)代入式(17)，得到上覆均勻土層和下臥基巖層各自的阻尼矩陣CS和CR；通過(guò)有限元直接集裝方法，得到雙層土離散化模型的阻尼矩陣C。上述過(guò)程均可通過(guò)ABAQUS軟件實(shí)現(xiàn)。

在得到雙層土模型的質(zhì)量矩陣M、剛度矩陣K和阻尼矩陣C后，利用所編制的計(jì)算程序，計(jì)算了SH波以不同入射角入射情形下雙層土模型的自由場(chǎng)；所取的入射角變化范圍為0°～90°，并以15°為間隔。選取y軸上自由表面節(jié)點(diǎn)A2、中部節(jié)點(diǎn)B2以及底部邊界靠?jī)?nèi)側(cè)節(jié)點(diǎn)C2作為觀測(cè)點(diǎn)(見(jiàn)圖5)，對(duì)計(jì)算結(jié)果進(jìn)行分析討論。同樣地，由于無(wú)阻尼情形下一維化算法的計(jì)算結(jié)果已得到驗(yàn)證，這里也主要討論考慮土層介質(zhì)阻尼情形下一維化算法的計(jì)算結(jié)果。

圖6繪出了SH波以30°角斜入射情形下，采用一維化算法計(jì)算得到的無(wú)阻尼和有阻尼雙層土模型在3個(gè)觀測(cè)點(diǎn)處的位移響應(yīng)。從圖6中可以發(fā)現(xiàn)，有阻尼雙層土模型中均勻土層和基巖層的位移時(shí)程曲線與無(wú)阻尼雙層土模型的相比，兩者波形一致；有阻尼與無(wú)阻尼雙層土模型相比，前者在各個(gè)時(shí)刻的位移響應(yīng)絕對(duì)值均減小了，尤其是在自由表面觀測(cè)點(diǎn)A2處，位移幅值衰減幅度最大。這說(shuō)明基巖與土層介質(zhì)阻尼也僅影響到位移幅值，而對(duì)位移時(shí)程的波形沒(méi)有影響。

根據(jù)計(jì)算結(jié)果，無(wú)阻尼與有阻尼雙層土模型相比，自由表面觀測(cè)點(diǎn)A2處的位移幅值從2.88cm減小至1.98cm，衰減幅度達(dá)31.3%；中部觀測(cè)點(diǎn)B2的位移幅值從1.48cm減小至1.33cm，衰減幅度達(dá)10.1%；底部觀測(cè)點(diǎn)C2的位移幅值從1.28cm衰減至1.20cm，衰減幅度達(dá)6.3%?？梢?jiàn)，SH波斜入射下土層介質(zhì)阻尼顯著影響到雙層土自由場(chǎng)的位移幅值，且其影響程度同樣從底部邊界向自由表面逐漸增大。

圖7繪出了SH波以不同入射角入射情形下，無(wú)阻尼和有阻尼雙層土在表面觀測(cè)點(diǎn)A2處的位移幅值隨入射角的變化情況。從圖7中可以發(fā)現(xiàn)，無(wú)阻尼和有阻尼雙層土的位移幅值均隨著SH波入射角的增大而減??；在入射角為90°時(shí)，位移幅值減小為0，這也驗(yàn)證了SH波沿底部邊界水平方向入射時(shí)不會(huì)在自由表面上引起位移的理論認(rèn)識(shí)。此外，與無(wú)阻尼雙層土模型相比，相同入射角下有阻尼雙層土模型的位移幅值均減小了。計(jì)算結(jié)果表明，在SH波入射角為75°時(shí)，有阻尼雙層土模型在自由表面觀測(cè)點(diǎn)A2處的位移幅值從1.89cm(無(wú)阻尼時(shí))減小至1.14cm，衰減幅度最大，達(dá)39.7%。從衰減幅值看，在SH波入射角為0°～45°時(shí)，有阻尼雙層土模型的位移幅值皆減小了0.90cm；在SH波入射角超過(guò)45°后，位移幅值相對(duì)衰減較小。

(a) 觀測(cè)點(diǎn)A1位移時(shí)程的計(jì)算結(jié)果

(b) 觀測(cè)點(diǎn)B1位移時(shí)程的計(jì)算結(jié)果

(c) 觀測(cè)點(diǎn)C1位移時(shí)程的計(jì)算結(jié)果

圖7 觀測(cè)點(diǎn)A2位移幅值隨SH波入射角的變化曲線Fig.7 Peak displacements at joint A2by incident SH wave with different angles

4 結(jié) 論

本文基于劉晶波等的算法，提出了考慮土層介質(zhì)阻尼影響的水平成層彈性介質(zhì)出平面自由場(chǎng)的一維化時(shí)域算法，并編制了相應(yīng)的數(shù)值計(jì)算程序。通過(guò)算例分析，討論了土層介質(zhì)阻尼對(duì)SH波斜入射下水平成層彈性介質(zhì)出平面自由場(chǎng)的影響。分析結(jié)果表明：

(1)介質(zhì)阻尼對(duì)水平成層彈性介質(zhì)出平面的自由場(chǎng)有重要的影響。SH波斜入射下有阻尼與無(wú)阻尼情形相比位移幅值有明顯衰減，對(duì)單層土，30°斜入射下衰減幅度可達(dá)40.5%；對(duì)雙層土，最大衰減幅度可達(dá)39.7%。

(2)土層介質(zhì)阻尼僅影響到SH波斜入射下自由場(chǎng)的位移幅值，而對(duì)位移時(shí)程的波形沒(méi)有影響。

(3)不管是單層土還是雙層土，與無(wú)阻尼土層相比，有阻尼土層位移幅值的衰減幅度從底部邊界向自由表面逐漸增大。

[ 1 ] WOLF J P, OBERNHUBER P. Effects of horizontally traveling waves in soil-structure interaction[J]. Nuclear Engineering and Design, 1979, 57(2): 221-244.

[ 2 ] 李山有, 廖振鵬, 周正華. 大型結(jié)構(gòu)地震反應(yīng)數(shù)值模擬中的波動(dòng)輸入[J]. 地震工程與工程振動(dòng), 2001, 6(2): 1-5. LI Shanyou, LIAO Zhenpeng, ZHOU Zhenghua. Wave motion input in numerical simulation of seismic response for 1arge-scale structure[J]. Earthquake Engineering and Engineering Vibration, 2001, 6(2): 1-5.

[ 3 ] 潘旦光, 樓夢(mèng)麟, 范立礎(chǔ). 多點(diǎn)輸入下大跨度結(jié)構(gòu)地震反應(yīng)分析研究現(xiàn)狀[J]. 同濟(jì)大學(xué)學(xué)報(bào), 2001, 29(10): 1213-1219. PAN Dan’guang, LOU Menglin, FAN Lichu. Status of seismic response analysis of long-span structures under multiple support excitations[J]. Journal of Tongji University, 2001, 29(10): 1213-1219.

[ 4 ] 杜修力, 陳維, 李亮, 等. 斜入射條件下地下結(jié)構(gòu)時(shí)域地震反應(yīng)分析初探[J]. 震災(zāi)防御技術(shù), 2007, 2(3): 290-296. DU Xiuli, CHEN Wei, LI Liang, et al. Preliminary study of time-domain seismic response for underground structures to obliquely incident seismic waves[J]. Technology for Earthquake Disaster Prevention, 2007, 2(3): 290-296.

[ 5 ] 傅淑芳, 劉寶誠(chéng). 地震學(xué)教程[M]. 北京: 地震出版社, 1991.

[ 6 ] 廖振鵬. 工程波動(dòng)理論導(dǎo)論[M]. 2版.北京: 科學(xué)出版社, 2002.

[ 7 ] HASKELL N A. The dispersion of surface waves in multi-layered media[J]. Bulletin of the Seismic Society of American, 1953, 43: 17-34.

[ 8 ] BREKHOVSHIKH L M. Waves in layered media[M]. New York: Academic Press Inc., 1980.

[ 9 ] 李山有, 王學(xué)良, 周正華. 地震波斜入射情形下水平成層半空間自由場(chǎng)的時(shí)域計(jì)算[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 2003, 33(3): 372-376. LI Shanyou, WANG Xueliang, ZHOU Zhenghua. The time-step numerical simulation of free field motion of layered half-space for inclined seismic waves[J]. Journal of Jilin University(Earth Science), 2003, 33(3): 372-376.

[10] 廖振鵬, 黃孔亮, 楊柏坡, 等. 暫態(tài)波透射邊界[J]. 中國(guó)科學(xué)(A輯), 1984(6): 556-564. LIAO Zhenpeng, HUANG Kongliang, YANG Baipo, et al. A transmitting boundary for transient wave[J]. Science in China(Series A) , 1984(6): 556-564.

[11] 劉晶波, 王艷. 成層半空間出平面自由波場(chǎng)的一維化時(shí)域算法[J]. 力學(xué)學(xué)報(bào), 2006, 38(2): 219-225. LIU Jingbo, WANG Yan. A 1D time-domain method for 2D wave motion in elastic layered half-space by antiplane wave oblique incidence[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(2): 219-225.

[12] 劉晶波, 王艷. 成層介質(zhì)中平面內(nèi)自由波場(chǎng)的一維化時(shí)域算法[J]. 工程力學(xué), 2007, 24(7): 16-22. LIU Jingbo, WANG Yan. A 1D time-domain method for in-plane wave motion of free field in layered media[J]. Engineering Mechanics, 2007, 24(7): 16-22.

[13] LIU Jingbo, WANG Yan. A 1D time-domain method for inplane wave motions in a layered half-space[J]. Acta Mechanica Sinica, 2007, 23(6)：673-680.

[14] 趙密, 杜修力, 劉晶波, 等. P-SV波斜入射時(shí)成層半空間自由場(chǎng)的時(shí)域算法[J]. 地震工程學(xué)報(bào), 2013, 35(1): 84-90. ZHAO Mi, DU Xiuli, LIU Jingbo, et al. Time-domain method for free field in layered half space under P-SV waves of oblique incidence[J]. China Earthquake Engineering Journal, 2013, 35(1): 84-90.

[15] 范留明, 趙欽. 成層半空間出平面自由波場(chǎng)的界面子波算法[J]. 地震工程學(xué)報(bào), 2014, 36(4): 765-772. FAN Liuming, ZHAO Qin. Interfacial wavelet superposition method for out-of-plane wave motion of a free field in a layered half space[J]. China Earthquake Engineering Journal, 2014, 36(4): 765-772.

[16] 王篤國(guó), 趙成剛. 地震波斜入射時(shí)二維成層介質(zhì)自由場(chǎng)求解的等效線性化方法[J]. 巖土工程學(xué)報(bào), 2016, 38(3): 554-561. WANG Duguo, ZHAO Chenggang. Two-dimensional equivalent linear seismic analysis of free field in layered half-space due to oblique incidence[J]. Chinese Journal of Geotechnical Engineering, 2016, 38(3): 554-561.

[17] LYSMER J, KULEMEYER R L. Finite dynamic model for infinite media[J]. Journal of Engineering Mechanics, 1969, 95(4): 859-877.

[18] WOLF P J. A comparison of time-domain transmitting boundaries[J]. Earthquake Engineering and Structure Dynamics, 1986, 14(4):655-673.

[19] DEEKS A J, RANDOLPH M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of Engineering Mechanics, 1994, 120(1): 25-42.

[20] 王振宇, 劉晶波. 成層地基非線性波動(dòng)問(wèn)題人工邊界與波動(dòng)輸入研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 2004, 23(7): 1169-1173. WANG Zhenyu, LIU Jingbo. Study on wave motion input and artificial boundary for problem of nonlinear wave motion in layered soil[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(7): 1169-1173.

[21] 劉晶波, 王振宇, 杜修力, 等. 波動(dòng)問(wèn)題中的三維時(shí)域黏彈性人工邊界[J]. 工程力學(xué), 2005, 22(6): 46-51. LIU Jingbo, WANG Zhenyu, DU Xiuli, et al. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems[J]. Engineering Mechanics, 2005, 22(6): 46-51.

[22] 劉晶波, 谷音, 杜義欣. 一致黏彈性人工邊界及黏彈性邊界單元[J]. 巖土工程學(xué)報(bào), 2006, 28(9): 1070-1075. LIU Jingbo, GU Yin, DU Yixin. Consistent viscous-spring artificial boundaries and viscous-spring boundary elements[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(9): 1070-1075.

[23] 杜修力, 趙密, 王進(jìn)廷. 近場(chǎng)波動(dòng)模擬的人工應(yīng)力邊界條件[J]. 力學(xué)學(xué)報(bào), 2006, 38(1): 49-56. DU Xiuli, ZHAO Mi, WANG Jinting. A stress artificial boundary in FEA for near-field wave problem[J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(1): 49-56.

[24] 盧華喜, 梁平英, 尚守平. 地基非線性波動(dòng)問(wèn)題中黏-彈性人工邊界研究[J]. 巖土力學(xué), 2008, 29(7): 1911-1916. LU Huaxi, LIANG Pingying, SHANG Shouping. Research on viscoelastic artificial boundary for problem of nonlinear wave motion in soil[J]. Rock and Soil Mechanics, 2008, 29(7): 1911-1916.

[25] 張波, 李術(shù)才, 楊學(xué)英, 等. 三維黏彈性介質(zhì)人工邊界研究[J]. 巖土力學(xué), 2009, 30(11): 3469-3475. ZHANG Bo, LI Shucai, YANG Xueying, et al. Study of three dimensional viscoelastic medium artificial boundary[J]. Rock and Soil Mechanics, 2009, 30(11): 3469-3475.

[26] 趙密, 杜修力, 劉晶波. 一種高階精度人工邊界條件：出平面外域波動(dòng)問(wèn)題[J]. 工程力學(xué), 2012, 29(4): 7-14. ZHAO Mi, DU Xiuli, LIU Jingbo. A high-order accurate artificial boundary condition: out-of-plane exterior waver problem[J]. Engineering Mechanics, 2012, 29(4): 7-14.

[27] CHOPRA A K. Dynamics of structures: theory and applications to earthquake engineering [M]. 3rd ed. New Jersey: Pearson Education Inc., 2007.

A 1D time-domain method for free field motion in layered mediawith damping under obliquely incident SH wave

GAO Zhineng1, ZHUO Weidong1,2, GU Yin1,2

(1. College of Civil Engineering, Fuzhou University, Fuzhou 350116, China; 2. Key Laboratory for Multi Disaster Prevention andGovernance of Civil Engineering of Fujian Province , Fuzhou University , Fuzhou 350116, China )

The influence of obliquely incident seismic wave on the seismic responses of large and long-span structures need to be considered in seismic design. Based on a 1D time-domain method for the out-of-plane wave motion of the free field in a layered half space proposed by Liu and Wang, formulas of an explicit time-stepping method to solve the out-of-plane wave motion of the free field in layered media with damping under obliquely incident SH wave were established by using the finite difference method, and its numerical program was developed with Fortran programming language. A cases study of out-of-plane wave motions of the free field in one-layer soil and two-layer soil with damping under obliquely incident SH wave were carried out respectively. The results show that: soil damping only has effect on the peak displacements of free field motions, and has no effect on the waveforms of the displacement time histories; the displacement amplitudes in layered soil with damping reduce significantly compared with that of without damping under obliquely incident SH wave, and the displacement amplitudes can be reduced by 40.5% and 39.7% for one-layer soil and two-layer soil, respectively. The case study also indicates that the attenuation of displacement amplitudes in layered soil with damping gradually increased from the bottom boundary to the free surface.

layered media; out-of-plane wave motion; oblique incidence; 1D time-domain method; damping; time-stepping method

國(guó)家自然科學(xué)基金資助項(xiàng)目( 51108088)

2016-01-29 修改稿收到日期： 2016-07-01

高智能 男，博士生，1984年生

卓衛(wèi)東 男，博士，教授，1966年生

TU311.3; P315.9

A

10.13465/j.cnki.jvs.2017.16.006