◎陈斯祺 深圳市深水水务咨询有限公司

滑坡是一种全球性的地质灾害之一，具有随机性、突发性和瞬时性等特征。岩质边坡的断裂破坏很大程度上是由于受到力的作用后，产生了变形和裂纹。由于裂纹的存在，在荷载的作用下，裂纹会不断地扩展，最终导致滑坡断裂破坏。通过有限元程序ANSYS求解材料应力强度因子，并将计算结果和理论分析结果进行对比分析，检验得出强度因子的精度，以此来判断利用ANSYS来求解材料的应力强度因子的准确性和可行性。

对于滑坡断裂破坏机制研究，有利于从根本了解滑坡的断裂破坏过程及机理，从而减少滑坡灾害意外事故的发生，降低滑坡灾害造成的经济损失和保障人民的生命安全，对滑坡的预测和治理具有重大意义。

1.研究现状

随着断裂力学在工程应用中的高速发展，断裂力学理论已经逐渐渗透到了工程各个领域。Price(1966)首次将断裂力学理论应用到岩石材料中，并提出了岩石裂隙的形成是源于表面能增加的结果；陈洪凯针对危岩体的稳定性，研究了其断裂主控结构面的强度参数计算方法，对危岩演化力学机理建立了比较系统的理论，并从断裂力学角度分析了危岩体的稳定性和破坏机理；赵海涛对应力强度因子在线弹性断裂力学以及ANSYS中的计算方法进行了分析，通过算例表明，用ANSYS求解材料的断裂强度因子是准确的、可行的。

基于ANSYS有限元软件对断裂力学分析的研究成果较多，但只是单纯利用ANSYS有限元软件对断裂力学进行分析，并没有考虑软件对断裂强度因子的分析精度。实际工程中的许多裂纹，其走向通常与边坡坡面大体一致，因此可将其作为二维平面问题来研究。经过断裂理论长期的发展，二维问题的研究相比三维裂纹问题要较为成熟，本文将其作为二维平面问题分别对Ⅰ型裂纹和Ⅱ型裂纹进行数值模拟，通过从ANSYS有限元软件中提取的Ⅰ型和Ⅱ型应力强度因子分别与其理论计算对应得到的Ⅰ型和Ⅱ型应力强度因子进行对比，检验从ANSYS有限元软件中提取的Ⅰ型和Ⅱ型应力强度因子的精度。

2.Ⅰ型裂纹数值模拟

薄板裂纹几何尺寸如下：正方形其边长L=1 m，裂纹长度2 a，板厚W=0.0001 m，Ⅰ型薄板中心裂纹的数值模拟采用2 维8 节点高阶PLANE183单元进行模拟。Ⅰ型薄板中心裂纹的材料参数：弹性模量为2.1MPa，泊松比为0.35。荷载为：双向均匀拉应力σ。

该数值模拟主要是为了分析在不同裂纹尺寸及不同应力状态下对Ⅰ型裂纹从A NS YS 有限元软件中提取出来的应力强度因子的精度影响。裂纹长度分别设置为0.04m、0.08m、0.12m、0.16m和0.20m，双向均匀拉应力σ分别设置为2MPa、4MPa、6MPa、8MPa和10MPa。共计25组数值模拟工况，分别进行如表1编号。

表1 Ⅰ型不同数值模拟工况的编号

首先，在ANSYS有限元软件的工作界面中建立Ⅰ型薄板中心裂纹模型，再将建立好的裂纹薄板模型图进行网格划分，将裂纹尖端处进行网格加密处理，有利于提高从ANSYS有限元软件中提取的应力强度因子的精度（见图1）。最后，通过ANSYS有限元软件将已建立好的裂纹薄板模型设置求解类型并定义裂纹薄板模型的荷载，进行后处理，形成裂纹尖端等效应力云图（见图2）。

图1 网格划分图

图2 Ⅰ型裂纹尖端等效应力云图

围绕裂纹尖端定义三组节点，且其中每组节点必须相邻。通过ANSYS有限元软件中的相互作用积分法计算出裂纹尖端应力强度因子。为了减小误差，取三组的平均值作为数值模拟的Ⅰ型应力强度因子与理论结果进行对比。

当裂纹长度为0.04m、0.08m、0.12m、0.16m和0.20m时，Ⅰ型裂纹薄板的双向均匀拉应力σ分别设置为：2MPa、4MPa、6MPa、8MPa和10MPa。可分别得表2至表6所示数据。

表2 A组-裂纹长度2a=0.04m

表3 B组-裂纹长度2a=0.08m

表4 C组-裂纹长度2a=0.12m

表5 D组-裂纹长度2a=0.16m

表6 E组-裂纹长度2a=0.20m

通过数值模拟的平均Ⅰ型应力强度因子与理论公式计算得到的Ⅰ型应力强度因子，可得到如图3所示关系图。

图3 理论值与模拟平均值关系图

由关系图可知，当裂纹长度分别为0.04m、0.08m、0.12m、0.16m和0.20 m时，Ⅰ型裂纹薄板的双向均匀拉应力σ分别设置为：2MPa、4MPa、6MPa、8MPa和10MPa。数值模拟的平均Ⅰ型应力强度因子分子与理论公式计算得到的Ⅰ型应力强度因子十分接近，相对误差结果分别均为0.9%、2.0%、3.5%和5.5%。由此可证明数值模拟值产生的误差大小与其应力状态无关，仅与数值模拟网格划分有关。

3.Ⅱ型裂纹数值模拟

薄板裂纹几何尺寸如下：正方形其边长L=1 m，裂纹长度2 a，板厚W=0.0001m，Ⅱ型薄板中心裂纹的数值模拟采用2 维8 节点高阶PLANE183单元进行模拟。Ⅱ型薄板中心裂纹的材料参数为：弹性模量为2.1MPa，泊松比为0.35。荷载为：双向均匀剪应力τ。

该数值模拟主要是为了分析在不同裂纹尺寸及不同应力状态下对Ⅱ型裂纹从A NS YS 有限元软件中提取出来的应力强度因子的精度影响。裂纹长度分别设置为0.04m、0.08m、0.12m、0.16m和0.20m，双向均匀剪应力τ分别设置为2MPa、4MPa、6MPa、8MPa和10MPa。共计25组数值模拟工况，分别进行编号（见表7）。

表7 Ⅱ型不同数值模拟工况的编号

首先，我们在ANSYS有限元软件的工作界面中建立Ⅱ型薄板中心裂纹模型。将建立好的裂纹薄板模型图进行网格划分，将裂纹尖端处进行网格加密处理，有利于提高从ANSYS有限元软件中提取的应力强度因子的精度（见图4）。最后，将已建立好的裂纹薄板模型设置求解类型并定义裂纹薄板模型的荷载，进行后处理形成等效应力云图（见图5）。

图4 网格划分图

图5 Ⅱ型裂纹尖端等效应力云图

当裂纹长度分别为0.04 m、0.08m、0.12m、0.16m和0.20m时，其双向均匀剪应力τ 分别设置为：2MPa、4MPa、6MPa、8MPa和10MPa。可得表8至表12所示数据。

表8 F组-裂纹长度2a=0.04m

表9 G组-裂纹长度2a=0.08m

表10 H组-裂纹长度2a=0.12m

表11 I组-裂纹长度2a=0.16m

表12 J组-裂纹长度2a=0.20m

通过数值模拟的平均Ⅱ型应力强度因子与理论值，可得到如图6所示关系图。

图6 理论值与模拟平均值关系图

由关系图可知，当裂纹长度分别为0.04m、0.08m、0.12m、0.16m和0.20m时，其双向均匀剪应力τ分别设置为：2MPa、4MPa、6MPa、8MPa和10MPa。数值模拟的平均Ⅱ型应力强度因子与理论公式计算得到的Ⅱ型应力强度因子十分接近，相对误差结果分别均为0.1%、0.5%、1.1%、2%和3.2%。由此可证明数值模拟值产生的误差大小与其应力状态无关，仅与数值模拟网格划分有关。

4.结语

（1）通过数值模拟得到的平均Ⅰ型应力强度因子分子与理论公式计算得到的Ⅰ型应力强度因子十分接近，数值模拟精度较高，相对误差结果分别均为0.9%、2.0%、3.5%和5.5%，相对误差最大为5.5%。由此可证明数值模拟值产生的误差大小与其应力状态无关，仅与数值模拟网格划分有关。

（2）通过数值模拟得到平均Ⅱ型应力强度因子与理论公式计算得到的Ⅱ型应力强度因子十分接近，数值模拟精度较高，相对误差结果分别均为0.1%、0.5%、1.1%、2%和3.2%，相对误差最大为3.2%。由此可证明数值模拟值产生的误差大小与其应力状态无关，仅与数值模拟网格划分有关。