西部某水电站覆盖层地基离心机动力试验数值模拟

邹德高，刘京茂，汪玉冰，李俊超，李 多，陈 涛，王 锋

（1.大连理工大学水利工程学院，辽宁省大连市 116024；2.大连理工大学海岸和近海工程国家重点实验室，辽宁省大连市 116024；3.浙江大学软弱土与环境土工教育部重点实验室，浙江杭州市 310058；4.四川华能泸定水电有限公司，四川省成都市 610072；5.中国电建集团成都勘测设计研究院有限公司，四川省成都市 610072）

0 引言

土石坝因其良好的地基适应性，是深厚覆盖层上建坝的首选坝型[1]。根据我国水电能源战略布局的需求，我国建设和规划了一批坐落在覆盖层上的高土石坝工程，例如：瀑布沟（覆盖层深度h=70m）、长河坝（h=76.5m）、阿尔塔什（h=90m）、冶勒（h=420m）、旁多（h=424m）等[1，2]，这些高坝很多位于地震频发、地质复杂的强震区且覆盖层存在易液化土层，强震作用下覆盖层地基变形和稳定是评价深厚覆盖层上高土石坝抗震安全的重要内容。

有限元数值分析方法是目前研究深厚覆盖层上高土石坝抗震安全的重要手段。有限元数值分析主要是通过计算获得土工构筑物的地震加速度、动应力或振动孔压等变量的空间分布，从而为土工构筑物的抗震稳定性评价或地基液化分析提供依据。有限元分析土体或坝体的地震动力响应时，土的动力非线性性质可以采用等价黏弹性模型及弹塑性模型，相应地将这些本构模型与有限元结合，发展了等价线性分析方法以及弹塑性动力分析方法[3]。其中基于等价黏弹性模型的等价线性分析方法理论简单、应用方便，而且在参数确定和应用方面积累了比较丰富的试验资料和工程经验，已为工程界普遍接受；另一种基于弹塑性模型的非线性分析方法能够较好地模拟土体的实际反应，并能够直接计算土工构筑物的地震变形和孔压发展过程，在理论上更为严密和合理[4，5]。但已有的商业软件集成的大都是理想弹塑性等经典弹塑性本构模型[6]，这类本构模型在屈服面内一般认为是弹性的，即不发生塑性累积变形的，并且往复循环加载时也无法计算塑性变形，因此这类型模型无法反映土体卸载体缩引起的孔压增加以及循环往复加载引起的孔隙水压力增加（剪缩引起的正孔压）或减小（剪胀引起的负孔压）[6]，导致这类模型难以合理的描述地基变形和孔压发展规律，不利于准确评价深厚覆盖层上高土石坝的抗震安全，因此必须开发适用于深厚覆盖层上高土石坝抗震安全评价的分析软件和本构模型。

本文基于大连理工大学自主开发的岩土工程高性能非线性有限元分析软件GEODYNA，考虑饱和多孔介质Biot动力固结理论，采用大连理工大学静动统一实用化弹塑性本构模型[9]，对地震作用下西部某水电站覆盖层地基加速度和孔隙水压力进行了数值模拟，验证了开发的本构模型和软件的合理性和可靠性。研究成果可为评价深厚覆盖层上高土石坝抗震安全提供可靠的评价技术。

1 离心机振动台试验

1.1 模型简介

试验以西部某水电站工程为研究对象，试验研究采用浙江大学研制的土工离心机[10]。试验模型箱采用叠环式，其内部尺寸为：740 mm×340 mm×425 mm（长度 × 宽度 × 高度），试验时离心加速度为50g。由于覆盖层地基断面材料比较复杂，为了便于离心机试验研究，对地基材料分区进行了简化，简化后的地基共2层（见图1）。试验上层砾石层选用配置的福建标准砂，厚度为100 mm，试验下层砂土选用现场砂，厚度为300mm。模型内部均设置了9个水平向加速度传感器（编号A1-A9）和9个孔压计（编号为P1-P9），加速度传感器和孔压计的位置相同（见图1）。

图1 模型试验情况Figure 1 Illustration of model test information

离心机振动台试验输入地震是硬梁包水电站场地谱地震波，输入地震波仅为水平向。如图2所示为A0测点实测的输入加速度时程（模型试验中实测值为其50倍）。在后文数值模拟中，以实测地震波作为计算输入地震动。

图2 模型试验中实测水平地震波Figure 2 Time history curve of seismic acceleration

1.2 计算模型简介

如图3所示是本次计算采用的计算网格图（基于模型实际尺寸建模），单元采用四边形等参单元。在后文进行动力试验时，地震动采用底部输入，为了模拟叠环式模型箱，对地基两侧节点的位移自由度绑定进行模拟分析。试验过程中地基土均为饱和土，为了模拟孔压的时程变化过程，因此在后文进行动力分析时采用基于动力固结方程和弹塑性模型的有效应力方法进行数值模拟。

图3 模型试验有限元网格（单位：mm）Figure 3 Finite element model （Unit：mm）

2 数值模拟方法与理论

2.1 动力流固耦合分析方法简介

基于Biot动力固结方程的有效应力分析方法是分析饱和土液化变形问题的有效途径。将Biot动力固结方程进行空间域离散，并引入系统阻尼并写成矩阵形式[11，12]，可得：

式中：M——土体的质量矩阵；

C——阻尼矩阵；

Q——耦合矩阵；

f(1)——土体的荷载向量；

Mf——流体的质量矩阵；

S——流体的压缩矩阵；

H——流体的渗透矩阵；

f(2)——流体的荷载向量。

本文在计算阻尼矩阵时阻尼比取0.05。大连理工大学在饱和多孔介质弹塑性动力框架上，将岩土工程有限元分析过程进行了类的抽象和封装，集成了土工构筑物和地基的动力有效应力分析方法，自主开发了岩土工程高性能非线性有限元分析软件GEODYNA。

2.2 地基土材料本构简介

土体本构模型的合理性直接影响数值模拟结果的可靠性。大连理工大学在压力相关和状态相关静动统一本构模型的基础上[13，14]，针对土体关键力学特性—应力历史相关性模拟存在的问题对塑性模量以及本构模型实用性等进行了进一步的发展[13，14]。改进模型基于比例记忆理论实现了不规则加载过程中重要应力历史的记忆，更合理的描述了复杂循环加载条件下土体应力应变关系。该模型参数的物理意义大都明确，且可根据常规试验直接确定。作者等（2018）[9]采用已有复杂应力路径条件下（包括，等小主应力、等平均主应力、等大主应力、等应力比条件下的单向循环和双向循环）的试验结果对模型进行了验证，结果表明，改进模型可以较好地反映不规则加载条件下土体的循环滞回、循环硬化等变形特性，尤其是可以考虑复杂循环不规则加载历史对后期加载变形的影响。同时基于第十五届国际大坝数值分析标准研讨会发布的Menta面板堆石坝的动力反应分析案例，验证了本构模型在模拟筑坝堆石料力学特性方面的合理性[15]。提出模型同样适用于砂土，表1给出了提出模型的本构模型参数，其是根据离心机模型试验中现场砂和福建标准砂土的三轴试验成果标定获得。

表1 本构模型参数Table 1 Model parameters of soil

3 计算与试验结果对比分析

3.1 地震加速度

以下数值模拟和试验结果，均是根据缩尺规则将模型结果等效到原型。图4给出了A2（距离模型箱底部的高度为0.15m）、A5（距离模型箱底部的高度为0.25m）、A8（距离模型箱底部的高度为0.35m）三个测点处的加速度时程与实测的对比。从图4中可以看出：埋深最深测点吻合的最好，埋深最浅的测点吻合较差；位于下侧现场砂土层中的A2、A5测点，计算和实测加速度幅值大小变化规律基本吻合；位于上侧福建标准砂土层的A8测点，在约6s以前计算和实测加速度幅值大小基本一致，在6s以后计算加速度幅值明显小于实测值。

图4 计算和实测加速度时程对比Figure 4 Comparison between measured and simulated acceleration time histories

图5给出了A0、A8点计算和实测傅里叶谱的结果，可以看到：实际输入地震动的主频约为1.6Hz，其与模型基频（震前约2.0Hz）存在一定的差异，地震过程中因孔压上升导致地基模量降低导致模型基频略有降低，震后的基频约为1.7Hz。在频率1.7Hz位置（与地震动主频接近），计算和实测的傅里叶谱幅值均出现了峰值，但两者傅里叶谱值相差不大。在1.0Hz附近区域，计算的傅里叶谱极值大于实测值，这与计算和实测的加速度时程峰值规律并不相同。在大于3.0Hz的范围内，计算的傅里叶幅值要明显的小于实测值。综上可知，计算加速度的峰值较小是因为高频含量较小导致的。造成这种差异的原因可能一方面与数值分析中未考虑模型箱底部和侧边界地震波反射有关[16]；另一方面与数值分析模拟接近液化低有效应力时土体阻尼比大小选择有关；此外有学者认为这种现象与土体的剪胀有关，但从后文孔压时程可以看到，计算孔压和实测值吻合较好[17]。限于篇幅，本文并未对上述因素进行数值分析探讨。

图5 实测和计算傅里叶谱（单位：g）Figure 5 Comparison between measured and simulated acceleration Fourier spectrum （Unit：g ）

3.2 孔隙水压力

图6 给出了P2、P5和P8测点的计算和试验孔压发展过程对比。从图6中可以看出：位于下侧现场砂土层的P2测点计算和实测孔压值发展过程和量值均比较接近；位于下侧现场砂土层的P5测点实测孔压略小于计算值，但两者的变化规律是一致的；位于上侧福建标准砂土层的P8测点实测孔压略大于计算值，但孔压累积和消散的过程是基本一致的。总体来讲，实测与计算的孔压无论在量值和发展规律上都吻合较好。

图6 孔压时程对比Figure 6 Comparison between measured and simulated pore water pressure time histories

4 结论

深厚覆盖层上建设高土石坝工程是坝工界正在面临的难题之一，其中强震作用下覆盖层地基变形和稳定是深厚覆盖层上高土石坝的抗震安全评价重要内容。本文应用大连理工大学自主研发的岩土工程高性能非线性有限元分析软件GEODYNA和静动统一实用化弹塑性本构模型，对西部某水电站覆盖层地基离心机试验进行了动力弹塑性有效应力数值分析。结果表明：计算和实测加速度幅值大小变化规律吻合较好；计算的孔隙水压力随地震时间的上升斜率、幅度大小变化及消散趋势与试验数据基本一致，说明本文采用的本构模型能合理描述地震作用下地基土的孔压发展规律。

本文研究成果验证了开发的本构模型和软件的合理性和可靠性，在此基础上联合已开发的土—界面—结构体系跨尺度计算理论和方法[18-20]，可为实际深厚覆盖层高土石坝地震安全评价和抗震措施设计提供精细化数值分析工具。