赵杰　徐剑　王桂萱　尹训强

摘要：从岩石材料细观非均匀性特点出发，运用连续介质损伤力学理论，推导了岩石材料细观损伤本构关系式。以ANSYS二次开发平台为主体，基于细观损伤力学，提出基于宏观均质假定的等效损伤单元来模拟隧道洞口段岩质边坡的损伤开裂过程。计算结果表明：在细观方面边坡损伤破坏是个不断损伤的过程，在地震荷载作用下，边坡结构首先会出现微裂纹，随着地震动时间增加，已出现的微裂纹逐渐扩展形成可见的裂缝，最后形成连续的损伤区域；尤其在隧道拱顶上方形成连续贯通的损伤区域，可看出拱顶已完全破坏，给隧道安全运行带来不利影响，因此在设计施工时要对拱顶上方进行加固处理。该计算方法和分析成果可用于同类复杂结构的损伤破坏分析。 关键词：岩质边坡；非均质性；损伤力学；等效损伤单元；损伤破坏；ANSYS

中图分类号：U45 文献标识码：A 文章编号：1000-0666（2016）01-0053-07

0 引言

近年来发生的重大岩石动力灾害事件表明：矿山巷道、地下洞室、隧道工程、边坡工程垮塌和破坏，多数是由于岩体内的节理、裂隙及其扩展造成（罗国煌等，2004；李明田等，2003）。岩体内的节理、裂隙就是其损伤的实质性表现。所以，从细观损伤理论出发，研究围岩体细观损伤扩展机理和演化特征对岩体工程稳定性的评价与支护设计是非常有意义的。

Kachanov在1958年提出损伤力学思想，至今已经过50多年的发展，一些理论已日益成熟并在各学科和工程技术领域得以应用。周维垣等（1998）应用连续介质损伤力学理论，从岩体内部微裂纹产生和扩展的损伤机理出发，建立了坚硬岩体的弹性一损伤耦合的各向异性弹脆性损伤本构模型。朱万成等（2003）从岩石的细观非均匀性特点出发，应用弹性损伤力学对岩石材料细观的损伤演化过程进行了描述。赵宝友等（2003）基于损伤力学原理，通过引入损伤变量的方法详细推导了ABAQUS中的损伤塑性本构模型屈服阶段及其卸载后的损伤机制和力学行为，并以某水电站地下厂房洞室结构为例，进行二维动力时程非线性分析，得出强震下采用此损伤塑性模型是合理和必要的。隋斌等（2008）利用FLAC3D软件模拟了地震作用下某地下洞室群围岩的动态响应，采用新的劈裂破坏判据对震后可能出现的劈裂损伤区范围进行了预测。左双英和肖明（2009）基于深埋大型地下洞室群围岩损伤机理，对映秀湾水电站洞室群在三向地震荷载耦合作用下的动力响应进行数值模拟，得出地震作用下围岩损伤分布的特征。张志国等（2010）采用动力时程分析法研究地下洞室地震响应特性，结合三维弹塑性损伤有限元理论，编写了动力时程计算程序，分析了渔子溪水电站地下厂房洞室群的围岩稳定。

笔者在上述研究的基础上，基于岩石材料的宏观均质假定，运用连续介质弹性损伤力学模型来描述单元的材料特性，采用最大拉应力与摩尔一库伦破坏准则作为损伤阈值，以通用有限元软件ANSYS开发平台为主体，运用UPFs的二次开发特点，基于用户自定义单元的功能将损伤单元法引入到ANSYS中，用来分析地震作用下隧道洞口段岩质边坡从微裂纹产生、扩展到形成可见裂缝，直至最终的损伤过程。

1 损伤演化本构模型

1.1 岩石材料非均质性特征

岩石作为一种自然地质体，由于各种因素使岩石内部存在大量不同的微缺陷，所以岩石介质的构成非常复杂，通常对其进行数学描述也是非常困难的。唐春安（1993）认为可以利用统计的方法对离散后的岩石介质进行近似描述。建立材料细观非均质特性随机模型时，所谓非均质性是指材料的力学参数，如强度、弹性模量、泊松比等。在给定的空间上对模型进行网格划分，为了描述材料参数的这种空间变化，需要解决两个问题：第一，如何使材料参数在数值上满足非均匀分布规律；第二，如何将这些不同材料参数随机地分布到划分的单元中去。本文将岩质边坡看成一个大的样本空间，采用等效损伤单元对岩石材料进行离散，使得离散后的单元形式具有类似于普通有限单元的特性，并且认为其本身为均匀、连续的介质，然后在样本空间内引入概率统计的方法，假定各单元的材料力学性质存在差异且服从Weibull概率统计分布，从而可表现出一定的非均匀型和离散性。杨天鸿等（2004）认为采用Weibull随机分布可以较好地模拟材料性质的非均匀性，该分布可以用分布密度函数来定义：式中，μ为岩石材料随机分布力学性质参数（强度、弹性模量、泊松比等）；μ⁰为岩石材料基元体力学性质平均值；m为岩石力学性质分布函数的密度形状函数，其物理意义反应了岩石材料的均质性；f（μ）是岩石力学性质μ的统计分布函数。

式（1）很好地反映了岩石材料的细观力学性质非均匀性分布情况。随着均匀性系数m的增加，基元体的力学性质将集中于一个狭窄的范围之内，表明单元强度趋于均匀；当m值减小时，基元体的力学性质分布范围变宽，表明岩石介质的性质趋于均匀。可见m值反映了统计模型中材料结构的离散程度，能够较为真实地反映岩石材料的非均质性。

1.2 岩体本构关系与屈服准则

余大庆（1993）认为岩石应力-应变曲线是由于其受力后的不断损伤引起微裂纹萌生和扩展而造成的，而不是由于其塑性变形。尤其是在拉伸应力作用下，其脆性更加明显。因此，用弹性损伤力学本构模型来描述岩石材料细观单元力学性质是合适的。在本文的研究方法中，等效单元应力一应变本构关系亦采用了弹性损伤模型。在施加地震荷载前，用弹性模型和泊松比来定义每一个等效损伤单元，并且将其假定为线弹性、均值、无损伤的。当单元进入损伤阶段后，依据Le-maitre应变等价原理，认为应力作用在受损材料上引起的应变与有效应力作用在无损材料引起的应变等价。据这一原理，受损材料单元的本构关系可通过无损材料中的名义应力得到：式中，E⁰为材料初始弹性模量；E为材料损伤后的弹性模量，D为损伤变量。D=0表示材料处于无损伤状态，D=1表示处于完全损伤（断裂或者破坏）状态，00

在岩石材料损伤演化中，随着等效损伤单元应力会不断增加，当单元的应力或应变状态满足某个给定的损伤阈值（准则）时，单元开始损伤。这里选择两个损伤阈值准则，第一个采用最大拉应变准则，当细观单元的最大拉伸主应变达到其给定的应变阀值时，该单元开始发生拉伸损伤；第二个采用摩尔一库仑准则，当细观单元的应力状态满足摩尔一库仑准则时，该单元发生剪切损伤。需要说明的是，岩石材料在外部荷载作用下，其抗拉强度远低于抗压强度的特点，且单元的拉伸损伤和剪切损伤不会同时发生于同一荷载步下。换句话说拉伸准则具有优先权，若细观单元满足拉伸准则，则不需再判断其是否满足摩尔一库仑准则。只有在损伤单元未满足拉伸准则之后才判定其是否满足摩尔一库仑则。图1、2分别给出了单元材料拉伸本构关系曲线与压缩本构关系曲线图。单元满足拉损伤模型以及摩尔一库仑准则对应的损伤变量表达为式中，f_tr为单元参与强度；ε_t0为拉损伤应变阈值；ε_tu为单元材料极限拉伸应变；ε_c0为单元应力达到单元材料单轴抗压强度时的主要应变；ε_cu为单元材料的极限应变；ε_c为单元的单轴压缩应变；N为与下降段形状有关的参数。

2 等效损伤单元在ANSYS中的实现

本文等效损伤单元采用的本构关系与屈服准则是其与普通有限单元最大的区别之处。基于等效单元特殊性，需要以ANSYS为开发平台，利用UPFs为开发工具进行二次开发。图3给出了基于ANSYS为二次开发平台的等效损伤单元建立流程，其中虚线方框部分由UPFs的功能实现。详细的程序流程如下：

首先，建立隧道洞口段岩质边坡有限元模型，对用户进行单元编译和链接，施加位移边界条件，控制求解流程以及结构的输出（此部分利用GUI和APDL完成）；然后，进入地震动力分析求解阶段，在每一荷载步施加地震荷载向量；最后，在ANSYS数据库中将建立单元的所有信息输入到接口子程序UserElem.f中，并完成以下步骤（此部分利用子程序接口UserElem.f完成）：

①检验是否第一次进入等效损伤单元子程序，若“是”，则利用蒙特卡罗方法生成服从Weibull分布的材料参数随机数，然后生成最终的单元材料力学参数，将此数据存储在数组中并使用FOR-TRAN语言中的COMMON将其定义为全局变量，以便以后的平衡迭代步和荷载步调用，进入②；若“否”，则将直接获得上一迭代步的材料力学参数，进入②。

②通过节点真实位移{u}_t+{△u}^it计算单元的应变s；由损伤变量Dⁱ_l（Dⁱ_lt表示拉损伤；Dⁱ_lc表示压损伤，其中i表示当前迭代步，l表示当前荷载步）判断该单元是否发生过损伤，即Dⁱ_l是否大于D⁰_l，其中，0是上一荷载步的最后迭代步，若“否”，则说明该单元此前未发生损伤，进入③；若“是”，则说明该单元此前发生过损伤，进入④。

③由②所求应变s检验是否发生拉损伤，如发生了则进入⑤，若未发生拉损伤，则检验是否为压损伤，如发生了则进入⑥，未发生ANSYS数据库中，计算整体结构响应。

④判断发生损伤的类型：若为拉损伤，计算压损伤变量Dⁱ_lc。

⑤计算相应的拉损伤变量Dⁱ_lt以及弱化后弹性模量E。

⑥计算相应的压损伤变量Dⁱ_lc以及弱化后弹性模量E。

⑦计算单元矩阵，包括单元刚度阵、阻尼阵、质量阵以及单元等效荷载向量，然后将上面得到数据返回到ANSYS数据库中，计算整体结构响应。如果计算结果不满足收敛条件，则进入下一平衡迭代步，若满足则进入下一荷载步，直至地震分析结束。

3 岩质边坡损伤演化模拟

3.1 边坡计算参数及有限元模型

计算参数选取：弹性模量E=3000GPa，泊松比v=0.23，质量密度ρ=2300kg/m³，动抗拉强度f_t=1.5×10⁵。选取材料的随机力学参数为：弹性模量和抗拉强度的均质度m=2，泊松比的均质度m=100。选取瑞利阻尼进行计算，计算时取第1阶与第10阶的固有频率，阵型阻尼比为0.05。本文所建立的边坡模型坡度为45°，坡高为20m，坡前计算长度取20m，坡后计算长度取44m，坡脚以下取20m。为研究隧道洞口段岩质边坡在地震作用下的损伤破坏过程，在平面应变状态下，采用等效损伤单元（几何形状如四节点平面应变单元）对边坡模型进行有限元网格离散，离散后的单元尺寸约为0.4m×0.4m，模型共分18589个单元、18941个节点。模型左右两边施加水平约束，底部施加水平和竖向约束，边坡有限元模型见图4。分水平与竖直2个方向输入地震波，其中水平峰值加速度为0.474g，竖向加速度峰值为0.312g，地震动总的持续时间为10.56s，加速度时程曲线见图5。

3.2 边坡地震损伤演化分析

本文从地震荷载作用整个过程中，选取6个典型的时间点（2.64s、2.92s、5.78s、6.96s、8.98s、9.42s）来分析隧道洞口段岩质边坡从开始出现微裂纹到形成清晰可见的裂缝，裂缝不断扩展至形成最终损伤区域总过程，其损伤演化过程及其主拉应力的分布见图6。隧道洞口段岩质边坡损伤演化分析过程分析中：在静力荷载（本文仅施加了重力荷载）作用的基础上施加地震荷载。地震荷载激励的初期阶段，边坡绝大部分区域处于受压状态下，虽然结构没有出现严重的裂缝，但由于单元材料随机分布的影响，少数材料强度较弱的单元依旧出现了不同程度的损伤，坡面开始出现微裂纹，但边坡结构整体上处于完好的状态。随着地震激励时间的增加，已出现的微裂纹不断扩展，新的微裂纹也不断出现，这些微裂纹逐渐扩展、延伸，最后在坡后形成数条清晰可见的裂缝，如图6a、b所示；随着地震动的继续震荡，坡后裂缝范围继续扩展、延伸，坡面发生破坏区域逐渐扩大；洞口拱顶上方和坡面中间开始出现裂纹并扩展成较大裂缝，形成多处贯通的损伤区域，同时造成了主拉应力的重新分布，如图6c、d所示；最后，随着地震动的继续震荡，拱顶上方裂缝继续扩展，最终形成了贯穿整个拱顶上方的损伤区域，此时可以认为拱顶上方已完全破坏，会给隧道的安全运行带来影响；坡面中间形成的自上而下裂缝会继续扩展、延伸，最终形成了一条清晰可见的裂缝，可以近似将该裂缝看成边坡在强震作用下产生的滑坡趋势，如图6e、f所示。

4 结论

本文通过引入服从Weibull函数的细观损伤微元体，运用连续损伤介质理论，理论上推导出了可反映岩石非均质性和宏观破裂与失稳过程的岩石细观损伤本构关系式。

地震作用下岩质边坡的破坏是材料细观损伤不断产生、累积的结果。在地震荷载激励的初期阶段，边坡绝大部分区域处于受压状态下，虽然结构没有出现严重的裂缝，但由于单元材料随机分布的影响，少数材料强度较弱的单元依旧出现了不同程度的损伤，坡后开始出现微裂纹，这些微裂纹逐渐贯通并形成数条清晰可见的裂缝。随着地震动的继续震荡，在洞口拱顶上方和坡面中间开始出现新的裂缝。拱顶上方裂缝最终形成一个连续贯通的损伤区域，此时可以认为拱顶上方已完全破坏，这会给隧道的安全运行带来影响，因此在洞口设计、施工时要对拱顶上方进行必要的加固处理；坡面中间形成的自上而下的裂缝，可以看成是边坡结构在强震作用下产生的滑坡趋势，设计、施工也要给与重视。