基于有限元方法研究岫巖羅圈里隕石坑地形對地震反應的影響

肖 遙，丁 浩，黃 河

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

?

基于有限元方法研究岫巖羅圈里隕石坑地形對地震反應的影響

肖 遙，丁 浩，黃 河

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

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

關鍵詞：隕石坑地形；地震反應；粘彈性人工邊界；等效節點力輸入

0 引言

羅圈里位于遼寧省鞍山市岫巖滿族自治縣，是一個由地外天體撞擊形成的隕石坑。該坑形態呈碗狀，直徑約1.8公里，坑底與坑唇山脊高差為135～230 米，平均150米，地貌形態上表現為一個近圓形盆地。隕石撞擊坑形成后，坑內曾積水形成了小湖泊并沉積了上百米厚的湖泊沉積物［1］。約在3.9萬年前，地勢較低的東部形成了缺口，積水流出并侵蝕山體，最后湖泊消失，演變為如今所看到的碗形凹地（圖1）。

圖1　岫巖羅圈里隕石坑全景照片Fig.1 Panoramic photograph of crater in Xiuyan Luoquanli

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

1 有限元分析中的幾個技術問題

1.1 人工邊界問題

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

圖2　半無限空間系統分析模型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　遠置人工邊界數值模擬結果Fig.4 Numerical simulation results of the far-set artificial boundary

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

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

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

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

等價于一阻尼系數ρcj的阻尼器并聯上一個剛度系數為的線性彈簧。若取=0，則粘彈性邊界退化為粘性邊界。

二維粘彈性人工邊界等效物理系統的彈簧系數KB和阻尼系數CB分別為：

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

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

圖5　均勻半空間不同觀測點位移時程Fig.5 Displacement time history of different observation points in uniform half- space

1.2 地震動輸入問題

地震動輸入的具體實現方法與所采用的人工邊界條件密切相關，設置粘彈性邊界后，加速度形式的地震動輸入方法不適用于外源波動問題［3］，需要將人工邊界入射場轉化為等效邊界力來解決粘彈性人工邊界外源輸入問題，并在二維問題中應用，這里稱該方法為等效節點力法。

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

底邊界：

側邊界：

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

根據上述計算等效地震荷載的方法編制了計算程序ABEF，用于計算邊界結點的等效荷載，從而實現粘彈性邊界單元的波動輸入。下面通過一個數值算例對以上方法進行驗證。

圖6　等效節點力驗證模型Fig.6 Equivalent node force verification model

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

2 羅圈里地震地質環境和土層結構介紹

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

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

圖7　輸入地震動時程Fig.7 Ground motion time history

圖8　各觀測點加速度時程Fig.8 Acceleration time history of each observation point

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

3 有限元計算和對比

3.1 盆地模型概況和基巖輸入地震動

圖10　L1 和 L2 測線位置以及速度結構圖Fig.10 Line position and velocity structure diagram of L1 and L2

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

圖11　剖面結構示意圖Fig.11 Sketch map of profile structure

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

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

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

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

3.2 計算結果與分析

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

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

圖14　各觀測點地表加速度時程Fig.14 Acceleration time history of each observation point

圖15　水平場地地表加速度時程Fig.15 Acceleration time history of horizontal ground

圖16　加速度反應譜對比分析Fig.16 Comparative analysis of acceleration response spectrum

4 結論

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

參考文獻：

［1］陳鳴，肖萬生，等. 岫巖隕石撞擊坑的證實［J］.科學通報，2009，22（54）：3507-3511.

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

［3］何建濤，馬懷發，等.黏彈性人工邊界地震動輸入方法及實現［J］.水利學報，2010，41（8） ：960-968.

［4］樊華，袁一凡，等.地震走時層析成像方法及在隕石坑探測中的應用［J］.CT理論與應用研究，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

中圖分類號：P315

文獻標志碼：A

DOI:10.13693/j.cnki.cn21-1573.2016.02.001

文章編號：1674-8565（2016）02-0001-09

基金項目：遼寧省第二批科學技術計劃項目（項目編號2013231028）

收稿日期：2016-01-06

修訂日期：2016-02-28

作者簡介：肖遙（1986-），男，黑龍江省賓縣人，2011年畢業于中國地震局工程力學研究所，碩士，工程師，現主要從事地震工程、地震風險評估、巖土地震反應分析等方面的研究工作。