肖 遙，丁 浩，黃 河

（遼寧省地震局，遼寧 沈陽 110034）

?

基于有限元方法研究岫巖羅圈里隕石坑地形對(duì)地震反應(yīng)的影響

肖 遙，丁 浩，黃 河

（遼寧省地震局，遼寧 沈陽 110034）

摘要：基于有限元方法對(duì)岫巖羅圈里隕石坑地形對(duì)地震反應(yīng)的影響進(jìn)行了探索。結(jié)合文獻(xiàn)資料、巖土勘察工作建立了岫巖羅圈里隕石坑二維有限元分析模型。選用粘彈性人工邊界和等效節(jié)點(diǎn)力輸入方法解決了有限元模擬半空間無限域時(shí)的邊界問題和地震動(dòng)輸入問題。為了突出地形因素的影響，采用彈性材料介質(zhì)進(jìn)行數(shù)值模擬分析，并與水平場(chǎng)地模型的分析結(jié)果進(jìn)行了對(duì)比。數(shù)值模擬結(jié)果和對(duì)比表明：羅圈里隕石坑地形對(duì)于隕石坑中部的地震動(dòng)有顯著的放大作用，遠(yuǎn)高于水平場(chǎng)地模型。此外，盆地中部短周期部分的反應(yīng)譜也明顯高于盆地邊緣和水平一維場(chǎng)地模型，這是導(dǎo)致海城地震中盆地中部自振周期較短的單層房屋破壞明顯嚴(yán)重的主要原因。

關(guān)鍵詞：隕石坑地形；地震反應(yīng)；粘彈性人工邊界；等效節(jié)點(diǎn)力輸入

0 引言

羅圈里位于遼寧省鞍山市岫巖滿族自治縣，是一個(gè)由地外天體撞擊形成的隕石坑。該坑形態(tài)呈碗狀，直徑約1.8公里，坑底與坑唇山脊高差為135～230 米，平均150米，地貌形態(tài)上表現(xiàn)為一個(gè)近圓形盆地。隕石撞擊坑形成后，坑內(nèi)曾積水形成了小湖泊并沉積了上百米厚的湖泊沉積物［1］。約在3.9萬年前，地勢(shì)較低的東部形成了缺口，積水流出并侵蝕山體，最后湖泊消失，演變?yōu)槿缃袼吹降耐胄伟嫉兀▓D1）。

圖1　岫巖羅圈里隕石坑全景照片F(xiàn)ig.1 Panoramic photograph of crater in Xiuyan Luoquanli

在海城地震中，羅圈里距1975年海城7.3級(jí)地震震中約65km，位于VI度區(qū)范圍內(nèi)，但顯示為高烈度異常點(diǎn)，烈度高達(dá)VIII度。震害嚴(yán)重的地方集中于盆地中心羅圈里小隊(duì)（當(dāng)時(shí)稱呼）。該隊(duì)共有房屋120多間，其中倒塌39間，嚴(yán)重破壞23間，其余輕微破壞，并有一人在地震中遇難。盆地以外的附近地區(qū)，震害都很輕微，基本上不超過VI度。地形條件是造成羅圈里地震反應(yīng)異常的主要原因之一，羅圈里是由隕石撞擊形成的微型盆地。很多研究表明盆地地形對(duì)于地震動(dòng)具有明顯的放大作用，但對(duì)于隕石坑這類微型盆地的放大作用尚沒有基于定量數(shù)值模擬分析的研究。本項(xiàng)工作將基于有限元方法對(duì)羅圈里隕石坑地形對(duì)地震反應(yīng)的影響進(jìn)行研究。

1 有限元分析中的幾個(gè)技術(shù)問題

1.1 人工邊界問題

在分析羅圈里隕石坑場(chǎng)地地形對(duì)地震反應(yīng)的作用時(shí)，考慮波動(dòng)能量向無限域地基的逸散影響是重要的。目前解決半無限域空間問題最常用的有限元數(shù)值方法是在截取的有限域上設(shè)置人工邊界，對(duì)分析場(chǎng)地地震作用下反應(yīng)設(shè)置合理的人工邊界對(duì)正確反映模型的整體動(dòng)力特性非常重要。

圖2　半無限空間系統(tǒng)分析模型Fig.2 Analysis model of semi-infinite space system

圖3　均勻半空間粘彈性人工邊界模型及加載示意圖Fig.3 A viscoelastic artificial boundary model of uniform half-space with loading sketch

圖4　遠(yuǎn)置人工邊界數(shù)值模擬結(jié)果Fig.4 Numerical simulation results of the far-set artificial boundary

圖4中給出了遠(yuǎn)置人工邊界的數(shù)值模擬結(jié)果，可以看出，應(yīng)力波隨著時(shí)間向遠(yuǎn)處傳播，完全投射出分析區(qū)域，用遠(yuǎn)置人工邊界方法模擬半無限空間問題可以得到相對(duì)精確的結(jié)果。但對(duì)于分析區(qū)域范圍來講遠(yuǎn)置人工邊界需要更大的分析空間，帶來計(jì)算成本的上升，對(duì)于算例中簡(jiǎn)單的二維問題，與使用粘彈性人工邊界相比，使用遠(yuǎn)置人工邊界劃分的網(wǎng)格數(shù)多出了100倍，計(jì)算時(shí)間多出了十幾倍。對(duì)于復(fù)雜問題和三維問題遠(yuǎn)置人工邊界付出的計(jì)算代價(jià)是不可承受的。

最簡(jiǎn)單的局部人工邊界是Lysmer與Kuhlemeyer提出的粘性人工邊界，該邊界通過沿人工邊界設(shè)置一系列阻尼器來吸收射向人工邊界的波動(dòng)能量，從而達(dá)到模擬波射出人工邊界的透射過程。表達(dá)式為：

其中，c為阻尼系數(shù)；ρ為介質(zhì)密度；cb為介質(zhì)波速。對(duì)一維問題，粘性邊界為精確邊界，對(duì)于多維問題為近似解，不能吸收大角度入射的外傳波能量。雖然，粘性邊界精度較低，僅為一階計(jì)算精度，但沿人工邊界設(shè)置的阻尼系數(shù)與頻率無關(guān)，處理方法簡(jiǎn)單，物理概念清晰，在工程中獲得較多應(yīng)用。

粘彈性邊界是基于衰減散射波表達(dá)，采用與粘性邊界相類似的推導(dǎo)過程建立起來的一種局部人工邊界條件。它通過沿人工邊界設(shè)置一系列由線性彈簧和阻尼器組成的簡(jiǎn)單物理元件來吸收射向人工邊界的波動(dòng)能量和反射波的散射，從而達(dá)到模擬波射出人工邊界的透射過程。在任意距點(diǎn)源半徑為rb的人工邊界點(diǎn)j上，邊界應(yīng)力與速度和位移的關(guān)系為：

等價(jià)于一阻尼系數(shù)ρcj的阻尼器并聯(lián)上一個(gè)剛度系數(shù)為的線性彈簧。若取=0，則粘彈性邊界退化為粘性邊界。

二維粘彈性人工邊界等效物理系統(tǒng)的彈簧系數(shù)KB和阻尼系數(shù)CB分別為：

式中：KBN、KBT分別為彈簧法向與切向剛度；R為波源至人工邊界點(diǎn)的距離；cs和cp分別為S波和P波波速；G為介質(zhì)剪切模量；ρ為介質(zhì)質(zhì)量密度；A為人工邊界代表的面積（二維問題中為長(zhǎng)度）；Nα與Tα分別為法向與切向粘彈性人工邊界參數(shù)。由于粘彈性人工邊界是基于全空間波動(dòng)理論推導(dǎo)得出的，直接用于半空間問題時(shí)，粘彈性動(dòng)力人工邊界的剛度系數(shù)會(huì)偏大，計(jì)算表明，粘彈性人工邊界有良好的魯棒性，人工邊界參數(shù)Nα與Tα可以在一定范圍內(nèi)取值，均可以給出良好的計(jì)算結(jié)果，在二維問題中推薦Nα的取值范圍［2］為［0.8，1.2］，Tα的取值范圍為［0.35，0.65］。

圖5展示了有限元遠(yuǎn)置人工邊界擴(kuò)展解、固定邊界、采用粘性人工邊界、粘彈性人工邊界的計(jì)算結(jié)果。從圖5給出的結(jié)果可以清楚地看出：如不采用人工邊界，地震動(dòng)將在邊界處產(chǎn)生明顯反射；粘彈性邊界與有限元遠(yuǎn)置人工邊界擴(kuò)展解結(jié)果更加相近，尤其是對(duì)于A、B、D三點(diǎn)具有較好的計(jì)算精度。因此本文中將采用粘彈性人工邊界來分析羅圈里隕石坑地形對(duì)于地震反應(yīng)的影響。

圖5　均勻半空間不同觀測(cè)點(diǎn)位移時(shí)程Fig.5 Displacement time history of different observation points in uniform half- space

1.2 地震動(dòng)輸入問題

地震動(dòng)輸入的具體實(shí)現(xiàn)方法與所采用的人工邊界條件密切相關(guān)，設(shè)置粘彈性邊界后，加速度形式的地震動(dòng)輸入方法不適用于外源波動(dòng)問題［3］，需要將人工邊界入射場(chǎng)轉(zhuǎn)化為等效邊界力來解決粘彈性人工邊界外源輸入問題，并在二維問題中應(yīng)用，這里稱該方法為等效節(jié)點(diǎn)力法。

以平面問題S波從人工邊界底部垂直入射為例，等效荷載計(jì)算公式為：

底邊界：

側(cè)邊界：

式（3）～（5）中：?1t= l cS，?2t=（ 2 H - l） cS， ρ、cS、 λ分別為介質(zhì)密度、S波波速、拉梅常數(shù)；H為底邊界到地表的距離，l為B結(jié)點(diǎn)到底邊界的距離；1t?、?2t 分別為B結(jié)點(diǎn)處底邊界入射S波和地表反射S波的時(shí)間延遲；AB為人工邊界面上邊界結(jié)點(diǎn)B的影響面積；等效地震荷載的下標(biāo)代表結(jié)點(diǎn)號(hào)和分量方向，上標(biāo)代表結(jié)點(diǎn)所在人工邊界面的外法線方向，與坐標(biāo)軸方向一致為正，相反為負(fù)。

根據(jù)上述計(jì)算等效地震荷載的方法編制了計(jì)算程序ABEF，用于計(jì)算邊界結(jié)點(diǎn)的等效荷載，從而實(shí)現(xiàn)粘彈性邊界單元的波動(dòng)輸入。下面通過一個(gè)數(shù)值算例對(duì)以上方法進(jìn)行驗(yàn)證。

圖6　等效節(jié)點(diǎn)力驗(yàn)證模型Fig.6 Equivalent node force verification model

如圖6所示的二維彈性半空間，半空間介質(zhì)密度ρ=1900kg/m3，剪切波速VS=300m/s，根據(jù)經(jīng)典的彈性理論，剪切模量G=ρ× VS2，彈性模量E=2（1 +μ）G ，泊松比μ=0.35。從底部垂直入射S波，入射波加速度、速度和位移時(shí)程見圖7。如圖8所示：可以看出，地表反應(yīng)最大值接近入射波幅值的2倍，而且向地表方向加速度幅值逐步增加，這些結(jié)論均與理論解相吻合，因此可以說明本文采用的粘彈性邊界單元及波動(dòng)輸入方法是可以用于數(shù)值模擬的。

2 羅圈里地震地質(zhì)環(huán)境和土層結(jié)構(gòu)介紹

陳鳴等［1］在岫巖羅圈里隕石坑中心位置實(shí)施了科學(xué)鉆探，鉆孔深度為307m，鉆孔揭露了坑內(nèi)上部107m厚的湖泊相沉積物，以及下部厚度為188m的角礫巖堆積透鏡體，坑內(nèi)填充的角礫巖是不同程度沖擊變質(zhì)巖石的混合堆積物（圖9）。為了進(jìn)一步研究岫巖羅圈里隕石坑的盆地效應(yīng)，查明岫巖羅圈里隕石坑第四系土層結(jié)構(gòu)，黃河等①黃河，肖遙等，海城地震程度異常研究，沈陽：遼寧省地震局，2016在岫巖羅圈里開展了鉆探工作。鉆探采用了XY-150型鉆機(jī)，共完成鉆孔4個(gè)，深度50.0至100.0m，總進(jìn)尺250m，每個(gè)孔均進(jìn)行了波速測(cè)試，測(cè)試間距為2.0m。根據(jù)鉆探結(jié)果揭示，場(chǎng)地主要有雜填土，含礫粉質(zhì)粘土，有機(jī)質(zhì)土，粘性土等組成，其中鉆孔1（100m）的土層結(jié)構(gòu)和波速列于表1。

樊華等（2013）利用中國(guó)地震局地球物理勘探中心在岫巖羅圈里隕石坑取得的淺層地震折射資料，采用地震走時(shí)層析成像法進(jìn)行反演，得到兩條近垂直交叉剖面的P波速度結(jié)構(gòu)［4］。根據(jù)以上資料，可以獲得羅圈里隕石坑的剖面結(jié)構(gòu)（圖11），通過該結(jié)構(gòu)可以建立羅圈里隕石坑的有限元分析模型。

圖7　輸入地震動(dòng)時(shí)程Fig.7 Ground motion time history

圖8　各觀測(cè)點(diǎn)加速度時(shí)程Fig.8 Acceleration time history of each observation point

圖9　岫巖羅圈里隕石坑深部科學(xué)鉆位置及柱狀圖示意Fig.9 Deep scientific borehole position and column of Xiuyan Luoquanli crater

3 有限元計(jì)算和對(duì)比

3.1 盆地模型概況和基巖輸入地震動(dòng)

圖10　L1 和 L2 測(cè)線位置以及速度結(jié)構(gòu)圖Fig.10 Line position and velocity structure diagram of L1 and L2

表1　羅圈里鉆孔1（100m）橫波波速

圖11　剖面結(jié)構(gòu)示意圖Fig.11 Sketch map of profile structure

圖12顯示了根據(jù)以上文獻(xiàn)資料和鉆探結(jié)果，在有限元軟件ABAQUS中建立了羅圈里盆地的數(shù)值分析模型，在兩側(cè)和邊界添加粘彈性人工邊界單元以模擬無限域?qū)Φ卣饎?dòng)波動(dòng)的影響，加載時(shí)采用施加等效結(jié)點(diǎn)力的方法。為盡可能的減少計(jì)算時(shí)間并保證一定的計(jì)算精度有限元單元網(wǎng)格劃分程度L＜Vs/10。以剪切波速大于500m/s的沖擊角礫巖作為地震動(dòng)輸入界面。其中湖泊相沉積層密度ρ=1850kg/m3，彈性模量E=6.0E9Pa，泊松比μ=0.35；沖擊角礫巖密度ρ=2400kg/m3，彈性模量E=1.5E10Pa，泊松比μ=0.25。基于羅圈里盆地的鉆孔資料建立了觀測(cè)點(diǎn)A的一維波動(dòng)模型，并利用專業(yè)土層反應(yīng)分析軟件ESE對(duì)場(chǎng)地地震反應(yīng)進(jìn)行了一維頻域分析。

地震波從遠(yuǎn)處的震源經(jīng)巖土介質(zhì)傳播到地表，需經(jīng)歷多次的反射與折射，根據(jù)折射定律，地震波傳播到地表時(shí)將近似于豎直方向。為此，本文將基巖輸入地震動(dòng)視為豎直向上傳播的 SH 波。按照震中距相近的原則選擇了San Fernando地震中的Cedar Springs臺(tái)站記錄到的一條天然地震波，獲得了19.2gal加速度峰值水平的地震波作為地震動(dòng)輸入，如圖13所示。

圖12　羅圈里隕石坑數(shù)值分析模型示意圖Fig.12 Numerical analysis model sketch map of Luoquanli crater

圖13　輸入地震動(dòng)加速度、速度、位移時(shí)程圖Fig.13 Acceleration， velocity and displacement time history diagram of input ground motion

3.2 計(jì)算結(jié)果與分析

在彈性介質(zhì)假設(shè)的情況下，突出了場(chǎng)地地形因素即盆地效應(yīng)的影響，從而避免了土介質(zhì)非線性反應(yīng)對(duì)于地形因素即盆地效應(yīng)的干擾。因此盡管彈性介質(zhì)假設(shè)的模擬結(jié)果會(huì)與實(shí)際地震反應(yīng)分析情況有較大的偏差（這是因?yàn)榉浅＼浀?07m的湖泊相沉積將對(duì)地震動(dòng)產(chǎn)生非常重要的影響），但對(duì)于定性分析地形因素的影響仍具有較強(qiáng)的意義。從圖14可以看出，盆地地表不同位置對(duì)于地震反應(yīng)的放大作用有顯著的不同，最大的放大作用出現(xiàn)在觀測(cè)點(diǎn)B，觀測(cè)點(diǎn)A的放大作用也較大，而觀測(cè)點(diǎn)C、D的放大作用相對(duì)較小，觀測(cè)點(diǎn)D的放大作用最小，這與實(shí)際海城地震中羅圈里的震害情況相符合。并且本文對(duì)于土層厚度100m的水平場(chǎng)地模型（對(duì)應(yīng)于隕石坑中部）進(jìn)行了計(jì)算，結(jié)果見圖15。對(duì)比圖14和圖15可以看出二維隕石坑模型對(duì)于隕石坑中部地表加速度峰值的放大作用明顯更強(qiáng)，達(dá)到水平場(chǎng)地模型的1.5倍左右。這是由于地震波在盆地特殊地形的作用下反射疊加，形成“聚焦效應(yīng)”使得地震作用明顯增強(qiáng)。另外，也可以看出二維隕石坑模型中地震波在15秒至30秒之間仍有較強(qiáng)的作用，而水平一維場(chǎng)地模型中15秒之后地震波很快衰減，即隕石坑地形會(huì)明顯加強(qiáng)地震動(dòng)的持時(shí)，這也與在多次震害實(shí)際中地震波在盆地地形中反射會(huì)延長(zhǎng)地震波的作用時(shí)間這一特征相符合。

圖16對(duì)比了各觀測(cè)點(diǎn)以及水平地形的結(jié)果的反應(yīng)譜。從圖16中可以看出在短周期部分（0.1s至0.2s）地震反應(yīng)譜大小順序基本為A、B、C、D，該周期段也恰好為單層磚房結(jié)構(gòu)的自振周期段，單層磚房結(jié)構(gòu)的破壞主要受自振周期T1點(diǎn)對(duì)應(yīng)的反應(yīng)譜加速度大小Sa （T1）控制，這一模擬結(jié)果與海城地震中反應(yīng)出來的實(shí)際震害情況相符合。

圖14　各觀測(cè)點(diǎn)地表加速度時(shí)程Fig.14 Acceleration time history of each observation point

圖15　水平場(chǎng)地地表加速度時(shí)程Fig.15 Acceleration time history of horizontal ground

圖16　加速度反應(yīng)譜對(duì)比分析Fig.16 Comparative analysis of acceleration response spectrum

4 結(jié)論

本研究基于有限元方法對(duì)岫巖羅圈里隕石坑地形對(duì)地震反應(yīng)的影響進(jìn)行了探索。結(jié)合文獻(xiàn)資料、巖土勘察工作建立了岫巖羅圈里隕石坑二維有限元分析模型。選用粘彈性人工邊界和等效節(jié)點(diǎn)力輸入方法解決了有限元模擬半空間無限域時(shí)的邊界問題和地震動(dòng)輸入問題。為了突出地形因素的影響，采用彈性材料介質(zhì)進(jìn)行數(shù)值模擬分析，并與水平場(chǎng)地模型進(jìn)行了對(duì)比。數(shù)值模擬結(jié)果和對(duì)比表明：羅圈里隕石坑地形對(duì)于隕石坑中部的地震動(dòng)有顯著的放大作用，遠(yuǎn)高于水平一維場(chǎng)地模型。此外，盆地中部短周期部分的加速度反應(yīng)譜也明顯高于盆地邊緣和水平一維場(chǎng)地模型，這是導(dǎo)致海城地震中盆地中部自振周期較短的單層房屋破壞明顯嚴(yán)重的主要原因。

參考文獻(xiàn)：

［1］陳鳴，肖萬生，等. 岫巖隕石撞擊坑的證實(shí)［J］.科學(xué)通報(bào)，2009，22（54）：3507-3511.

［2］劉晶波，谷音，等.一致粘彈性人工邊界及粘彈性邊界單元［J］.巖土工程學(xué)報(bào)， 28（9）：1070-1075.

［3］何建濤，馬懷發(fā)，等.黏彈性人工邊界地震動(dòng)輸入方法及實(shí)現(xiàn)［J］.水利學(xué)報(bào)，2010，41（8） ：960-968.

［4］樊華，袁一凡，等.地震走時(shí)層析成像方法及在隕石坑探測(cè)中的應(yīng)用［J］.CT理論與應(yīng)用研究，2013，22（4）：605-613.

FEM-BASED RESEARCH INTO THE INFLUENCE OF CRATERTERRAIN ON SEISMIC RESPONSE IN XIUYAN LUOQUANLI

XIAO Yao， DING Hao， HUANG He
（Earthquake Administration of Liaoning Province， Liaoning Shenyang 110034， China）

Abstract：The paper presents a research， which is based on the finite element method， into the influence of Crater-terrain on seismic response in Xiuyan Luoquanli. Xiuyan Luoquanli crater two-dimensional finite element analysis model was established by literature and geotechnical investigation work. The viscous-spring artificial boundary and the equivalent nodal force input method are used to solve the boundary problem and the ground motion input problem in finite element simulation. Elastic material medium and the comparation with the result of horizontal site model are used to focus on the influence of terrain factors. The results and comparison show that Luoquanli crater terrain has significant amplification effect for crater central vibration，and much higher than the level of horizontal site model. Moreover， the acceleration response spectrum values during short periods in centre are significantly higher than that the level of in the edge of the basin and the horizontal site， that is the primary reason for the serious damage of one-storey house with short natural period of vibration in centre of crater.

Key words：crater topography； seismic response； visco-elastic artificial boundary； equivalent nodal force input

中圖分類號(hào)：P315

文獻(xiàn)標(biāo)志碼：A

DOI:10.13693/j.cnki.cn21-1573.2016.02.001

文章編號(hào)：1674-8565（2016）02-0001-09

基金項(xiàng)目：遼寧省第二批科學(xué)技術(shù)計(jì)劃項(xiàng)目（項(xiàng)目編號(hào)2013231028）

收稿日期：2016-01-06

修訂日期：2016-02-28

作者簡(jiǎn)介：肖遙（1986-），男，黑龍江省賓縣人，2011年畢業(yè)于中國(guó)地震局工程力學(xué)研究所，碩士，工程師，現(xiàn)主要從事地震工程、地震風(fēng)險(xiǎn)評(píng)估、巖土地震反應(yīng)分析等方面的研究工作。