FLAC3D在土石坝地震反应分析中的应用

杨 星 ， 余 挺 ， 杨 贵

（1.河海大学岩土力学与堤坝工程教育部重点实验室，江苏 南京 210098；2.河海大学岩土工程科学研究所，江苏 南京 210098；3.中国水电顾问集团成都勘测设计研究院，四川 成都 610072）

0 引言

目前，土石坝地震反应分析普遍采用的拟静力法无法反映坝体的动力特性、变形特性、输入地震波的频率及持续时间等特性，且得到的最小安全系数不能定量地反映大坝的抗震安全度。基于等价粘弹性模型的等效线性分析的有限单元法在参数确定和工程应用方面积累了较丰富的经验，得到了广泛应用。然而，该方法的缺点是加荷和卸荷的模量相同，不能直接计算土体在周期荷载作用下的永久变形；所用割线模量在小应变时与非线性切线模量相近，但在大应变时两者相差很大，偏于不安全；塑性屈服模拟不合理以及不能计算大变形。FLAC3D(Fast Lagrangian Analysis of Continua)有限差分程序采用完全非线性动力分析方法，可以遵循任何指定的非线性本构模型；可以直接计算永久变形；可以分析岩土体在动力荷载作用下发生的大变形，具有强大的后处理能力等[1]。

本文利用FLAC3D动力弹塑性模型对土石坝进行了地震反应分析，着重讨论了用FLAC3D进行土石坝地震反应分析时地震波滤波和基线校正、动力边界条件设置、瑞利阻尼参数确定等关键问题，并用FISH语言实现了坝体初始剪切模量和体积模量随平均主应力的非线性变化。在此基础上，对土石坝进行了地震反应分析，为利用FLAC3D进行土石坝地震反应分析提供了范例。

1 计算模型的建立

1.1 几何模型和网格

某心墙堆石坝坝高100 m，坝顶宽10 m，上、下游坝坡坡比1∶2，心墙顶宽6 m，心墙上、下游坡比为1∶0.2。针对目前 FLAC3D前处理能力较弱，而ABAQUS有限元软件具有较强建模能力的特点，土石坝模型在ABAQUS有限元软件中建立。主体模型采用C3D8六面体单元，坝坡边缘处采用三棱柱单元。运用FORTRAN语言编写了模型转换程序，直接提取ABAQUS生成的数据文件（*.inp）中的节点（Node）、单元 （Element）和集合 （Set）信息， 转化成为FALC3D相应的节点（Gridpoint）、 单元（Zone）和组（Group）信息，并保存为网格数据文件（*.flac3d）。坝体几何模型和转换后的FLAC3D网格见图1。计算中没有考虑上游水荷载和坝体内水的渗流作用。

图1 坝体几何模型和有限差分网格

1.2 本构模型和材料参数

坝体堆石料采用摩尔-库伦（Mohr-Coulomb）弹塑性本构模型，当剪应力水平在弹性范围内时，通过粘滞阻尼达到能量耗散；当剪应力达到屈服强度时，通过塑性流动实现能量耗散[2]。该模型能模拟土石坝地震时的塑性剪切变形，并且由于采用了弹塑性本构模型，可以直接计算坝体的永久变形。坝体堆石料动力计算的初始剪切模量随平均有效主应力的变化关系为

式中，C为动剪切模量系数；n为动剪切模量指数；两者值由试验确定。σm为平均有效主应力；Pa为大气压。

在FLAC3D中，根据式（1）通过FISH语言编程实现初始剪切模量随平均有效主应力的非线性变化。

堆石料的动泊松比υd一般在0.3～0.4之间。动泊松比υd取为0.35，则体积模量K=3G。计算中采用的材料参数见表1。坝体初始剪切模量计算结果见图2。

1.3 地震波滤波和基线校正

模型动力计算的地震波采用1989年美国加利福尼亚洛马普列塔（Loma Prieta）地震期间列克星顿（Lexington）大坝的场地实测波，地震波峰值加速度（PGA）为0.17 g，持续时间约为40 s。输入的地震波加速度时程曲线见图3。本计算只考虑水平向地震波输入。

表1 堆石料的计算参数

图2 坝体初始剪切模量分布（单位：Pa）

图3 输入的地震波加速度时程曲线

模型的网格尺寸会影响地震波传播的准确性。库勒迈耶（Kuhlemeyer）和赖斯默（Lysmer）的研究表明[3]，要想精确描述模型中地震波的传播，模型网格尺寸 必须小于输入波最高频率对应的波长的1/10～1/8。对加速度做谱分析，其功率谱见图4。由图4分析可知，地震波95%以上的能量集中在5 Hz范围内。因此，地震波滤波的截断频率确定为5 Hz。根据截断频率和初始剪切模量，划分网格时控制模型主体网格最大尺寸为5 m。

在FLAC3D地震计算中，若对时程加速度积分得到的最终速度和位移不为0，则在动力计算结束后模型会出现继续的残余位移，此时需要对加速度时程曲线作基线校正。使用SeismoSignal地震波处理软件对上述地震波进行滤波和基线校正，修正后的速度和位移时程曲线分别见图5、6。

1.4 动力边界条件设置

图4 加速度时程的功率谱

图5 修正后的速度时程曲线

图6 未修正和修正后的位移时程曲线

在动力计算中，模型边界条件的选取是一项重要内容，因为边界上波的反射会对动力计算结果产生影响。FLAC3D提供了静态（粘滞）边界和自由场边界两种边界条件。自由场边界条件要求把模型设置得足够远，模型的范围取得足够大，才能使模型边界上波的反射尽可能的小，并且自由场边界要求模型4个侧面必须垂直。对于本文单独以坝体为研究对象的模型，只能采用静态（粘滞）边界条件。

静态（粘滞）边界是Lysmer和Kuhlemeyer提出来的[4]，通过在模型边界的法向和切向分别设置自由的阻尼器来实现吸收反射波。若在静态（粘滞）边界上施加动荷载，则只能施加时程应力。时程速度转换为时程应力的公式如下

式中，σn、σs分别为施加在静态（粘滞）边界上的法向应力和切向应力；ρ为密度；Cp和Cs为通过介质的p波和s波的波速；υn和υs为边界上的法向和切向速度分量。

1.5 阻尼的选取

FLAC3D动力计算提供了瑞利阻尼、局部阻尼和滞后阻尼等3种阻尼形式。瑞利阻尼其理论与常规动力分析方法类似，计算得到的加速度响应规律比较符合实际[1]。因此，本文动力计算采用瑞利阻尼。

在FALC3D中，使用瑞利阻尼必须确定最小中心频率fmin和最小临界阻尼比ξmin。瑞利阻尼是与频率相关的，对于土石坝这种材料组成较为复杂的结构，中心频率既不是模型的自振频率，也不是输入地震波的主频，而与两者都有关系。本文采用对模型做无阻尼动力计算的方法确定中心频率。对土石坝进行无阻尼地震动力计算，得到坝顶点x方向的时程速度（见图7）。对时程速度作傅里叶变换，得到该点时程速度的功率谱（见图8）。确定该土石坝的中心频率fmin约为1.05 Hz。对于岩土材料，临界阻尼比的范围一般是2%～5%，在FLAC3D动力计算中一般可直接选取为5%[5]。

图7 速度时程曲线

图8 速度时程功率谱

2 计算结果分析

坝顶中心点的竖向永久变形见图9。从图9可知，坝体竖向永久变形达到22.3 cm。

输入加速度和坝顶中心点地震反应加速度的时程见图10。从图10可知，坝顶加速度最大值为5.03 m/s2，与输入地震波的PGA相比，放大倍数达到2.9。在实际工程中，坝顶部位应采取抗震加固措施，并预留足够的安全超高，以保证大坝的安全。

图9 坝顶竖向永久位移时程曲线

图10 加速度时程曲线

3 结语

利用FLAC3D动力弹塑性模型对土石坝进行了地震反应分析，结果表明，在地震荷载作用下，土石坝竖向永久变形达到22.3 cm，坝顶加速度最大值为5.03 m/s2，与输入地震波的PGA相比，放大了2.9倍，计算结果符合已有工程研究规律。在实际工程中，坝顶部位应采取抗震加固措施，并预留足够的安全超高，以保证大坝的安全。

［1］陈育民，徐鼎平.FLAC/FLAC3D基础与工程实例［M］.北京：中国水利水电出版社，2009.

［2］Noorzad R,Omidvar M.Seismic displacement analysis of embankment dams with reinforced cohesive shell［J］.Soil Dynamics and Earthquake Engineering,2010,30（11）:1149-1157.

［3］Kuhlemeyer R L,Lysmer J.Finite element method accuracy for wave propagation problems［J］.Journal of Soil Mechanics and Foundations,1973,99（SM4）:421-427.

［4］Lysmer J,Kuhlemyer R L.Finite dynamic model for infinite media［J］.Journal of Engineering Mechanics,1969,95（EM4）:859-877.

［5］Chakraborty D,Choudhury D.Investigation of the behavior of tailings earth dam under seismic conditions［J］.American Journal of Engineering and Applied Sciences,2009,2（3）:559-564.