裂隙岩体的三维非线性有限元分析

摘要：讨论天然裂隙岩体的机械力学行为及力学问题，建立了相应的有限元格式，介绍了裂隙岩体刚度方程的建立及弹塑性状态下的应力——应变关系是非线性的表达式，阐述了非线性问题的分类及解法。

关键词：裂隙岩体；接触单元；非线性问题

工程岩体经受过复杂的地质作用后，其内部存在不同规模、性质各异的地质结构面（节理、裂隙、断层等），使得岩体的结构特征与力学性质错综复杂。一方面，不连续面的存在使得岩体介质类型不同于其他任何一种力学材料。在进行裂隙岩体力学问题分析时，准确选择岩体力学模型成为必须解决的首要问题。另一方面，岩体中的不连续面使得岩体性质表现出不连续、非均匀、各向异性及尺寸效应等特征。这就使得传统弹性有限元法不能更好地模拟真实情况。

近年来，节理单元被广泛应用于模拟岩体中的断层、剪切带和不连续区域。在节理单元中，由于引入了岩石燕子节理的法向和切向的刚度，使得节理单元的应用受到了一定的限制。这是因为刚度系数通常都是从野外实测得到的，测试存在困难和不确定性。Katona于1983年提出了一种不用刚度系数的简单接触摩擦单元，它能模拟两物体之间的滑动摩擦、张开和闭合过程，适用于不计厚度影响的节理、断层和位移不连续场的有限元分析。由于这种单元采用的是常接触力二节点简单单元，因而很难适应接触面复杂的问题。同时，由于该单元选取节点接触力作为基本未知量，接触应力是由接触力平均得到的，这不仅降低了计算接触应力的精度，而且使得这种简单求应力的方法很难推广至三维问题。雷晓燕于1994年在此基础上提出了一种新的接触摩擦单元。在这种单元中，直接选取节点接触应力作为基本未知量，同时采用六节点的等参单元模拟接触面复杂的几何形状。接触单元的几何和静力约束是作为补充方程包含在刚度方程中的。运用虚位移原理导出接触摩擦问题的等效单元刚度——约束矩阵和等效载荷向量。

1、非线性问题的分类

引起结构产生非线性变形的原因很多，大致可以分为三类：材料非线性问题、几何非线性问题、状态非线性问题。

1）材料的非线性问题：指材料的物理定律即其应力——应变关系是非线性的，但只考虑小位移和小应变的情况，也就是指结构整个几何形状的变化及结构材料内部的应变与结构尺寸相比，是无限小的，这样可以忽略微元体的局部应变。例如，计算应力时可以采用原来的、未变形的微元体面积。至于应变——位移关系则采用线性的应变——位移关系式。各种小应变、小位移问题的结构弹塑性分析及岩土工程中的问题一般都属于这一类。

2）几何非线性问题：虽然假设线性的应力——应变关系，但非线性的应变——位移关系和几何形状的有限变化却引起几何非线性问题。大多数的几何非线性问题是小应变和大位移，当然，也包括大位移和大应变的情况。

3）状态非线性问题：例如，一根只能拉伸的电缆可能是松弛的，也可能是张紧的；轴承套可能是接触的，也可能是非接触的；冻土可能是冻结的，也可能是非冻结的；岩土体中的断层结构面间相对运动。由于系统状态的改变，从而引起系统的刚度在不同的值之间突然变化。状态改变也许和载荷直接有关（如电缆情况），也可能由某种外部原因引起（如冻土情况）。

岩土体材料的一个重要特征是其应力——应变关系具有明显的非线性性质。因此，对于岩土体通常都按材料非线性问题来考虑。材料非线性问题又可以分为两种情况：一种是非线性弹性问题；另一种就是非线性弹塑性问题，系材料超过屈服极限以后就呈现出非线性性质，各种结构和岩土介质若采用弹塑性的本构模型进行分析就是这类问题。在加载过程中，这两种非线性问题在本质上相同，但是卸载过程就会出现不同的现象，非线性弹性问题是可逆过程，卸载后结构或介质会恢复到加载前的位置，非线性弹塑性问题是不可逆的，它将会出现残余变形。

对于材料非线性问题，由于应力——应变关系是非线性的，对于这类问题的微分方程的求解在数学上有一定的困难。然而有限元法处理非线性问题却是十分有效的。用有限元法处理非线性问题的基本思想是用一系列线性问题的解来逐步逼近非线性问题的解。因而，非线性问题可以理解为一系列线性解进行迭代过程的结果。

2材料非线性问题的解法

用有限元解材料非线性问题，通常采用三种基本的方法来解析，即增量法、迭代法以及增量迭代法（或混合法）。在说明问题时，可以只考虑单个单元的非线性平衡方程：

[k （δ）]e{δ}={R} （1）

这是是非线性出现在单元刚度矩阵[k]e 中，k [δ]e是非线性材料性质[D （ε）]的函数。可以把[k]e 写成[k]e =[k{δ}，{R}]的形式，因而[k]e 中的材料参数不再是常数。图1中的载荷{R}与节点位移{δ}相对的非线性应力——应变关系。根据这个应力——应变关系或本构关系，进而确定非线性分析时变化的矩阵[D （ε）] 。本文主要介绍迭代法。

图1 非线性的载荷——位移曲线图

2.1常刚度迭代法

对于方程（1）：[k （δ）]e{δ}={R}，假设R= R0已知，则有限元解法的实质就是在R-δ曲线与R= R0直线的交点的横坐标。为此可用试探法，如果已经找到一个试探值δn，但不满足式（1），如图1所示，它在曲线上的纵坐标值Rn低于R0值，用Pn（R n，δn）表示。现在可以通过点Pn（R n，δn）绘一条直线p，其斜率为kn ，则

kn （δn+1-δn ） = R0 - Rn （2）

这样，直线p与R0 - Rn的交点的横坐标δn+1就是一个改进了的试探解。如果按δn+1 由式（1）计算得的Rn+1仍小于R0，就可以从 Pn+1（R n+1，δn+1）出发，重复前面过程，双可以得到进一步的改进试探解R n+2。这样，循环迭代，直到前后两次的位移值接近为止（图2）。

所谓常刚度，是指取kn=k0为某个常数，即在迭代时直线p的斜率不变，故称为常刚度法。

图2 常刚度迭代法图

2.2割线变刚度迭代法

用一系列的割线，即直线p与R= R0直线的交点P1S，P2S，P3S，……取逐次逼近R-δ曲线与R= R0直线的交点PS时，这些直线p的斜率并不一定相同。如对第n条直线，其斜率为kn，则迭代的方程为

R0= Rn+Kn（δn+1-δn） （3）

因而，可得第n+1次近似解为

δn+1=δn+ kn （δn+1-δn） （4）

这里的割线刚度kn=kns= Rn /δn。其中，kns为从原点出发到所迭代时与曲线交点割母的斜率。因而这个过程称为割线迭代法（图3）。

图3 割线变刚度迭代法图

2.3 切线变刚度迭代法

变切线刚度的表达式可表示为

kn=kn’= （5）

这里kn在物理上代表切线刚度，因而这个过程称为切线变刚度迭代法（图4）。

图4 切线变刚度迭代法图

2.4 New-Raphson 法

对应于非线性方程式（1），可以考虑一个非线性的初始试探解δ0，它与真实解δ*之差为一个小量，任何一个具有一阶导数的连续函数Ψ（δ），在δ0附近光滑，则可按一阶泰勒级数展开为

Ψ（δ）= Ψ（δ0）+ （δ-δ0）+…= 0 （6）

如仅取其线性项，则

Ψ（δ）= Ψ（δ0）+ （δ-δ0）= 0 （7）

所以，δ=δ0- （8）

这个迭代过程，直至在某个允许精度内 <ε。这种迭代方法就称为New-Raphson 法，简称N-R法。

3结语

节理、裂隙、断层及软弱夹层是地质中常见的现象，并且它们都对岩体的强度和稳定性有极大的影响。通常，对于工程岩体而言，有以下两种处理方法：①把节理岩体视作不连续介质，由岩石结构体和结构面组成，因此分别研究岩石和结构面的力学性质及岩石和结构面共同作用时的耦合原理，这种方法的代表有离散无法、块体理论和刚弹性法；②把节理岩体视作宏观上的连续体，建立岩体的等效本构关系，这种方法的代表有当量体法、断裂力学法和损伤力学法。本文根据岩体中含有节理、层理、断层及软弱夹层等不连续结构面，研究了塑性有限元分析中常采用的迭代法的详细过程。