DEEPSOIL和FLAC3D一维土层地震反应分析对比研究

2022-07-06梁腾龙骆星宇周东覃钊

广西大学学报（自然科学版） 2022年2期

梁腾龙，骆星宇，周东*，覃钊

(1.广西大学土木建筑工程学院，广西南宁530004；2.广西防灾减灾与工程安全重点实验室，广西南宁530004)

0 引言

土层地震反应分析是工程结构抗震设计的重要环节。在场地地震安全性评价中，场地土层地震反应分析软件将基岩地震动时程和场地地质条件勘测转化为场地设计反应谱，为工程建筑物提供合理的抗震设计依据。地震反应分析的计算方法可以分为频域分析和时域分析2类。工程抗震中常用的等效线性化方法就是一种频域分析方法。等效线性化方法的基本思想就是在总体动力学效应大致相当的意义上，用一个等效剪切模量及等效阻尼比替代所有不同应变幅值下的剪切模量和阻尼比，将非线性问题转化为线性问题来求解。土层地震反应分析的另一种计算方法为时域分析方法。时域分析方法是将场地简化为剪切运动的多自由度集中质量体系，每个质量单元间通过弹簧和阻尼器连接，通过建立刚度、阻尼和质量矩阵求解各单元的动力反应。目前，基于频域和时域计算原理开发的土层地震反应分析软件有SHAKE、LSSRLI-1、SOILQUAKE、DEEPSOIL等[1]。这些软件的开发和应用基本解决了一维土层非线性地震反应的计算问题，在硬土场地上的计算结果可靠且与实测地震记录相接近。

在对复杂工程场地进行地震反应分析时，借助于数值模拟软件的计算方法得到了广泛地使用。陈健云等[2]通过FLAC3D对浅埋软土场地地铁车站地震荷载作用下的动力响应进行数值研究。边金等[3]借助FLAC3D对地铁车站进行了动力分析。得益于自由场边界和多种动力荷载输入方式，FLAC3D在一维黏弹性土层场地地震反应的模拟上可以得到合理准确的结果，与传统地震反应分析软件具有良好的一致性。Mejia等[4]利用SHAKE程序与FLAC3D程序进行了半无限弹性土层模型的地震反应计算，结果表明FLAC3D的计算结果与目标曲线非常接近。苏建锋[5]验证了用FLAC3D软件方法进行土层地震反应计算的可行性，与等效线性化程序RSLEIBM计算结果的对比表明，在输入地震动强度较小时，2种方法计算得到的结果差别不大。FLAC3D软件提供了瑞利阻尼、局部阻尼和滞后阻尼3种阻尼形式来模拟土体的滞回特征，但与前2种阻尼相比，目前利用滞回阻尼的完全非线性分析在实际工程的应用还较少，与其他一维土层地震反应分析软件计算结果的对比特别是与实测地震记录的对比还亟待研究。

为了分析FLAC3D在一维土层地震反应分析中与现行主流地震计算程序之间的异同点，理解和掌握FLAC3D在非线性分析中阻尼参数的选择，依据日本国家地球科学与防灾研究所建立的KiK-net强震动观测网提供的硬土场地地震记录和地质数据，对DEEPSOIL频域、时域分析和 FLAC3D动力分析进行对比研究。通过对比DEEPSOIL和FLAC3D在不同类型场地和不同地震荷载输入条件下计算的地表地震动峰值和反应谱，检验FLAC3D在硬土场地上进行一维土层地震反应分析的适用性，为复杂场地的地震模拟提供依据。

1 土层地震分析软件简介

1.1 DEEPSOIL简介

DEEPSOIL软件是由美国伊利诺伊大学香槟分校( UIUC)专家研发的一维土层地震响应计算程序，可以在频域(线性和等效线性)和时域(线性和非线性)进行分析[6]。DEEPSOIL有一个图形用户界面，可以通过选择适当的阻尼公式和非线性土层参数，借助自动化的拟合曲线程序计算土体的非线性行为。DEEPSOIL中包含线性、等效线性和时域非线性3种分析方法，其中时域非线性分析用Newmark β方法(隐式)或Heun方法(显式)求解时域运动方程。本文所使用的DEEPSOIL版本是DEEPSOIL v7.0.26 Released版。

1.2 FLAC3D简介

FLAC3D是Itasca公司开发的有限差分数值模拟软件，通过动力模块实现三维土层模型的地震反应模拟，在岩土领域中得到了广泛运用。FLAC的动力计算采用完全非线性方法，除了常用的瑞利阻尼和局部阻尼外，还可以采用滞后阻尼的形式嵌入到程序当中的来模拟土体的滞回现象。由于滞后阻尼中的控制参数是通过对动剪切模量比随剪应变的关系曲线进行拟合得到的，因此在应用中可以很好地与等效线性方法形成对应关系。

2 计算土层模型及输入

2.1 简化土层剖面及输入

图1 简化土层场地剖面Fig.1 Physical properties and shear wave velocity of subsoil profiles

为对比DEEPSOIL频域、时域分析与FLAC3D动力分析的异同点，按非线性程度设计了3个简化土层剖面，以内部运动和露头运动2种方式输入地震波。所设计的简化土层剖面如图 1所示，3个土层剖面分别为非线性程度依次递增的等波速土层剖面、变波速土层剖面和非线性土层剖面。

等波速土层剖面的土层厚度H设为30 m，土层剪切波速设为200 m/s，土层重度γ设为18 kN/m3,基岩层剪切波速设为500 m/s，基岩重度γ设为22 kN/m3,阻尼比λ设为5%，土层卓越频率为1.667 Hz。

变波速土层剖面采用与等波速剖面相同的土层厚度和土层密度。由于实际工程中土层剪切模量通常与有效平均应力呈正相关，因此变波速土层剖面中设定土层剪切波速随埋深的增加而增大。参考文献[7]提出的描述剪切波速与深度之间拟合关系公式，变波速土层剖面的剪切波速Us计算公式为

图2 剪切模量与阻尼比拟合曲线Fig.2 Shear modulus and damping ratio curves

vs=a+bz，

(1)

式中：a、b为拟合参数；z为土层深度。为使变波速土层剖面与等波速土层剖面具有相同的卓越频率，拟合参数分别取值a=135，b=4.71。

非线性土层剖面采用与变波速土层剖面相同的剪切波速设定，通过动剪切模量比和阻尼比随剪应变衰减的关系曲线(简称模量衰减曲线)体现土体在动力作用下的非线性现象。为了描述土体在剪切应变增大时剪切模量衰减和阻尼比增加的现象，Seed利用模量衰减曲线概括了土体在动力作用下呈现的非线性，并给出了不同类型土壤的模量比和阻尼比随剪应变的实验值[8]。非线性土层剖面中采用Seed提出的砂土模量比和阻尼比实验值见表 1。DEEPSOIL和FLAC3D分别采用了不同的方式来拟合模量衰减曲线。DEEPSOIL程序中内设的修正双曲线模型允许通过适当定义模型参数的方式来考虑土体的非线性行为，在选择Seed砂土模量衰减实验值后可以通过拟合曲线程序自动填入非线性参数[9]。FLAC中采用滞后阻尼的方式模拟岩土材料的滞后特性，通过输入拟合参数的方式将模量衰减曲线引入模型动力计算中。采用4参数模型(Sig4)输入FLAC3D中的滞后阻尼，计算使用的阻尼参数分别为1.001，-0.609，-1.195和-0.001 2。DEEPSOIL和FLAC3D对Seed砂土实验值的拟合效果如图 2所示，拟合参数见表 1。从图 2中可以看出，为了同时满足模量比和阻尼比的拟合精度，DEEPSOIL拟合的剪切模量比在小应变时略小于目标值，而在大应变时略大于目标值；阻尼比在应变范围内均略大于目标值。FLAC3D的拟合结果与DEEPSOIL基本一致。

表1 土层非线性参数Tab.1 Nonlinear calculation data of soil layers

DEEPSOIL中动力荷载输入的方式分为内部运动和露头运动2种。DEEPSOIL手册中建议刚性基岩模型对应土层内部记录的地震波，而柔性基岩模型对应露头地震波。FLAC3D中动力荷载的输入也分为这2种方式：当模型为刚性基岩时，直接在模型底部网格施加加速度时程；当模型为柔性基岩时，通过安静边界和应力时程的方式施加地震荷载。

选择1940年Imperial Valley地震中 El Centro 记录的南北分量作为输入波。该地震记录共有1 560个数据点，采样频率为50 Hz，在2.04 s时达到水平加速度峰值0.296 g，主导频率为1.172 Hz，有效持续时间为31.16 s。计算时将加速度时程调幅至0.1 g作为输入。调幅后El Centro地震波的加速度时程和傅里叶幅值谱如图 3所示。

青浦区是上海市唯一与江苏、浙江接壤的郊区，加强流域联动合作显得尤为重要。为了解决边界水事纠纷，确保一方平安，青浦区积极开展与流域机构及兄弟省市的工作交流，初步与江苏、浙江相邻市县建立了水务执法联席会议制度，2011年5月18日会同苏州市水利局、嘉兴市水利局、上海水务行政执法总队水利支队、江苏省太湖地区水利工程管理处等单位共同开展太湖流域“一湖两河”联合巡查活动，共同维护流域正常水事秩序及防汛安全。此外，在省市边界河湖治理、供水工程管理等方面能够积极沟通、主动协调。2012年3月，积极做好了浙江嘉善县作为紧急备用水源的省际湖泊长白荡相关协调工作。

2.2 台站场地选择及输入

李兆焱等[10]对巨厚场地上地震记录进行了地震反应分析，指出虽然DEEPSOIL在强震中表现不理想，但在较弱地震动下，DEEPSOIL计算的地表 PGA与加速度反应谱均接近实际情况。王鸾等[11]指出在中强地震动和较弱地震动情况下，DEEPSOIL在硬场地上的计算结果与实测反应谱较为相近。基于以上结论，从KiK-net地震台网中选取2处硬土场地，筛选地表峰值在0.05 g左右的地震记录用于计算和分析。

KiK-net地震台网由日本国家地球科学与防灾研究所建立，包含日本地区上百个地震台站的记录，可以根据时间、震级、地表峰值、台站名称等信息来筛选地震波，同时部分台站在井下基岩与地表均设立了加速度传感器，可以获得包括场地土层剖面、地震波输入与地表响应的完整数据[12]。在所选择的硬土场地IWTH12台站和IBRH11台站的地震记录中，分别选择了地表峰值为0.056 g和0.053 g的2次地震用于计算。将地震记录中的地表加速度记录作为实测值，井下加速度记录作为刚性基岩的内部运动输入波。同时挑选距离计算土层距离较近且坐落于坚硬岩层上的IWTH02台站和IBRH12台站作为露头运动参考点，取参考点所记录同期地震的地表加速度记录幅值的一半作为计算场地的露头运动输入波。对有2个方向的地震动记录选取EW方向的记录参加计算。由于KiK-net中未提供场地土层剖面的密度信息，参考文献[13]中提出的压缩波速vp与土体密度ρ的统计关系式进行计算:

ρ=0.31×vp0.25，

(2)

式中vp为压缩波速。

计算台站场地的台站场地土层剖面和剪切波速分布图如图4所示，台站场地及地震波基本信息见表 2。

(a) IWTH12台站

根据场地剪切波速分布图及土层剖面图，分别选取深度-74.00、-30.00 m作为地震动输入界面。FLAC3D土层反应模型的单元网格尺寸设为1.00 m，采用自由场边界。土层模型所需土层参数如剪切模量G和体积模量K由剪切波速进行换算：

G=ρvs，

(3)

(4)

式中v为泊松比。

表2 台站场地及地震波基本信息Tab.2 Basic information of hard site stations for contrast test

DEEPSOIL中提供了各类土的非线性参数。根据KiK-net公开资料中的土性描述来确定不同土类的非线性曲线。土层非线性参数见表3。

表3 土层非线性参数Tab.3 Nonlinear calculation data of soil layers

(5)

式中：Sa1、Sa2分别代表计算加速度反应谱和实测加速度反应谱；lgT为反应谱中周期从0.01～10 s的对数。易知反应谱计算值越接近于实测值，则谱面积比SR越接近于1。

3 计算结果与分析

3.1 简化土层剖面对比分析

在土层地震反应分析结果的多个指标当中，加速度反应谱能够综合反映地震荷载的强弱和频谱特性，是工程结构抗震中常用的指标之一。通过FLAC3D动力计算、DEEPSOIL频域和时域分析3种方法计算了简化土层剖面地表加速度反应谱，以此来对比分析3种方法计算结果的异同。各类型场地计算的地表加速度反应谱如图5—7所示。

从图中可以明显看出，在简化土层的等波速剖面和变波速剖面上，FLAC3D模型的地表加速度反应谱计算结果与DEEPSOIL频域及时域分析结果基本一致，说明所建立的FLAC3D模型对粘弹性土层地表响应的计算结果是准确的。在简化土层的非线性剖面中，FLAC3D模型的计算结果与DEEPSOIL时域分析结果较为接近，而与DEEPSOIL频域分析结果在0.3～0.6 s周期范围内存在一定差异，其余周期范围变化趋势基本一致。当输入地震波的主导频率与场地卓越频率相近时，输入地震波中的强振分量将在频域计算中被放大，被称为虚假共振现象。DEEPSOIL频域分析的计算结果在0.3～0.6 s周期范围内的谱峰值大于FLAC3D与DEEPSOIL时域分析结果，可能是由上述虚假共振现象造成的。

土层的地表反应受地震动输入方式的影响。在同一类型的场地上，图5—7中内部运动和露头运动的计算结果表现出明显的差异，这体现了基岩层与土层波阻抗比对地表地震反应的影响。另外通过不同类型土层剖面之间的对比也可以看出，土层非线性程度对地表响应有较大的影响，随着土层非线性程度的增加，地表反应谱的形态、峰值以及峰值所在的周期都存在明显差异。

(a) 内部运动

不同场地类型地表地震动幅值(peak ground acceleration, PGA)见表 4。由表 4可以看出，在相同场地类型和地震波输入方式下，FLAC3D与DEEPSOIL的计算结果较为接近。随着土层非线性程度的增加，采用露头运动输入的柔性基岩场地的地表PGA变化不大，而刚性基岩场地的地表PGA变化明显。

表4 不同场地类型地表PGA计算值对比Tab.4 Surface PGA of various layers

3.2 台站场地对比分析

IWTH12台站和IBRH11台站所在的场地等效剪切波速分别为312、242 m/s，覆盖层厚度分别为20、30 m，属于我国规范中的II类场地。在2个台站数据中挑选2组地震波输入来计算反应谱，得到结果如图8(a)—11(a)所示。同时挑选距离计算土层距离较近且坐落于坚硬岩层上的台站IWTH02和IBRH12作为柔性基岩计算的露头运动参考点，结果如图 8(b)—11(b)所示。为了对比线性和非线性土层地震分析方法对计算结果的影响，分别按照阻尼比为5%的黏弹性场地和引入模量衰减曲线的非线性场地进行计算，不同计算方法反应谱面积比对比见表5。

(a) 内部运动

当输入类型为内部运动时，从图 8(a)可以看出，DEEPSOIL计算的黏弹性场地地表反应谱曲线与实测值在形状上较为相似，但在谱峰值上存在明显差异。对比图 8(a)和图 9(a)可以看出，在引入模量衰减曲线后，非线性场地的计算结果与实测值在峰值和形状上都较黏弹性场地更为相近。IBRH11台站也呈现了相似的现象。同时，DEEPSOIL时域和频域计算结果在0.03～0.1 s周期范围内一定程度上再现了实测地震记录中的高频成分。由表5可见，谱面积比SR同样显示了这一点，与黏弹性场地相比，当考虑土体的非线性时，IWTH12台站计算的谱面积比从50%增加至60%左右，IBRH11台站计算的谱面积比从50%增加至80%左右，表明了输入地震波类型为内部运动时，考虑土体的非线性的计算结果与实测值更为接近。

当输入类型为露头运动时，对比图 8(b)和图 10(b)可以看出，IWTH12台站的计算结果与实测值较为接近，而IBRH11台站的计算结果与实测值有较明显的差异。这可能是由于所选择的露头运动参考点与计算场地的距离不同导致的。类似地，与黏弹性场地的计算结果相比，非线性场地的计算结果在0.1～0.3 s周期范围内与实测值更为相似。由表5可见，谱面积比SR显示，IWTH12台站非线性反应谱面积比从45%增加至50%左右，IBRH11台站反应谱面积比从40%增加至60%左右，表明输入地震波类型为露头运动时，非线性场地的计算效果要略优于黏弹性场地。

地表PGA实测值与计算值对比见表 6。由表 6可以看出，FLAC3D的计算结果与DEEPSOIL十分接近，但均明显小于实测地表PGA。非线性场地的计算结果要比黏弹性场地更为接近实测值，但仍然有明显差距。