基于有限元方法研究岫岩罗圈里陨石坑地形对地震反应的影响

肖 遥，丁 浩，黄 河

（辽宁省地震局，辽宁 沈阳 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年毕业于中国地震局工程力学研究所，硕士，工程师，现主要从事地震工程、地震风险评估、岩土地震反应分析等方面的研究工作。