高面板堆石坝三维真非线性动力分析方法及模型试验

1、高面板堆石坝三维真非线性地震反应分析方法及模型试验验证基金项目：“九五”国家科技攻关项目(96-221-02-03-02)和国家电力青年基金项目(97Q02和SPQKJ06)作者简介：赵剑明（1970），男，山东肥城人，博士/博士后，高级工程师，岩土及土工抗震研究*赵剑明 汪闻韶 常亚屏 (中国水利水电科学研究院岩土工程研究所，北京，100044)摘要基于所建立的土石料三维粘弹塑性动力本构模型，并采用新型三维各向异性有厚度薄单元来模拟面板和堆石的接触面特性，结合考虑孔隙水压力的消散和扩散的有效应力分析方法，开发了高面板堆石坝地震反应分析的三维真非线性有效应力分析方法。通过典型坝例分析，比较了真

2、非线性分析与等效线性分析结果的差异，表明真非线性方法较真实地反应了结构的地震反应，而且能够直接计算出坝体的残余变形，在理论上更为合理。采用模型坝土石料低应力状态下实测的静力和动力特性参数，对三维模型坝进行了地震反应分析，并与相应的大型振动台模型试验结果进行了对比分析，两者结果相当一致，验证了本研究所提出的计算方法的可靠性。关键词：高面板堆石坝；动力本构模型；地震反应；真非线性；模型试验1 前言混凝土面板堆石坝以其安全性、经济性、适应性等优点受到国内外坝工界的普遍重视。目前相当多的面板堆石坝位于强震区，使得面板堆石坝抗震问题研究的重要性和迫切性越来越突出。面板堆石坝地震反应的二维分析已有了较多的

3、计算方法和相应程序，但由于面板堆石坝多修建在狭谷之中，具有明显的三维效应，仅做二维分析是不够的，因此面板堆石坝地震反应的三维动力分析成为必要；随着计算技术和计算机的快速发展，也使得三维分析成为可能，所以近期面板堆石坝的三维动力分析引起了国内外学者的重视。地震反应分析方法从基于的本构模型来分可分为两大类，一类是基于等价粘弹性模型的等效线性分析方法，另一类是基于（粘）弹塑性模型的真非线性分析方法。等价粘弹性模型概念明确，使用方便，在参数的确定和应用方面积累了较丰富的试验资料和工程经验，目前应用较多。然而，等价粘弹性模型的缺陷是明显的：等价粘弹性模型在加荷与卸荷时模量相同，因而不能直接计算土体在周期

4、荷载连续作用下发生的残余变形；不能考虑应力路径的影响；不能考虑土的各向异性；较大应变时误差大，偏于不安全。由此可见，基于等价粘弹性模型的等效线性分析方法得到的地震反应并不是土体真实的地震反应，要想得到土体真实的地震反应，应采用基于（粘）弹塑性模型的真非线性分析方法。2 高面板堆石坝三维真非线性地震反应分析方法2.1 土石料的三维真非线性动力本构模型“八五”期间，汪闻韶院士和李万红1建立了一个二维的适用于土石坝的粘弹塑性动力本构模型，初步应用于工程分析，取得了一些较好的结果。本研究根据面板坝的特点，结合试验资料，对该模型进行推广和完善，建立了适用于面板坝的三维粘弹塑性动力本构模型，并开发了相应三

5、维真非线性分析方法。本研究采用的粘弹塑性模型将土视为粘弹塑性变形材料，模型由初始加荷曲线、移动的骨干曲线和开放的滞回圈组成。这种真非线性模型的特点是：（1）与等效线性粘弹性模型相比，能够较好地模拟残余应变，用于动力分析可以直接计算残余变形；在动力分析中可以随时计算切线模量并进行非线性计算，这样得到的动力响应过程能够更好地接近实际情况。（2）与基于Masing准则的非线性模型相比，增加了初始加荷曲线，对剪应力比超过屈服剪应力比时的剪应力应变关系的描述较为合理；滞回圈是开放的；考虑了振动次数和初始剪应力比等对变形规律的影响。模型的数学表达式如下：初始加荷曲线：(1)骨干曲线：(2)滞回圈：(3)式

6、(2)和式(3)中，在加荷时取（）、（），在卸荷时取（）、（）。在此非线性动力模型中，骨干曲线和滞回圈的原点不断移动产生残余变形，即有：(4)以上诸式中：和为剪应力和剪应变，max为极限剪应力，max=f /Rf ；Rf为破坏比；f为破坏剪应力；' 为有效内摩擦角；' 为有效正应力；0为骨干曲线和滞回圈原点相应的剪应变或称塑性剪应变；h是以0为零点的剪应变；A、B为模型参数；DRSd为动剪应力比幅值；DRS动剪应力比，DRS=RS-RS0，RS=/'，RS0为初始剪应力比。模型参数A和B以及0可以用剪应力比控制的循环三轴试验来测定1，主要受振次、动剪应力比幅值和初始剪应

7、力比的影响比较大。模型参数A和B也可由等效线性粘弹性模型参数换算近似得到，换算原则是使两变形模型的骨干曲线重合和滞回圈包围的面积相等。2.2 单元型式堆石坝体及地基主要采用三维八结点六面体等参单元来模拟，在边界不规则处采用六结点五面体三棱柱单元来填充。混凝土面板可用薄板单元来模拟，不过根据已有的研究结果，在二维问题中，面板采用梁单元或四边形单元来模拟，计算结果相差不大，表明面板随坝体一起运动时弯曲应力不大；而用板单元模拟面板时存在着板单元与六面体等参元位移模式不协调的误差。另一方面，尽管面板厚度小，与其它两个方向的尺寸相差大，但面板单元形态规则，不至于在有限元计算中造成大的误差，面板用六面体单

8、元离散能够满足计算精度要求。基于以上原因，并考虑到计算方便，本研究选用了同类型的三维八结点六面体等参元来模拟。为模拟混凝土面板三维受力状态下的应力应变关系的非线性，并考虑计算方便，对面板采用分段线弹性模型。2.3 接触面的模拟本研究中，在文献2的基础上，建立了一种三维各向异性有厚度薄单元3来模拟面板和堆石的接触面特性，如图1所示。这种单元刚阵与一般的三维六面体等参单元在形式上是相同的。yzx图1 三维有厚度接触面薄单元接触面上的变形可以分为基本变形和破坏变形两部分。基本变形与其它土体的变形一样，不管破坏与否都是存在的；破坏变形包括滑动破坏和拉裂破坏，只有当剪应力达到抗剪强度产生了沿接触面的滑动

9、破坏，或接触面受拉产生了拉裂破坏时才存在。在三维接触面中，有三个方向的可能破坏位移：接触面法向的张裂和沿接触面的两个方向的滑移。对于图1所示的三维有厚度接触面薄单元，则：(5)2.4 接缝的模拟为了模拟面板周边缝和垂直缝的特性，在三维有限元计算中设置了无厚度的六面体缝间连接单元4。2.5 动水压力关于动水压力，本文采用了附加质量法，具体是用一种广义边界元法来处理库水无界性的问题，给出单位加速度时，不同坝坡条件下的坝面动水压力分布系数，在计算时将其转化为相应的附加质量进行分析。2.6 基于真非线性模型的动力反应分析算法鉴于本文采用的真非线性粘弹塑性模型的特点，为了更有效地进行动力反应分析，本文采

10、用了增量法和全量法交替进行的算法以控制增量法的误差积累。根据非线性粘弹塑性模型及有限元原理，推导出结构的增量和全量方程分别为：（6）（7）式中，和分别为结点位移、速度和加速度，为弹性位移，代表增量，M为质量矩阵，Ct和Cs分别为切线和割线阻尼矩阵，Kt和Ks分别为切线和割线刚度矩阵，Fa为地震力，Fe为应力超过强度时加以修正的等价结点力（超越力）。具体求解按增量步进行。对每一增量步，先求解增量方程（6），然后如果为奇数增量步，则在假定不变的条件下，由全量方程（7）计算弹性位移；如果为偶数增量步，则在假定不变的条件下计算加速度，并用此加速度校正方程（6）中的，以减少用增量法解方程产生的误差积累。

11、尽管混凝土面板坝中坝体大部分是非饱和的，但对于建造在深厚覆盖层上的面板坝以及下游水位较高时，还是必须考虑振动孔隙水压力的影响。因此，本文采用了考虑孔压消散和扩散的三维非线性有效应力地震反应分析方法5。2.7 真非线性分析与等效线性分析的若干比较为了验证本文提出的真非线性分析方法的特点，针对一修筑于梯形河谷中坝高100m的典型面板坝，分别采用等效线性分析方法和真非线性分析方法进行了地震反应分析6，在顺河向输入El Centro波，最大加速度取为0.2g。Disp.(cm)等效线性分析真非线性分析图2 典型结点动位移时程曲线图3 典型单元动剪应变时程曲线Strain(10-3)等效线性分析真非线性

12、分析图2是两种分析方法得到的同一典型结点(698号结点)的动位移时程曲线。图3是两种分析方法得到的同一典型单元（715号单元）的动剪应变时程曲线。由图2和图3可见，等效线性方法和真非线性方法得到的应变和位移地震反应有着明显的区别：等效线性分析得出的动应变和位移围绕零点振动，没有偏移，无残余变形产生；真非线性分析得出的动应变和位移在振动过程中偏离零点，产生残余变形，并且地震过程中残余变形不断积累和增长。可见，本文提出的面板坝真非线性动力分析方法和等效非线性方法在概念上有着本质的区别，在计算结果上存在差异，真非线性方法较真实地反应了结构的地震反应，而且能够直接计算出坝体的残余变形，在理论上更为合理

13、。3 三维真非线性地震反应分析方法的模型试验验证3.1 模型试验及模型坝计算网格目前，面板堆石坝实际震害和记录资料很少，特别是缺乏强震区高面板堆石坝的资料。为了研究高面板堆石坝的动力特性和在地震作用下的动力反应性状，同时为验证和改进计算方法及程序，中国水科院利用具有世界先进水平的5×5m大型三向模拟地震振动台，对实际面板坝进行了三维大型振动台模型试验6。模型试验中采用的三维模型坝的几何比尺为123.5，坝高100cm，坝轴线长315.44cm，上游坝坡坡比为1：1.55，下游坝坡设三级马道，分别宽2.43cm，坝坡呈上缓下陡，坝坡比分别为1：1.50、1：1.40、1：1.40和1：

14、1.40。在模型试验中，对模型坝分别输入了人工合成原波（简称人工原波）、压缩人工合成原波（简称压缩人工波）和压缩松潘波，地震波特性及波形曲线见文献6。研究内容包括在不同地震波作用下模型坝的反应加速度、面板应力和应变、坝体变形等地震反应性状。为了进一步验证本研究提出的三维真非线性分析方法的可靠性，采用模型土石料低应力状态下实测的静力和动力特性参数，具体参数见文献6，对模型试验中采用的三维整体模型进行了三维真非线性地震反应分析，并与相应的试验结果作了比较验证。在对模型坝进行三维真非线性地震反应分析中，整个模型坝共划分了1996个结点和1583个单元，包括面板单元、坝体单元、接触面单元和接缝单元。模

15、型坝横剖面和纵向单元划分情况分别如图4a和图4b所示。x、y及z分别表示顺河向、坝轴向及竖直向。输入地震波采用了试验输入的地震波，即分别为人工原波、压缩人工波和压缩松潘波。计算中的输入形式和方式与试验中的保持一致，动力分析中的时间步长取为0.002s。100.00155.01152.98面板xz0+69.040+384.48yz图4 三维模型坝单元剖分图 (cm)图4b 纵向（上游立视）(cm)图4a 横剖面（0+234.67）3.2 动力分析结果 本研究计算了各种地震波作用下，三维模型坝的动力反应，并与试验结果进行了比较，取得了较为一致的结论，下面给出其中的典型结果。加速度反应 图5是三维模

16、型坝在压缩人工地震波单向（x）作用下，主断面沿中心线、上游坡和下游坡加速度放大倍数分布情况。图6是三维模型坝在不同地震波作用下坝顶（中心线）加速度放大倍数与输入地震加速度最大幅值的关系。图5 三维模型坝在压缩人工地震波单向（x）作用下，主断面沿中心线、上游坡和下游坡加速度放大倍数分布 B1109x1(0.110) 计算结果 B1109x2(0.186) 计算结果 B1109x7(0.554) 计算结果 B1109x1(0.110) 试验结果 B1109x2(0.186) 试验结果 B1109x7(0.554) 试验结果 中心线00.00.51.0Ratioh/H12345上游坡0123450.

17、00.51.0h/HRatio02468100.00.51.0h/HRatio下游坡图6 三维模型坝中心线上坝顶加速度放大倍数与输入加速度最大幅值的关系0.00.10.20.30.40.50.60.7Amax(g)0.02.04.06.08.0Ratio 1109x单向输入 计算结果 1109x单向输入 试验结果 SP1109x 单向输入 计算结果 SP1109x 单向输入 试验结果 人工原波x单向输入 计算结果 人工原波x单向输入 试验结果混凝土面板应变和应力反应 图7是三维模型坝在双向地震波作用下面板顺坡向应变沿坝高分布情况。三维模型坝在单向（x）地震波作用下面板应变与输入加速度最大值的关

18、系如图8所示。图9是三维模型坝在三向地震波作用下面板最大主应力沿坝高的分布情况。图8 三维模型坝在单向（x）地震波作用下面板应变与输入加速度最大值的关系Strain (10-4)Amax (g) 0.43H高程处x向应变 计算 0.43H高程处x向应变 试验 0.70H高程处x向应变 计算 0.70H高程处x向应变 试验 0.43H高程处y向应变 计算 0.43H高程处y向应变 试验图7 三维模型坝在双向地震波作用下面板顺坡向应变沿坝高分布情况Strainxmax (10-4) 1109xy1(0.127, 0.095) 计算 1109xy1(0.127, 0.095) 试验 1109xy2(

19、0.178, 0.146) 计算 1109xy2(0.178, 0.146) 试验 1109xz1(0.095, 0.076) 计算 1109xz1(0.095, 0.076) 试验 1109xz2(0.164, 0.117) 计算 1109xz2(0.164, 0.117) 试验0.000.050.100.150.200.250.300.00.20.40.60.81.0h/HStress 1max (100kPa) 1109xyz1(0.139, 0.143, 0.146) 计算 1109xyz1(0.139, 0.143, 0.146) 试验 1109xyz2(0.229, 0.227,

20、0.151) 计算 1109xyz2(0.229, 0.227, 0.151) 试验 1109xyz2(0.315, 0.298, 0.193) 计算 1109xyz2(0.315, 0.298, 0.193) 试验0.000.050.100.150.200.250.300.350.400.00.20.40.60.81.0h/H图9 三维模型坝在三向地震波作用下面板最大主应力沿坝高分布情况坝体残余变形 本研究计算了三维模型坝在各种地震波作用下的地震残余变形，其中整理得到的三维模型坝在地震动结束后坝顶中心线沉降变形沿坝轴线的分布情况如图10所示。图10 三维模型坝在地震动结束后坝顶中心线沉降变形

21、沿坝轴线的分布0.0040.0080.00120.00160.00200.00240.00280.00320.00B(cm)0.000.200.400.600.801.00S(mm) 1109x2(0.149) 计算 1109x2(0.149) 试验 1109xz2(0.164, 0.117) 计算 1109xz2(0.164, 0.117) 试验 1109xyz4(0.315, 0.298, 0.193) 计算 1109xyz4(0.315, 0.298, 0.193) 试验 b1109x8(0.664) 计算 b1109x8(0.664) 试验由各图可见，计算结果与试验结果相当一致。具体结

22、论见下一小节。3.3对比验证及小结本文应用本研究提出的三维真非线性地震反应分析方法，采用模型土石料低应力状态下实测的静力和动力特性参数，对三维模型坝进行了地震反应分析，并将计算结果对试验结果进行了对比分析，分析结果表明，计算得到的地震作用下坝体的加速度反应、面板的应力和应变及坝体残余变形与试验结果相当一致。从而，在计算和试验相互印证的同时，也进一步验证了本研究提出的三维真非线性地震反应分析方法的可靠性，并可以得出如下的一些结论：1）对相同高程，下游坡面上的加速度反应比上游面板上及中心线上的加速度反应大，表现出较明显的坝体表层加速度放大效应，在工程上应注意加强下游坡的抗震保护，尤其是上部表面；而

23、上游面加速度反应较下游为小，体现了面板的约束作用。2）地震波振幅对加速度放大倍数有较大影响，对于整个坝高，加速度放大倍数有随输入加速度最大幅值Agxmax的增大而减小的趋势，而且使加速度分布也有所变化。3）地震波类型等对加速度放大倍数分布有较大影响，其中人工原波得到的放大倍数显著地小于经过时间比尺压缩的人工波所得到的放大倍数，这主要是由于人工原波卓越频率远小于坝体基本频率的缘故；同样时间比尺Ct的压缩人工波与压缩松潘波由于它们卓越频率段不同，在模型上产生的加速度放大倍数和分布也有所不同。4）计算得到面板顺坡向应变及主应变和主应力最大值发生位置虽然有所变化，但最大值位置基本发生在0.50.7H的

24、范围内。5）地震波类型和加速度幅值对面板应变和应力有重要影响，其影响规律与影响加速度的规律类似。6）坝顶地震沉降在坝轴向中部较大，但总体上地震沉陷量并不大，在本研究输入的各种地震波作用下均不超过坝高的0.1%，表明面板堆石坝具有较好的抗震性能。4 结论本文基于所建立的土石料三维粘弹塑性动力本构模型，并采用新型三维各向异性有厚度薄单元来模拟面板和堆石的接触面特性，结合考虑孔隙水压力的消散和扩散的有效应力分析方法，开发了高面板堆石坝地震反应分析的三维真非线性有效应力分析方法。通过典型坝例分析，比较了真非线性分析与等效线性分析结果的差异，表明真非线性方法较真实地反应了结构的地震反应，而且能够直接计算

25、出坝体的残余变形，在理论上更为合理。采用模型坝土石料低应力状态下实测的静力和动力特性参数，对三维模型坝进行了地震反应分析，并与相应的大型振动台模型试验结果进行了对比分析，两者结果相当一致，验证了本研究所提出的计算方法的可靠性。根据模型试验和计算分析成果，得出了一些有关面板堆石坝地震反应性状、抗震设计方法及抗震措施方面的研究结论，可为工程建设提供参考。参考文献1 李万红、汪闻韶，无粘性土动力剪应变模型J。水利学报，1993年第9期。2 殷宗泽、朱泓、许国华，土与结构材料接触面的变形及其数学模型J。岩土工程学报，1994，16（2）。3 赵剑明、汪闻韶，混凝土面板堆石坝面板地震反应分析J。岩石力学

26、与工程学报，2001年，20（s2）。4 顾淦臣，土石坝地震工程M。河海大学出版社，1989。5 赵剑明、汪闻韶、张崇文，土石坝振动孔压影响因素的研究J。水利学报，2000年第5期。6中国水利水电科学研究院，强震区面板坝大型振动台模型试验及动力分析R。“九五”国家科技攻关专题研究报告，2000。3-D Authentic Nonlinear Dynamic Analysis Method of High CFRD and Verification by Model Test DataZhao Jianming Wang Wenshao Chang Yaping(Geotechnical Engineering Department, IWHR, Beijing, 100044)Based on the proposed 3-D authentic nonlinear dynamic constitutive model of soil and applying a new kind of 3-D anisotropic thin element to simulate t