考虑渗流体力的某大坝地震响应仿真分析

韩 彰,苏怀智,李 慧

(1.河海大学水文水资源与水利工程科学国家重点实验室,江苏 南京 210098;2.河海大学水利水电学院,江苏 南京 210098)

近年来,国家规划、建造了一批土石坝[1-2],这些土石坝多处于高地震烈度区,需进行地震动力响应分析,以确保在建和已建土石坝的安全[3-5]。目前,土石坝安全及动力响应分析的研究手段主要有极限平衡法、物理模型试验、数值仿真、剪切楔法及集中质量法等。沈振中等[6]推导了满足力和力矩平衡的坝坡稳定极限平衡水平条分法的计算公式;刘小生等[7]开展了面板堆石坝振动模型试验及动力分析研究;窦兴旺等[8]开展了深覆盖层上高堆石坝振动台试验与动力数值分析验证研究。由于数值仿真能直接反映力学演化过程,且计算成本低,因此,数值仿真已被广泛应用于土石坝的动力响应分析,如窦兴旺等[9]研究了人工边界法在土石坝动力分析中的应用;杨秀竹等[10]开展了地震力作用下软基土石坝的动力响应分析;丁陆军[11]基于有限元法和土体的等价线性分析,运用三维有限元动力分析程序对某面板堆石坝进行了动力分析,研究了该坝的整体受力情况; 曹学兴等[12]改进了传统的Hardin-Drnevich模型,提出了可以考虑土石料动力特性参数围压依赖性的改进动力本构模型;庞林祥等[2]提出了适用于高土石坝的地震动力响应分析计算方法;孔宪京等[13-14]、Piao等[15]、赵剑明等[16]、杨昕光等[17]都开展了相关研究。但上述研究主要集中在土石坝地形条件、动力本构模型、抗震性能、计算方法等方面,而考虑渗流体力的土石坝动力仿真研究国内外还鲜有报道。针对某水库大坝混凝土连接坝段、均质壤土挡水坝段静动力计算中须考虑渗透体力的问题,本文采用等效线性黏弹性模型,提出了一种在静动力计算中渗流体力的施加方法,并采用MSC.Marc有限元软件,通过Fortran语言编写用户子程序，对该坝地震响应过程开展了仿真分析。

1 工程概况

某水库工程(从左至右)主要由混凝土挡水坝段、混凝土取水坝段、泄洪排沙底孔坝段、表孔溢流坝段、混凝土连接坝段、均质壤土挡水坝段等组成。混凝土挡水坝段为现浇C20混凝土重力坝,坝体断面上游为垂直面,下游坝坡1∶0.8,分为2个小坝段,长度分别为26 m和33 m,总长59.0 m；坝顶高程1 068.50 m,坝基坐落在弱风化岩体上,坝基最低建基面高程1037.0 m,向左岸侧逐渐抬高至1061.0 m,最大坝高31.5 m。混凝土连接坝段的连接采用插入式,连接坝段分为4个坝段,每段坝段长25 m,其中半插入段长75 m,插入段长25 m,总长100 m；坝基坐落在弱风化基岩上,建基面高程1029.00 m,最大坝高39.50 m。均质壤土挡水坝段上游坝坡自上而下分别为1∶3.0和1∶3.5,下游坝坡自上而下分别为1∶2.75和1∶3.0,上游坝坡采用现浇C20混凝土板护面,混凝土板厚200～300 mm，混凝土板下铺设粒径20～40 mm的砂砾石垫层,厚800 mm;下游坝坡高程1045.0 m以上坡面采用预制C20混凝土网格草皮护坡。该区地震震动峰值加速度0.10g,对应地震基本烈度为Ⅶ度。根据该水库坝址工程场地地震安全性评价报告,确定该坝址场地50年内超越概率为10%的地震烈度值为7.23,超越概率为2%的地震烈度值为7.88,地震持续时间取30 s。

2 数值模拟方法

采用动力计算本构模型对该水库均质壤土挡水坝段及其混凝土连接坝段非线性有限元动力特性进行数值模拟。

2.1 动力计算本构模型

a. 均质壤土挡水坝段坝体土料和棱体排水料的动力计算模型。地震荷载是一种非等幅等周期的不规则荷载,在一次地震中,土料将经历数十次甚至上百次卸载和再加载的过程,并且它们之间无规律可循,为解决此问题,比较常用的方法是应用Masing规则,制定一个应力应变关系的骨架曲线,在此基础上,可采用双线性、黏弹性和弹塑性等多种本构模型[18-20]。由于等效线性黏弹性模型简单,结果精确,动力计算分析采用等效线性黏弹性模型,即假定坝体土料为黏弹性体,采用等效剪切模量G和等效阻尼比λ这两个参数来反映土体动应力应变关系的非线性和滞后性两个基本特征,并表示为剪切模量和阻尼比与动剪应变之间的关系。

b. 混凝土连接坝段计算模型。混凝土连接坝段动力计算分析时采用线性弹性模型。

c. 接触面的动力计算模型。接触面单元的动力参数采用河海大学的试验成果[18]。接触面单元的动力模型剪切劲度Kc与动剪应变γd、阻尼比λc的关系分别为

(1)

(2)

其中Kcmax=Γσn0.7τf=σntanδ

式中:σn为接触面单元的法向应力;δ为接触面的摩擦角;λcmax为最大阻尼比;Β、Γ为试验参数。

2.2 坝体结构网格、边界条件与荷载

图1 水库坝体三维有限元模型

图2 桩号位置示意图

根据初步设计图、地质勘测报告等资料,考虑坝体分区、施工程序及加载过程,对坝体进行剖分,建立三维有限元模型。网格采用六面体和四面体结构化网格,总共剖分6 560个单元,7 417个结点,其中混凝土连接坝段716个单元,均质壤土挡水坝段5 228个单元,排水棱体616个单元。坝体的三维有限元模型见图1,桩号位置见图2。动力计算中忽略坝体与基岩的动力相互作用,将坝体与基岩的接触面设为三向约束。将重力荷载、地震荷载及正常蓄水位下稳定渗流场所形成的坝体渗透压力施加于坝体。借助Abaqus商业软件进行渗流计算,为了完成Abaqus与Marc的对接,编写3个Fortran子程序处理渗流计算得到的流速,将流速转化为坡降数据,作为渗流力施加的输入数据,在Marc中以加体力边界的方式完成渗流压力的施加。编写3个子程序的目的为:①找出Abaqus与Marc结点坐标对应关系(Marc退化单元筛分),进行结点坐标转换以及Marc中8、6结点单元筛选;②找出Abaqus与Marc中单元对应关系,进行8、6、4结点单元转换;③提取渗透流速,计算出渗透坡降,导出dat输入文件。

坝体0+158断面渗透压力分布以及坝轴线断面渗透压力分布见图3和图4。蓄水期0+158断面渗透压力最大值为0.28 MPa,坝轴线断面渗透压力最大值为0.08 MPa。

图3 坝体0+158断面渗透压力分布(单位:MPa)

图4 坝轴线断面渗透压力分布(单位:MPa)

2.3 动力计算参数

通过试验测得动剪切模量比Gd/Gdmax和动阻尼比λd与动剪应变γd的关系曲线。动力计算时输入相应关系曲线的控制数据,根据应力应变值进行内插和外延取值。本工程坝料试验结果如表1所示。坝料模量系数K′=696.0,模量指数n′=0.5,固结比kc=1.0。

混凝土连接坝段动力计算分析时采用的线性弹性模型,静力弹性模量Es提高30%即为动力弹性模量Ed。接触面单元的动力模型采用河海大学试验成果得出的接触动本构模型[18]。接触面的摩擦角取32°，最大阻尼比取0.2，Β、Γ分别取2.0和22.0。

表1 不同动剪应变下坝料的动剪切模量比和动阻尼比

2.4 坝体动力非线性有限元计算步骤

在进行动力计算分析之前,必须先进行静力计算分析,以获得动力分析坝体的初始应力状态。静力分析方法比较简单,这里不再赘述。计算中采用MSC.Marc有限元软件进行建模与仿真,Marc具有较强的处理非线性问题的能力,并具有良好的二次开发接口,通过Fortran语言自编用户子程序，对该坝地震响应过程开展了仿真分析，计算的主要步骤如下:

步骤2求出土体单元的初始动剪切模量Gdmax,0,土体单元的初始阻尼比根据经验取为5%。

步骤3整个地震历程划分为若干个时段。

步骤4对每个时段的动剪切模量迭代求解。

步骤5用Willson-θ法建议的放大时间间隔ΔT=θΔt代替实际时间间隔Δt,对每个时段进行时程分析。

步骤6计算各单元的质量矩阵和刚度矩阵,形成总体质量矩阵M和刚度矩阵K,采用空间迭代法求出坝体基频ω,并计算单元阻尼矩阵,由各单元的变阻尼矩阵ce组装得到总体阻尼矩阵C。

步骤10根据求出的结点位移un+1计算各单元的动剪应变γn+1和动剪应力τn+1。

步骤11重复步骤 5～10,得到各单元在每个时段内的动剪应变γd时程。

步骤12求出各单元γd时程中的最大值γdmax，根据等效动剪应变γd eff=0.65γdmax,查Gd/Gdmax～γd和λc～γd曲线得到新的Gd和λc。

步骤13重复步骤4～12,直到前后两次采用的Gd相对误差小于10%。

步骤14重复步骤3～13,直到各个时段全部计算结束,即整个地震历程结束。

步骤15输出计算结果。

2.5 地震加速度的输入

坝体的动力反应计算需考虑“正常蓄水位+地震”工况。首先进行静力分析,并将水库水位蓄至正常蓄水位,随后施加地震荷载,忽略坝体与地基的相互动力作用,加速度直接施加于坝体与坝基接触面上,进行地震反应分析。

根据水库坝址工程场地地震安全性评价报告中的场地基岩水平加速度反应谱曲线(图5),采用人工合成法合成50年超越概率为10%的设计地震加速度,其中地震的持续时间取30 s。地震波输入方向为:x向沿原河流方向水平加速度输入;y向沿高程方向竖直加速度输入,依据NB 35047—2015《水电工程水工建筑物抗震设计规范》,将其峰值折减2/3;z向为沿坝轴方向横向加速度输入。图6为人工合成的设计地震三向加速度曲线。计算中将整个地震历程划分为30个大时段,每个大时段又划分为50个小时段,因此,积分计算的时间步长为0.02 s。

图5 场地基岩水平加速度反应谱

图6 地震加速度时程曲线

计算过程中,记录了坝体0+158断面12个结点(图7)的加速度反应和动位移反应及8个单元(图8)的最大、最小动主应力反应。

图7 0+158断面12个结点位置及编号

图8 0+158断面8个单元位置及编号

3 三维非线性动力计算结果与分析

3.1 设计地震反应

由于试验条件不允许,混凝土连接坝段和均质壤土挡水坝段永久变形计算参数参考了积石峡混凝土面板堆石坝[21]数据和相关经验。典型坝体断面(0-32断面)连接坝段最大加速度幅值和最大动位移等值线分布见图9和图10。

图9 0-32断面动位移等值线分布(单位:cm)

图10 0-32断面加速度等值线分布(单位:m/s2)

设计提供的模拟地震动加速度曲线历时达30 s,因此,在整理成果时给出了30 s的时程曲线,设计地震工况基岩输入加速度取50年超越概率为10%的动力有限元分析计算成果:①x向和y向最大加速度分别为5.0 m/s2和4.0 m/s2,分别出现在坝体最大断面下游坝顶附近和0-81断面上游坝顶附近,x向放大系数为7.70,y向放大系数为6.16;②x向和y向最大位移分别为9.0 cm和1.8 cm,分别出现在坝体最大断面坝顶附近和0-81断面上游坝顶附近;③土体第一和第三最大主应力都为200 kPa,分别出现在0-81断面上游坝底和0-32断面上游坝底;④x向和y向地震永久变形分别为0.7 cm和2.5 cm,x向出现在0-32和0-81断面下游坝顶以及坝体最大断面下游坝坡1/3坝高处附近，y向出现在坝体最大断面坝顶附近,地震沉陷量为坝高的0.06%。

图11 0-81断面接触面结点分布

由计算成果可知:①最大地震动加速度、动位移反应位于坝顶局部位置,坝顶存在明显的鞭鞘效应,需要在坝顶进行适当的抗震加固;②计算得到的最大放大系数为7.70,这可能是由于大坝很长,约束又加在两端,正如橡皮筋一样,导致计算得到的反应加速度放大倍数偏大,动力响应偏大,建议在坝顶加土工格栅。

图12 端头接触面结点分布

3.2 上下游接触面及端头接触面典型结点动位移

图11为0-81断面接触面结点分布,图12是端头接触面结点分布,由于这些典型结点的x、y、z向动位移反应历时曲线过多,这里不再附图展示。经统计分析,0-81断面接触面各结点在30 s设计地震过程中三向动位移极值范围见表2,端头接触面第1～4列各结点在30 s设计地震过程中的三向动位移范围见表3。由表2和表3可见,三向动位移反应均比较小,并且对称分布,因此动力工况下接触面工作状态基本正常。

表2 0-81断面接触面结点动位移统计

表3 端头接触面第1～4列结点动位移统计

3.3 端头接触面典型单元法向动应力

端头接触面地震历程30 s末法向动应力分布如图13所示,可见设计地震工况下,端头接触面法向动应力最大值为3.70 MPa,出现在接触面底部。经统计,端头接触面法向动应力均为压应力,因此混凝土连接坝段在设计地震工况下是稳定的。

图13 端头接触面典型单元法向动应力分布(单位:MPa)

4 结 论

a. 该大坝的最大地震动加速度、动位移反应位于坝顶局部位置,坝顶存在明显的鞭鞘效应,需要在坝顶进行适当的抗震加固。

b. 最大地震动应力反应位于混凝土连接坝段上游坝底附近位置;计算得到的最大放大系数为7.70。

c. 各分区的设计与填筑标准、坝体分层填筑方案合理,坝体抗震安全性较好。

[1] 吴中如.水工建筑物安全监控理论及其应用[M].北京:高等教育出版社,2003.

[2] 庞林祥,傅华.高土石坝的地震动力响应分析方法研究与应用[J].水利水电技术,2015(7):3-9.(PANG Linxiang,FU Hua.Study on method for analysis on seismic dynamic response of high earth-rockfill dam and application[J].Water Resources and Hydropower Engineering,2015(7):3-9.(in Chinese))

[3] 李阳,任亮,王瑞骏,等.基于ABAQUS的面板砂砾石坝三维动力有限元分析[J].水资源与水工程学报,2014(6):39-43.(LI Yang,REN Liang,WANG Ruijun,et al.Finite element analysis of 3D dynamics for concrete face rock fill dam based on ABAQUS[J].Journal of Water Resources & Water Engineering,2014(6):39-43.(in Chinese))

[4] 任亮.考虑流变特性的面板砂砾石坝静动力有限元分析[D].西安:西安理工大学,2013.

[5] 李炎隆,徐娇,李守义.沥青混凝土心墙堆石坝动力有限元分析[J].应用力学学报,2016,33(2):247-253.(LI Yanlong,XU Jiao,LI Shouyi.Finite element analysis of dynamic on rockfill dam with asphaltic concrete core wall[J].Chinese Journal of Applied Mechanic,2014(6):39-43.(in Chinese))

[6] 沈振中,崔娟,任华春.土石坝抗震稳定性分析的极限平衡水平条分法[J].南水北调与水利科技,2011,9(3):145-149(SHEN Zhenzhong,CUI Juan,REN Huachun.Quasi-static limit equilibrium horizontal slice method forseismic stability analysis of earth-rock dam[J].South to North Water Diversion and Water Science & Technology,2011,9(3):145-149(in Chinese))

[7] 刘小生,王钟宁,赵剑明,等.面板堆石坝振动模型试验及动力分析研究[J].水利学报,2002,33(2):29-35.(LIU Xiaosheng,WANG Zhongning,ZHAO Jianming,et al.Advancement of technology on shaking table model test and dynamic analysis of CFRD[J].Journal of Hydraulic Engineering,2002,33(2):29-35.(in Chinese))

[8] 窦兴旺,林义兴,夏颂佑,等.深覆盖层上高堆石坝振动台试验与动力数值分析验证研究[J].水利水电科技进展,2000,20(6):19-22.(DOU Xingwang,LIN Yixing,XIA Songyou,et al.Shaking-table test and dynamic numerical analysis of high rookfill dams constructed on deep overburdens [J].Advances in Science and Technology of Water resources,2000,20(6):19-22.(in Chinese))

[9] 窦兴旺,夏颂佑,许百立.人工边界法在土石坝动力分析中的应用[J].河海大学学报(自然科学版),2000,28(5):72-75.(DOU Xingwang,XIA Songyou,XU Baili.Application of artificial boundary method to dynamic analysis of embankment dams[J].Journal of Hohai University(Natural Sciences),2000,28(5):72-75.(in Chinese))

[10] 杨秀竹,李花珍,王星华.地震力作用下软基土石坝的动力响应分析[J].铁道科学与工程学报,2008,5(3):51-53.(YANG Xiuzhu,LI Huazhen,WANG Xinghua.Dynamical response analysis of the earth-rock dam on soft soil during earthquake[J].Journal of Railway Science and Engineering,2008,5(3):51-53.(in Chinese))

[11] 丁陆军.面板堆石坝三维动力反应分析[D].西安:西安理工大学,2005.

[12] 曹学兴,何蕴龙,熊堃,等.考虑围压效应的高土石坝动力响应分析[J].天津大学学报(自然科学与工程技术版),2014,47(2):116-123.(CAO Xuexing,HE Yunlong,XIONG Kun,et al.Dynamic response of high rockfill dam considering confining pressure effect [J].Journal of Tianjin University(Science and Technology),2014,47(2):116-123.(in Chinese))

[13] 孔宪京,周扬,邹德高,等.汶川地震紫坪铺面板堆石坝地震波输入研究[J].岩土力学,2012,33(7):2110-2116.(KONG Xianjing,ZHOU Yang,ZOU Degao,et al.Study of seismic wave input of Zipingpu concrete face rockfill dam during Wenchuan earthquake[J].Rock and Soil Mechanics,2012,33(7):2110-2116.(in Chinese))

[14] 孔宪京,邹德高,邓学晶.高土石坝综合抗震措施及其效果的验算[J].水利学报,2006,37(12):1489-1495.(KONG Xianjing,ZOU Degao,DENG Xuejing.Comprehensive earthquake resistant measure of high earth-rockfill dams and effectiveness verification[J].Journal of Hydraulic Engineering,2006,37(12):1489-1495.(in Chinese))

[15] PIAO R,RIPPE A,MYERS B.Earth dam liquefaction and deformation analysis using numerical modeling[J].Rock Mechanics and Geotechnical Engineering,2006,8(8):80-85.

[16] 赵剑明,温彦锋,刘小生,等.深厚覆盖层上高土石坝极限抗震能力分析[J].岩土力学,2010,31(增刊 1):41-47.(ZHAO Jianming,WEN Yanfeng,LIU Xiaosheng,et al.Study of maximum aseismic capability of high earth-rock dam on deep riverbed alluviums[J].Rock and Soil Mechanics,2010,31(Sup1):41-47.(in Chinese))

[17] 杨昕光,迟世春.土石坝坡抗震极限分析上限法[J].岩土力学,2013,34(3):721-726.(YANG Xinguang,CHI Shichun.Upper bound method for seismic stability limit analysis of earth-rock dams[J].Rock and Soil Mechanics,2013,34(3):721-726.(in Chinese))

[18] 顾淦臣.土石坝地震工程[M].南京:河海大学出版社,1989.

[19] 李子阳,谷艳昌,张磊.基于MSC.Marc的高拱坝蓄水期变形分析[J].水利水电科技进展,2008,28(4):15-19.(LI Ziyang,GU Yanchang,ZHANG Lei.MSC.Marc-based deformation analysis of high arch dam in water storage period [J].Advances in Science and Technology of Water resources,2008,28(4):15-19.(in Chinese))

[20] 杨贵,何敦明,刘汉龙,等.土石坝动力离心模型试验颗粒流数值模拟[J].河海大学学报(自然科学版),2011,39(3):296-301.(YANG Gui,HE Dunming,LIU Hanlong,et al.Advancement of technology on shaking table model test anddynamic analysis of CFRD [J].Journal of Hohai University(Natural Sciences),2011,39(3):296-301.(in Chinese))

[21] 赵剑明,常亚屏,陈宁.强震区高混凝土面板堆石坝地震残余变形与动力稳定分析[J].岩石力学与工程学报,2004,23(增刊 1):4547-4552.(ZHAO Jianming,CHANG Yaping,CHEN Ning.Study on earthquake-induced permanent deformation and dynamic stability of high concrete-faced rockfill dam [J].Chinese Journal of Rock Mechanics and Engineering,2004,23(Sup1):4547-4553.(in Chinese))