出平面剪切線源作用下沉積盆地地震響應(yīng)譜元法模擬

韓建彬，劉中憲，羅兆輝

(天津城建大學(xué) 土木工程學(xué)院，天津 300384)

土木工程

出平面剪切線源作用下沉積盆地地震響應(yīng)譜元法模擬

韓建彬，劉中憲，羅兆輝

(天津城建大學(xué) 土木工程學(xué)院，天津 300384)

由于盆地內(nèi)外介質(zhì)的阻抗差別及特殊的幾何形狀，盆地內(nèi)部將出現(xiàn)顯著的地震動放大效應(yīng)．針對出平面線源作用下沉積盆地地震響應(yīng)問題，利用譜元法研究了震源主頻率、空間位置和沉積盆地形狀等因素對沉積盆地地震響應(yīng)的影響．?dāng)?shù)值分析表明：點源的頻率和沉積形狀對地表的響應(yīng)影響較大；點源的不同位置對于邊緣效應(yīng)有一定影響．

沉積盆地；出平面剪切線源；譜元法；地震響應(yīng)

國內(nèi)外多次大地震震害調(diào)查表明，沉積盆地內(nèi)部往往會出現(xiàn)烈度異常區(qū)域．這主要由于盆地內(nèi)外介質(zhì)的阻抗差別及其特殊的幾何形狀，使得波動能量容易進入，但很難再向外部逸散，地震波會在盆地邊緣發(fā)生面波轉(zhuǎn)換，并在內(nèi)部某些區(qū)域產(chǎn)生聚焦效應(yīng)，最終使得地震動幅值增大，持時變長，對高層建筑、大跨橋梁等自振周期較長的工程結(jié)構(gòu)造成嚴(yán)重破壞[1]．如1985年墨西哥的8.1級地震，由于軟土盆地效應(yīng)，使得距震中400,km外的墨西哥城遭受嚴(yán)重破壞．我國也有許多大中城市位于沉積盆地之上，或者在城市中分布著古河道，或者曾經(jīng)經(jīng)歷過大體積土體的開挖及回填．盆地效應(yīng)的研究對于這些城市的防震減災(zāi)具有現(xiàn)實意義，同時也對認(rèn)識盆地這種大型不均勻的場地條件對地震動的影響、促進地震工程學(xué)的發(fā)展具有重要意義．

目前，已有幾種較為成熟的數(shù)值方法用于模擬地震波在沉積盆地中的傳播和散射，包括有限元方法[2]、邊界元法[3]、有限差分算法[4-5]、偽譜法以及譜元法[6-8]．

在上述方法當(dāng)中，有限元法可以模擬各種復(fù)雜介質(zhì)模型．該方法基于波動方程弱形式，能夠自動滿足自由邊界條件，同時在解決非均勻復(fù)雜彈性介質(zhì)模型時也很有優(yōu)勢[9]；然而有限元法需要比較強大的計算機資源，以及更多的計算時間，同時，低階有限元法會出現(xiàn)令人不滿意的頻散現(xiàn)象[10]，而傳統(tǒng)的高階有限元法會產(chǎn)生虛假波．對于有限差分方法而言，在處理有規(guī)則的邊界模型及介質(zhì)速度或密度差別不是很大時的模型有優(yōu)勢[11-12]，但是很難有效處理自由界面和復(fù)雜的模型．在流體力學(xué)中的偽譜法已用于解決彈性動力學(xué)問題[13]，這種方法是用矩形網(wǎng)格到彎曲網(wǎng)格的映射算法來處理起伏的自由邊界條件，因此需要大量的計算機資源．邊界積分和邊界元是另外兩種可用于模擬自由起伏邊界的數(shù)值方法，雖然這些方法計算效率高，但多限于線性均勻介質(zhì)問題．另外需指出的是，上述研究大多假定為平面波入射，對于近場情況，可能會產(chǎn)生較大誤差，因此宜采用點源(線源)解作為初步近似．

譜元法最早是由Patera[14]提出的，并用于流體力學(xué)中．該方法將有限元法和譜方法相結(jié)合，既有可以模擬任意復(fù)雜介質(zhì)模型的優(yōu)點，同時也保證了準(zhǔn)確度．筆者針對出平面剪切線源作用下沉積盆地二維響應(yīng)問題，采用譜元方法進行計算分析，討論沉積盆地形狀、不同震動頻率、出平面剪切線源不同位置等因素對沉積盆地地震響應(yīng)規(guī)律的影響．

1 計算模型

以出平面剪切線源為地震輸入，研究任意形狀二維沉積盆地的地震響應(yīng)，如圖1所示．采用譜元法，可方便地考慮地層剪切波速沿不同方向的變化．考慮土層的自重及密實度的影響，土體的密度及剪切模量隨著埋深的增加逐漸增大，當(dāng)垂直深度大于3,000,m時，假定為均質(zhì)基巖半空間．

圖1 計算模型

圖2 震源激勵函數(shù)曲線

2 計算方法

譜元法的基本思想是選取以正交多項式表示的基函數(shù)，基函數(shù)與Fourier級數(shù)中的正弦和余弦相類似，可以提高級數(shù)表示的解的收斂速度，其基本步驟是：①將計算區(qū)域分成許多子域(單元)；②在每個子域中將近似解表示成截斷的Chebyshev或Lagrange正交多項式展開；③用Galerkin方法求解正交問題的變分格式，得到全局的近似解．譜元法是彈性波方程在空間近似的高階變分方法，該方法能減少用于離散物理空間域的單元的數(shù)量，從而在很大程度上降低了計算成本[16]．

研究地震波在土體的傳播，應(yīng)該以半無限空間體為研究對象，但往往要取適當(dāng)?shù)挠邢迏^(qū)域來計算，所以邊界問題就需要考慮設(shè)置人工吸收邊界[17]．以二維有限空間域?∈R2為例，計算解決彈性波在不均一介質(zhì)中傳播問題．空間域?的邊界為Γ，包括內(nèi)部邊界Γint和人工外部邊界Γext．位移和速度分別用u(x,t)和v(x,t)來表示，包括計算域和邊界，t∈I=[0,T ]，T為地震波的傳播時間．

2.1 彈性波動譜元模擬的基本理論[18-19]

彈性波傳播的運動方程為

式中：ρ=ρ(x)為介質(zhì)密度；f(x,t)為介質(zhì)體力；u0(x)、v0(x)分別為初始位移和速度；σ(x,t)為應(yīng)力張量；T(x,t)為Neumann邊界條件中的外力；g(x,t)為Dirichlet邊界條件中的位移；分別表示速度和位移對時間的偏微分；div[σ]=σij,j，σ·n=σijnj．

外部邊界條件問題也早有研究．地球物理學(xué)在20世紀(jì)80年代由Bayliss和Turkel開始將數(shù)值近似算法應(yīng)用于邊界問題；20世紀(jì)90年代Givoli和Keller采用遠(yuǎn)場漸進展開，推導(dǎo)出非局部吸收邊界條件已經(jīng)被廣泛采用．在此采用基于Engquist和Majda (1977)、Clayton和Engquist(1977)推導(dǎo)的旁軸條件的變分形式，在沿著人工邊界部分采用彈性動力學(xué)方程有限波數(shù)擴張的方法，近似表示人工邊界的瞬態(tài)阻抗．這種方法由Lysmet和Kuhlemeyer(1969)首次提出，對于高頻波和較小入射角的波的反射有很高的精度．在人工邊界extΓ，邊界條件表示為

式中：t為邊界上的力；n為邊界指向外的法向量；vT=v-[v·n]表示速度在邊界上的投影；cp和cs分別為P波和S波的傳播速度．

2.2 控制方程的變分形式——弱形式

用數(shù)值的方法解決波動方程，可以采用弱積分表達式，使原始問題獲得更廣的數(shù)值解．首先引入測試函數(shù)空間?

在方程(1)的兩邊同時乘以測試函數(shù)w(x)，同時聯(lián)合邊界條件，用分部積分法積分可得

這里有初始條件

式中：a(.,.)表示初始應(yīng)力虛功的雙線性形式，表達式為

式中：c表示方程(9)中的彈性模量．

2.3 數(shù)值離散

像有限元中一樣，譜元中也把計算區(qū)域離散為互不重疊的四邊形單元：每個單元的測試函數(shù)表示為每個單元都映射到一個正方形參考域s，定義在2,[1,1] ΛΛ=-局部坐標(biāo)系(ξ ,)η中．同時生成可逆的單元映射矩陣，坐標(biāo)映射x=?e(ξ ,η)．

對于數(shù)值離散來說，應(yīng)該保證變分形式滿足分段多項式的近似空間，其中h表示基本單元網(wǎng)格的長度

式中：L2(?)2為定義在?上的平方可積函數(shù)空間；[PN(s)]2為所有N階多項式的張量空間．譜單元的離散空間可以用離散點(nel,N)表示，每個全局坐標(biāo)域中的積分都使用局部坐標(biāo)的張量積來表示．同時為了便于因式分解，(N+1)2個基點作為每個單元的積分點，這些積分點將單元劃分為配置網(wǎng)格，設(shè)為，也就是(N+1)個Gauss-Lobatto-Legendre積分點的二維張量積．

通過以上概述，它可以給出對于波動微分方程

的耦合方法，其中ut和vt表示所有節(jié)點在t時刻的位移和速度向量，F(xiàn)ext表示外力矩陣，F(xiàn)int表示內(nèi)部節(jié)點力矩陣．譜方法由于Lagrange插值和GLL積分采用相同的配置點，所以得到對角質(zhì)量矩陣M，大大降低了計算難度．

譜元法整合了有限元的幾何模型的靈活性和譜技術(shù)的快速收斂性．由于譜元法的穩(wěn)定性、GLL積分和插值的高精確度，通過設(shè)置Legendre多項式的階數(shù)和劃分單元數(shù)，可以使誤差達到很高精度．本文的模型劃分單元網(wǎng)格的長度為100,m，網(wǎng)格每條邊上有8個插值點．

該方法通過參數(shù)α和β，迭代出離散動量方程在時間tn+a的位移速度和加速度

3 模型計算結(jié)果與分析

假設(shè)出平面線源埋深設(shè)置在地表下5,km，考慮到邊界的影響，計算范圍取5,R(R為盆地半徑)，如圖3所示．從不同發(fā)震主頻率、不同剪切線源位置和不同的盆地形狀來研究盆地地震響應(yīng)，觀察參數(shù)變化對地表地震動強度的影響．為便于對比，圖4給出了無沉積盆地時，半空間地表在不同主頻率波源作用下x＝0點處的位移時程曲線．假設(shè)半空間介質(zhì)波速線性變化，波源在x＝0點正下方5,km處．圖5給出了主頻為1,Hz的出平面線源作用在不存在沉積(圖5a)、存在沉積(R＝1,000,m，圖5b)兩種地形時，地表不同接受點的位移時程曲線對比．

圖3 半圓形沉積盆地計算模型(R＝1,000,m)

圖4 無沉積盆地時x＝0點的位移時程曲線

圖5 地表不同接受點的位移時程曲線對比

3.1 不同主頻率的點源的影響

以從沉積的正下方5,km入射Ricker子波為點源，主頻率依次為0.5，1，2，5,Hz，其他參數(shù)相同．以R＝1,000,m的沉積盆地為研究對象，由于對稱分布，選取地表橫坐標(biāo)為0，500，900，2,000,m為記錄點，記錄各點位移時程，結(jié)果如圖6所示．

由圖6對比可知：①沉積盆地對地震動有很強的放大效應(yīng)，位移幅值最大達到?jīng)]有沉積的2.3倍；②沉積場地的卓越頻率在0.5～1.0,Hz，當(dāng)頻率＞1時，地表的位移幅值則隨著頻率的增加衰減很快，低頻位移幅值達到高頻位移幅值的2～4倍；③盆地的邊緣效應(yīng)比較明顯，低頻和高頻位移幅值之比在邊緣部分最大達到4倍；④沉積的存在對于沉積外的影響甚微，可以忽略．

圖6 地表不同接受點的位移時程曲線

3.2 不同點源位置對地表響應(yīng)的影響

取主頻為1,Hz的點源作用在埋深為5.0,km、水平坐標(biāo)分別為0，500，800，1,000，1,500,m、沉積盆地半徑為1,000,m的半圓形沉積上，分別做地表－2,000，－1,000，－900，0，900，2,000,m等處接受點位移時程曲線，如圖7所示．

由圖7可知，點源的水平位置的偏移對于地表各點的地震動影響較小，可能由于水平位移量相對于垂直埋深變化較小；但是由于沉積對地震波的阻抗，在x =－1,000,m處衰減最大，達到20%；x =－900,m處衰減了15%，其余點的衰減在10%以內(nèi)．

圖7 不同震源偏移量時地表不同接受點的位移時程曲線

3.3 不同形狀的盆地對其地震響應(yīng)的影響

通過設(shè)置不同形狀的盆地，模擬沉積地形對地震效應(yīng)的影響．主要研究了如下參數(shù)的地形：① 半徑為1,000,m的半圓形盆地；② a＝1,000,m，b＝500,m的淺橢圓形盆地；③ a＝1,000,m，b＝1,500,m的深橢圓形盆地，其示意見圖8．沉積部分的波速在300～800,m/s間均勻變化；外部在深度3,km內(nèi)，沿深度方向800～1,500,m/s均勻變化，震源采用1,Hz的主頻出平面點源波源．不同地形在地表的位移曲線如圖9所示．

圖8 地形示意圖

圖9 地表記錄站點在不同地形條件下合成的位移地震圖

由圖9模擬結(jié)果可知：①在地震波第一次到達中心點(x＝0)時，深橢圓盆地的位移振幅偏大，比淺橢圓盆地大20%；在第一次波谷到達中心之后，地震波能量由兩邊繼續(xù)向中心傳播，淺橢圓盆地能量消散相對較少，會再引起比較大的位移，而深橢圓盆地由于軟土的阻抗作用，位移振幅相對較??；②盆地的邊緣附近，淺橢圓盆地的響應(yīng)要大，而且邊緣內(nèi)外地震影響差異很明顯；內(nèi)部有較大位移，外部相對平穩(wěn)．這是由于邊界的反射作用，地震波在邊緣位置來回多次疊加，形成很明顯的邊緣效應(yīng)．

4 結(jié) 論

(1)沉積盆地對地震動有很強的放大效應(yīng)．

(2)地震波的頻率和介質(zhì)的頻率相近時，地震動的響應(yīng)強烈．

(3)盆地的邊緣效應(yīng)比較明顯，在沉積盆地的外側(cè)，基本不受沉積盆地的影響．

(4)對于半圓形沉積盆地來說，點源的相對位置的變化對盆地邊緣兩側(cè)附近的影響較大，對其余點的影響很?。?/p>

(5)不同形狀沉積盆地在出平面線源作用下，深橢圓盆地的地震響應(yīng)較大，能量消散相對較快，由于介質(zhì)邊界對波的反射和折射，圖像表現(xiàn)出明顯的邊緣效應(yīng)．

(6)選取恰當(dāng)?shù)膱龅貙ㄖ目拐鹩泻艽笞饔?，對于在沉積土層上的建筑物，比如天津的沖積土，應(yīng)加強抗震設(shè)防，同一建筑物應(yīng)避免建在不同的地質(zhì)條件上．

［1］ 李雪強. 沉積盆地地震效應(yīng)研究[D]. 哈爾濱：中國地震局工程力學(xué)研究所，2011：1-88.

［2］ 杜啟振，楊慧珠. 方位各向異性粘彈性介質(zhì)波場有限元模擬[J]. 地球物理學(xué)報，2003，52(8)：2,010-2,014.

［3］ KHAM M，SEMBLAT J F，BOUDEN-ROMDHANE N. Amplification of seismic ground motion in the tunis basin：numerical BEM simulations vs experimental evidences[J]. Engineering Geology，2013，155：80-86.

［4］ 王秀明，張海瀾. 用于具有不規(guī)則起伏自由表面的介質(zhì)中彈性波模擬的有限差分算法[J]. 中國科學(xué)，2004，34(5)：481-493.

［5］ 楊 瑩. 二維地震波場有限差分法數(shù)值模擬研究[D].北京：中國地質(zhì)大學(xué)，2009：1-48．

［6］ 劉有山，滕吉文，徐 濤，等. 三角網(wǎng)格譜元法地震波場數(shù)值模擬[J]. 地球物理學(xué)進展，2014，29(4)：1,715-1,726.

［7］ 劉啟方，于彥彥，章旭斌. 施甸盆地三維地震動研究[J]. 地震工程與工程振動，2013，33(4)：54-60.

［8］ 李洪建，韓立國，鞏向博. 復(fù)雜構(gòu)造網(wǎng)格化及高精度地震波場譜元法數(shù)值模擬[J]. 石油物探，2014，53 (4)：375-383.

［9］ TENG Y C. Three-dimensional finite element analysis of waves in an acoustic media with inclusion[J]. J Acoust Soc,Am，1988，86：414-422.

［10］ MARFURT K J．Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equation[J]. Geophysics，1984，49：533-549.

［11］ MADARIAGA R．Dynamics of an expanding circular fault[J]. Bull Seism Soc Am，1976，65：163-182.

［12］ LEVANDER A R. Fourth-order finite-difference P-SV seismograms[J]. Geophysics，1988，53：1,452-1,463．

［13］ ORSZAG S A．Spectral methods for problems in complex geometries[J]. J Comput Phys，1980，37：70-92.

［14］ PATERA A T. A spectral element method for fluid dynamics：laminar flow in a channel expansion[J]. J Comput Acoust，1994，2(4)：371-422.

［15］ 陳敬國. 波場模擬中的震源—Ricker子波淺析[EB/OL]. 北京：中國科技論文在線，2006[2014-09-10]. http：//www. paper. edu. cn/html/releasepaper/ 2006/07/105/

［16］ 王童奎，李瑞華，李小凡，等. 譜元法數(shù)值模擬地震波傳播[J]. 防災(zāi)減災(zāi)工程學(xué)報，2007，27(4)：470-477.

［17］ DIMITRI K，JEAN-PIERRE V．The spectral element method：an efficient tool to simulate the seismic response of 2D and 3D geological structures[J]. Bulletin of the Seismological of America，1998，88(2)：368-392.

［18］ DIMITRI K，JEAN-PIERRE V．The spectral element method for elastic wave equations application to 2-D and 3-D seisnic problems[J]. International Journal for Numerical Methods in Engineering，1999(45)：1,139-1,164.

［19］ DAUKSHER W，EMERY A F. The solution of elastostatic and elastodynamic problems with Chebyshev spectral fnite elements[J]. Comput Methods Appl Mech Engrg，2000(188)：217-233.

The Research of Seismic Response of the Sedimentary Basin Under the Anti-plane Shear Line Source Based on the Simulation of the Spectral Element Method

HAN Jian-bin，LIU Zhong-xian，LUO Zhao-hui
(School of Civil Engineering，Tianjin Chengjian University，Tianjin 300384，China)

Due to the difference of the medium impedance between the internal and external of basin,and the special geometric shapes of basin,there is a remarkable ground motion amplification effect in the basin. In this thesis,the spectral element method has been used to study the effect of the seismic response of the sedimentary basin caused by different main frequency of anti-plane shear line sources,different point source location,and different sedimentary shapes. The result shows that the frequencies of point source and sedimentary shape have greater influences on the surface response. A different position of the point source for the edge effect has a certain influence.

sedimentary basin；anti-plane shear line source；spectral element；seismic response

TU435；P315.31

A

2095-719X(2015)02-0077-07

2014-10-14；

2014-12-01

國家自然科學(xué)基金(51278327)

韓建彬（1988—），男，河北石家莊人，天津城建大學(xué)碩士生.