低山丘陵地形效应的三维有限元模拟

2023-12-15蒋其峰荣棉水许洪泰彭艳菊

地震工程学报 2023年6期

蒋其峰, 荣棉水, 魏玮, 许洪泰, 彭艳菊

(1. 山东省地震局, 山东济南 250014; 2. 北京工业大学城市建设学部, 北京 100124;3. 国家自然灾害防治研究院, 北京 100085)

0 引言

地形效应是地震学和地震工程学研究的重要内容。众多震害调查、地震观测研究表明,不规则地形对地震动的强度分布有非常重要的影响[1-3],山顶、陡坡等地形对地震动有较强的放大作用,这往往造成更大的破坏。2008年汶川地震时,单薄山脊、条形山体、坡型转折部位等地形放大效应明显,诱发或加重了地震破坏[4]。

数值模拟常被用于研究地形效应。Vincenzo DiFiore[5]曾建立二维局部场地模型使用有限元方法研究了斜坡对地震动的影响。李英民等[6]采用二维有限元数值分析方法,对黏弹性岩质坡地在白噪声输入下的响应进行计算,研究了坡地地形对竖向地震动反应谱特性的影响。杨宇等[7]以自贡西山公园为研究目标,建立二维有限元模型研究了山脊场地地形的影响。Mohammad Bararpour等[8]曾建立二维模型研究了山谷坡度等对地震动的影响。张季等[9]基于二维间接边界元法,研究了SH波入射时凸起场地的地形和土层放大效应。梁轩等[10]对不同高宽比的山体数值模型输入不同频率的Gabor小波,研究二维梯形山体地形效应及对周围场地地震反应的影响。

真实地形往往复杂多变,地形效应也呈现非常复杂的模式,与二维模拟相比,三维模拟更能够反映复杂地形对地震动的作用[11]。梁建文等[12]计算分析了三维凹陷地形地震响应,分析表明凹陷地形附近的地震动加速度反应谱与自由场的反应谱有明显的区别。蒋涵等[13]以庐山地区为研究对象,利用谱元法进行三维模拟,总结了山脊线与坡度和峰值速度放大系数的相关性规律。

2008年汶川地震时,自贡西山公园地形台阵记录到了众多的加速度时程,为研究地形效应提供了重要资料。本文基于自贡西山公园山体及附近地区的真实地形建立三维有限元模型,采用黏弹性人工边界,开展动力有限元计算,根据实际观测资料对比检验模拟结果,研究分析低山丘陵地形对地震动的影响。

在离震源较远时,地震波受到地层的影响,在接近地面时传播方向往往接近竖向。对于局部场地,地震波垂直入射所造成的影响一般比斜入射更大[14]。与压缩波相比,剪切波对地表建筑的危害更大,同时剪切波的影响也是建设工程震害防御工作中十分重要的问题。因此,本文地震波入射采用SH波垂直入射方式,运用ABAQUS有限元软件进行分析计算。

1 人工边界的建立

本文人工边界采用黏弹性人工边界,该边界较为成熟,精确度较高,在数值模拟中得到广泛应用。黏弹性人工边界需要在无限域的截断面上设置一系列的弹簧和阻尼元件(图1),用以模拟无限域中地震波的传播。

人工边界的弹簧和阻尼器参数取值如式(1)～(4)所示[15]。其中KNi为边界法向弹簧系数;KTi为边界切向弹簧系数;CNi为边界法向阻尼器系数;CTi为边界切向阻尼器系数;r取散射波源到人工边界面的距离;ρ为介质密度;λ、G分别为第一、第二拉梅常数;CP、CS分别为介质的纵波波速、横波波速。A、B的取值根据经验获得,分别取0.8和1.1。

(1)

(2)

CNi=BρCP

(3)

CTi=BρCS

(4)

地震动输入采用节点集中力输入的方式,与黏弹性人工边界相对应,采用杜修力等学者的剪切波垂直入射条件下地震载荷输入方法[16]。在这些理论的基础上,本文开发了黏弹性人工边界和集中力设置的有限元前处理程序。

2 人工边界设置验证

为检验人工边界设置方法和程序的准确性,本文建立了试验模型开展验证。试验模型为x轴、y轴、z轴方向长度分别为4 500 m、4 500 m、2 120 m的长方体,有限单元网格约90 m,考虑到三维地形变化复杂更适合用四面体网格进行划分,试验模型地表部分采用四面体单元,其余采用六面体单元(图2)。介质密度2 630 kg/m3,泊松比取0.22,弹性模量取32.5 GPa。输入剪切波位移时程采用主频1 Hz的雷克子波(图3),时长2 s,由xy面沿z轴方向垂直入射。

图2 试验模型及网格划分情况Fig.2 The experimental model and mesh division

图3 输入地震动位移时程与傅氏谱Fig.3 Displacement time history and Fourier spectrum of input ground motion

地表位移时程计算结果见图4。通过图4可知,地表位移时程有限元计算结果和理论值非常接近,说明所设置的黏弹性人工边界和地震载荷是合理的,将地表部分设置为四面体单元是可行的,这为研究更复杂的三维有限元地形模型奠定了基础。

图4 地表位移计算值与理论值的对比Fig.4 Comparison between calculated and theoretical values of surface displacement

3 自贡西山公园地震动模拟

自贡西山公园位于浅切割丘陵区,主要发育低山丘陵地貌,附近主要出露侏罗系基岩,第四系堆积物主要为残坡积物,厚度较薄,一般在2～4 m。自贡西山公园地形台阵主要由8个台站组成,记为S0～S7。S0台站位于山脚土层上,用于研究土层地震反应。S1～S7台站位于山体不同高程的基岩上面,可用于研究地形效应。其中S1台站位于山脚基岩,S6台站位于山顶基岩。每个台站都记录了EW、NS、UD三个方向的加速度时程。S1台站记录到的EW方向的加速度时程和傅氏谱见图5。根据传统谱比法,该山体的最大放大效应在0.1～20 Hz频带,水平向地震动大于竖向地震动;在0.1～10 Hz频带,EW方向地震动大于NS方向地震动[17]。由于EW方向地震动放大效应较为显著,所以本文使用EW方向振动的SH波作为输入。根据实际地震记录的傅氏谱,水平向地震动的频率主要集中在1 Hz左右的低频段。因此,数值模拟的输入波形采用低频地震动输入,共采用了三种主频的雷克子波波形,对应的加速度时程的主频依次为0.75 Hz、1 Hz、1.5 Hz(图6),三种时程的主要频率范围依次为[0.25 Hz,1.25 Hz]、[0.5 Hz,1.75 Hz]、[0.5 Hz,2.5 Hz]。输入加速度时程的最大值约为S1台站EW方向地震动峰值加速度的0.5倍,取为11.8 gal。

图5 S1台站EW方向加速度时程和傅氏谱Fig.5 The acceleration time history and its Fourier spectrum of S1 station in EW direction

图6 三种输入加速度时程的傅氏谱Fig.6 Fourier spectra of the three incident acceleration time histories

选取自贡西山公园及其附近地区作为研究区域,三维有限元模型根据研究区域的数字高程地形模型建立。数字高程模型采用ASTGTM数据,分辨率30 m。模型大小为2 010 m(长)×1 890 m(宽)×1 090 m(高)。考虑到自贡西山公园山体为基岩山体,且本文主要研究地形的影响,介质简化为均匀弹性介质。泊松比取0.22,横波波速取为800 m/s,密度取2 000 kg/m3[7]。网格大小取30 m左右,符合空间采样要求。模型及网格划分情况见图7(a),分布有台站的目标山体见图中白色方框。研究区域所采用的数字高程模型见图7(b),对应的坡度见图7(c)。

图7 自贡西山公园及附近地区三维有限元模型及高程、坡度分布Fig.7 The 3D finite element model of Xishan Park and nearby areas and the distribution of elevation and slope

4 结果检验与分析

为了检验结果的合理性,提取了山顶S6台站、山脚S1台站所在位置对应的节点的峰值加速度,计算了三种主频地震动输入下山顶S6台站峰值加速度与山脚S1台站峰值加速度之间的比值,本文称之为“放大系数”(表1)。从表1可知,输入加速度时程主频从0.75 Hz增大至1.5 Hz时,放大系数有随频率增加而增大的趋势,最高可达到1.2,总体上放大系数较小。

表1 不同输入地震动时程对应的峰值加速度放大系数Table 1 Amplification factor of PGA corresponding to different incident ground motions

地震动时程的峰值是多种频率地震动共同作用的结果,因此在保证地震动频带基本一致的情况下对比时程的峰值加速度较为合理。本文对S6台站、S1台站的东西向记录按照三种输入地震动的主要频带范围分别进行了带通滤波,滤掉了主要频带范围外的成分,并通过傅里叶逆变换得到了加速度时程,计算了滤波后的S6台站与S1台站的加速度时程峰值之间的比值,如表2所列。以[0.25 Hz,1.25 Hz]窗口滤波为例,加速度时程滤波前后对比见图8所示。对比表1和表2可知,模拟得到的峰值加速度的放大系数与相关频段的实际地震动时程峰值加速度放大系数基本一致,稍微偏小。偏小的原因可能是受地震波各频率成分强度占比、地层、构造等的影响。总体来说,本文的数值模拟是较合理的。

表2 按三种频带范围滤波后S6与S1加速度时程峰值比值Table 2 Peak value ratio of acceleration time history of S6 to S1 after filtering in three frequency bands

图8 加速度时程滤波前后对比(滤波窗口[0.25 Hz,1.25 Hz])Fig.8 Comparison between acceleration time histories before and after filtering (Window of [0.25 Hz,1.25 Hz])

为分析峰值加速度的强度及空间分布,本文对三维有限元计算结果进行了后处理,从地表各节点的加速度时程中提取了峰值加速度。将节点坐标与实际地理坐标进行了对应,绘制了不同主频的输入地震动作用下地表峰值加速度分布图(图9)。

图9 不同主频的输入加速度时程对应的地表PGA分布Fig.9 The PGA distributions corresponding to different input ground motions

对比三种不同频率地震波产生的地表峰值加速度分布情况可知,不同频率的地震波对山体的影响差别较大,地形效应与入射波频率密切相关。就研究区域内西北部的连续山体而言,在主频0.75 Hz的地震波影响下,地表峰值加速度普遍较高,但在主频1.0 Hz、1.5 Hz的地震波的影响下,该连续山体各山顶之间相对低洼位置的地表峰值加速度降低,这说明地震波频率对地表峰值加速度的分布有重要影响。

此外,三种主频入射波对应的计算结果显示,研究范围内的山体顶部存在不同程度的放大效应。随着输入加速度时程主频从0.75 Hz增加至1.5 Hz,较小规模山体的放大效应越来越明显。在波速相同时,频率增加意味着波长变短。在波长与山体空间尺度接近时,地震动更容易受到地形影响,这说明了较高主频的地震动可以凸显地形细节的影响。

低山丘陵的山体周围常分布有凹陷地形,通过对比分析可知,山脚凹陷位置地震动峰值加速度削弱明显。因为如果地面为平坦地形,那么理论上峰值加速度地表值为23.6 gal。而计算结果显示,许多凹陷地形的峰值加速度达不到这一水平(图9)。这一结果与梁建文等人对于三维凹陷地形的地震响应[12]的研究结论是一致的。

通过对地震动峰值加速度分布图的分析,可以看出,地表峰值加速度与地形有一定的联系,峰值加速度随高程的增加有增大的趋势,而坡度较陡的地方峰值加速度也较大。为进一步研究地形效应与地形要素的相关关系,提取了各节点对应的地表高程和坡度数据。为保证高程与所在山体位置的对应关系,本文以观测台站所在山体分析单个山体的高程、坡度与地表峰值加速度的对应关系,绘制了地表峰值加速度与高程、坡度之间的变化关系图(图10～12)。

根据地表峰值加速度与高程之间的关系图[图10(a)、图11(a)、图12(a)],地表峰值加速度有随高程增大而增大的趋势,两者有较好的正相关性,这与杨宇等人的研究成果一致[7]。这也意味着如果山体高度更高,放大系数可能会更大,比如Griffiths等[2]在对实际记录(由1 Hz地震仪记录)分析的基础上指出,高度236～394 m的山体顶部的放大系数达到1.7～3.4。在输入地震动主频较低时,两者相关性较强。随着输入地震动主频增大,高程中段值对应的点变得离散,两者相关性变弱。

图11 主频1.0 Hz的输入加速度时程对应的地表峰值加速度与高程、坡度关系Fig.11 Relationship between the peak acceleration of surface and the elevation and slope corresponding to the input acceleration time history with the dominant frequency of 1.0 Hz

图12 主频1.5 Hz的输入加速度时程对应的地表峰值加速度与高程、坡度关系Fig.12 Relationship between the peak acceleration of surface and the elevation and slope corresponding to the input acceleration time history with the dominant frequency of 1.5 Hz

根据地表峰值加速度与坡度之间的关系图[图10(b)、图11(b)、图12(b)],地表峰值加速度与坡度之间的相关性较弱,尤其是中低坡度值对应的峰值加速度分布较离散。经过分析,出现此类情况的主要原因是研究对象属于低山丘陵,山顶坡度比较平缓,坡度值较低,与山脚的坡度值接近[图7(c)],而两者对应的峰值加速度显然是不一样的,所以单靠坡度难以建立与峰值加速度之间的关系。这一点与高山峡谷地区的情况是不相同的。但两者并不是没有规律可言,从关系图中可知,坡度高值基本对应地表峰值加速度的中高值。而低山丘陵坡度高的地方大多位于山体中部位置[图7(b)、图7(c)],即坡度的高值一般对应高程的中段值,这启示我们如果要通过地形要素定量描述地形效应,那么综合高程和坡度建立地表峰值加速度预测模型比较合理。

为充分对比使用高程或坡度单个变量做峰值加速度的一元拟合与同时使用高程和坡度两个变量做峰值加速度的二元拟合的效果差异,本文以主频0.75 Hz的加速度时程输入下的峰值加速度作为实验对象,开展一元拟合与二元拟合的对比实验。拟合效果的优劣以拟合优度来描述,拟合优度是指回归公式对观测值的拟合程度,通常用确定系数R2来度量,R2最大值是1,越接近1,拟合效果越好。拟合方法采用多项式拟合,根据数据分布并经过测试筛选出合适的多项式次数。将高程与峰值加速度、坡度与峰值加速度、高度和坡度与峰值加速度的拟合依次记为拟合1、拟合2、拟合3,拟合采用的公式、拟合优度如表3所列,公式中系数的拟合结果如表4所列,拟合效果见图13。

表4 多项式系数拟合结果Table 4 Polynomial coefficients after fitting

图13 主频0.75 Hz的输入加速度时程对应的三种拟合结果Fig.13 Three fitting results corresponding to the input acceleration time history with the dominant frequency of 0.75 Hz

根据拟合结果可以得知,综合高程和坡度两个变量拟合峰值加速度,效果优于使用单个变量做拟合,二元三次多项式可以较好地拟合峰值加速度随高程和坡度的变化。如果定义某个地点为参考点,其他地点的峰值加速度除以参考点的峰值加速度即可得到放大系数,即本文放大系数是峰值加速度的同比例缩放。因此,二元三次多项式也可以较好地拟合地形放大系数随高程和坡度的变化。