基于CEM的加锚节理岩体动力学算法研究*

23

(1.长江勘测规划设计研究有限责任公司，湖北 武汉 430010；2.武汉大学土木建筑工程学院，湖北 武汉 430072；3.陕西省混凝土结构安全与耐久性重点实验室，陕西 西安 710123；4.长江岩土工程总公司(武汉)，湖北 武汉 430010)

大量岩土边坡地震崩塌与滑动主要发生在多震的西部地区，这些地区地形地质条件非常复杂，岩体中存在断层、节理、裂隙等各种结构面。这些地区复杂的地质构造和较高的地震烈度使得锚固岩体的抗震稳定也成为工程的关键问题之一。

数值计算方法[1]能够体现出复杂岩体的力学特性，已被成功地应用于岩质边坡的静力稳定性分析中，在岩质边坡的动力计算研究方面也取得了较大的进展，众多学者采用有限单元法(FEM)[2]、非连续变形分析法(DDA)[3]、离散单元法[4]等多种计算方法对节理岩体进行动力学计算分析，但在模拟加锚节理岩体方面大多处于初步发展阶段，考虑的形式比较简单，还需诸多努力。目前仍然缺少一种既能同时实现对结构面、加固措施的精细模拟，又能满足动态设计对于建模效率要求的数值模拟方法[5]。

为解决这个难题，相关专家提出了复合单元法(CEM)，该方法主要有以下三个优点：ⓐ准确：能够实现结构面和加固措施的离散模拟，具有细致精确的优点；ⓑ快速：在建模时无须考虑结构面和加固措施，通过引入复合单元，用复杂的计算问题来代替复杂的结构面和加固措施建模问题，拓扑计算得到整合后的复合单元均由程序自动实现；ⓒ动态设计：当设计方案优化调整时，不需要重新建模，根据新的支护措施，程序重新拓扑计算得到新的复合单元即可。

目前为止，静力学计算方面，复合单元法已经实现了对结构面(如节理、断层等)[5]、多物理场模拟(如温度、渗流)等和加固措施(如锚杆、锚索等)的模拟，但是，在动力学计算方面，仅实现了节理岩体模拟[6]，还未实现对加锚节理岩体动力学计算。

本文进一步发展了复合单元法，建立了加锚节理岩体复合单元法动力学算法，并首次将粘弹性边界在复合单元上进行了初步应用，通过算例证明了算法的合理性，为加锚节理岩体的动力学计算研究提供了新的思路。

1 复合单元法动力学算法步骤

复合单元法动力学模型的建立和有限元动力分析方法是类似的，通过位移空间离散化，可以最终得到系统求解方程为：

(1)

M——质量矩阵；

C——阻尼矩阵；

K——刚度矩阵；

P——节点载荷向量。

在得到上述系统运动方程之后，可以选择用直接积分法或者模态分解法来求解系统响应。

复合单元法动力学方法可以分为以下几步：

a.生成复合单元网格信息，一个复合单元可能会包括若干子单元和接触面。

d.根据b和c生成的复合单元刚度矩阵和复合单元质量矩阵进行模态分析，根据其自振频率和阻尼比组装形成Rayleigh阻尼矩阵Ce。

e.将复合单元区域的积分结果和有限元积分结果求和，即组装系统整体刚度矩阵K=∑Ke、整体质量矩阵M=∑Me、整体阻尼矩阵C=∑Ce和整体节点载荷向量P=∑Pe。

f.求解系统方程。

2 复合单元质量矩阵和积分方法

对复合单元采用一致质量矩阵为：

(2)

在复合单元法中，对于含各种不连续面的复合单元，对于子单元块体积分方法是先将组成该块体的顶点按照一定顺序和规则连接形成若干四面体，再将每个四面体用行波法自动剖分网格的方式分成4个六面体，具体方法见图1。

图1 子单元块体积分方法示意

在复合单元法中，对于锚杆杆体和砂浆子单元块体积分方法为是建立局部笛卡儿坐标系和局部柱坐标系，对其进行积分，坐标系见图2。

图2 局部笛卡尔坐标系和锚杆局部柱坐标系

3 复合单元质量刚度矩阵

加锚节理岩体的复合单元示意见图3。

图3 加锚节理岩体的复合单元示意

如图3所示，在复合单元内，分别定义不同材料区域为子单元，在复合单元引入5套虚节点加上实节点共6套节点，用于描述两个块体中的岩石、锚杆、砂浆等子单元，相应地定义6套位移分量分别用来插值6个子单元内部的位移。

根据虚功原理，可得

(3)

经过推导得到以下平衡方程组：

(4)

(5)

4 复合单元阻尼矩阵

阻尼矩阵采用Rayleigh阻尼假定，其阻尼矩阵是质量矩阵和刚度矩阵的组合：

[C]=α[M]+β[K]

(6)

式中α——Rayleigh阻尼的质量系数；

β——Rayleigh阻尼的刚度系数。

将上述复合单元的矩阵[M]、[C]和[K]代入式(1)，然后采用Wilson-θ法进行每一个时间步的求解。

5 复合单元黏弹性人工边界

进行边坡动力分析需要考虑边坡与周围基岩的动力相互作用，无限地基如何模拟是常常要面对的问题。黏弹性人工边界是一种十分理想的动力边界[7]，基本思想是在人工边界节点的法向和切向设置并联的弹簧和阻尼器来实现。

三维人工边界上法向和切向的弹簧刚度和阻尼系数取值如下：

(7)

(8)

式中kbN——单位面积的法向弹簧刚度；

kbT——单位面积的切向弹簧刚度；

CbN——单位面积的法向阻尼器的阻尼系数；

CbT——单位面积的切向阻尼器的阻尼系数；

r——波源至人工边界点的距离；

λ、G——介质Lame常数；

ρ——介质质量密度；

cP——介质的P波波速；

cS——介质的S波波速；

A、B——参数，取建议值[7]为A=0.8，B=1.1。

当黏弹性边界单元含有复合单元时，以复合单元内子单元为基本单位，依次得到每个子单元对应的虚实节点的法向与切向弹簧刚度和法向与切向阻尼器的阻尼系数，取值为单位面积的法向与切向弹簧刚度和单位面积的法向与切向阻尼器的阻尼系数均分别乘以其有效面积。

6 技术路线

技术路线见图4。

图4 技术路线

7 算例考证

7.1 算例一

取一个尺寸为1m×1m×10m(长×宽×高)的岩体，岩体中心内含有一条长度为10m砂浆锚杆，钢筋的半径16mm，砂浆的半径80mm，5m高处含有一个1m×1m水平节理。复合单元法计算网格见图5，单元数为480，节点数为820。模型底部全约束。表1为实体材料力学参数。表2为接触面力学参数。模型Y方向施加线荷载F(t)=100kN/m(0≤t≤1)，作用位置见图5。采用一阶振型的瑞利阻尼，分别考虑无阻尼、阻尼比0.05和阻尼比0.01的计算工况。

图5 复合单元法网格和荷载示意

表1 弹性材料参数

表2 接触面材料参数

7.2 算例一结果分析

图6～图8是以加锚节理岩体为对象不同阻尼比下线荷载作用中点位移历时曲线图，算例对比了复合单元法和有限单元法的计算结果，从图中可以看出，随着阻尼系数的增大，达到平衡位置的时间变短，并且，复合单元法和有限单元法这两种方法计算的结果变化规律相似，误差较小，量值也较接近，说明复合单元法的计算结果是合理的。

图6 无阻尼时线荷载作用中点位移历时曲线

图7 阻尼比0.01线荷载作用中点位移历时曲线

图8 阻尼比0.05线荷载作用中点位移历时曲线

7.3 算例二

从半无限空间截取一个尺寸为4m×1m×3m(长×宽×高)的弹性材料作为有限计算区域，3m高处还有一水平节理，顶端自由，法向为Y的面法向约束，法向为X和Z的面施加粘弹性边界。材料弹性模量为2.5Pa，泊松比为0.25，材料密度为1kg/m3，节理法向刚度为10Pa/m，切向刚度为5Pa/m。时间步长为0.01s，总计算时间为15s。动荷载为作用在顶面的压力荷载F(t)，数值大小见图9。网格模型和测点见图10，单元数为300，节点数为672。

图9 自由表面的压力荷载F(t)历时曲线

图10 网格模型和测点示意

为了研究半无限空间粘弹性边界计算结果的合理性，增加两组计算工况。工况一：采用固定边界，保持模型和荷载不变，将粘弹性边界改为法向约束；工况二：在工况一的基础上，将模型尺寸变为400m×1m×300m(长×宽×高)，顶面自由，其余为法向约束，此工况的结果可以作为准理论解。

7.4 算例二结果分析

图11和图12为测点1和测点2不同工况下的Z向位移历时曲线，从计算结果可以看出，粘弹性人工边界能够很有效地吸收散射波在边界处产生的能量，消除边界的反射，与准理论解从量值和规律上都比较吻合。固定边界由于不能吸收边界处散射波的反射，夸大了能量的输入，造成了振荡现象。由此，也可以证明粘弹性边界在复合单元法中的应用是合理的。

图11 测点1不同工况下的Z向位移历时曲线

图12 测点2不同工况下的Z向位移历时曲线

8 结论与展望

本文开展了加锚节理岩体复合单元法动力学研究，推导了加锚节理岩体复合单元法质量矩阵、刚度矩阵和阻尼矩阵，建立了复合单元粘弹性人工边界，建立了惯性力、阻尼力、动力荷载及弹性力作用下的复合单元动力系统的控制方程组，编制相应程序，并进行算例验证。

a.通过算例结果的对比可知，复合单元法的计算精度与有限单元法基本一致，粘弹性边界在复合单元法中的应用是合理的，验证了复合单元法能有效地模拟加锚节理岩体，进行动力学计算。

b.提供了一套完整的锚固岩体复合单元法动力学分析的技术路线，为加锚节理岩体动力学计算提供了一套高效的数值模拟方法。

c.本文建立的加锚节理岩体复合单元法动力学模型仅是初步研究成果，只能处理一些简单的算例，仅考虑了弹性状态，后续研究工作需要合理的锚固岩体动力学本构模型作为技术支撑。

d.复合单元法进行动力分析时，复合单元粘弹性边界条件下地震波的输入，垂直从地基底部入射或者考虑波动斜入射等问题，都值得未来更加深入地研究。