高等岩石力学岩石力学有限元法

高等岩石力学岩石力学有限元法第一页，共三十一页，2022年，8月28日岩石力学有限元法第一节概论第二节施工建造过程的模拟第三节节理及不连续面的模拟第四节多节理岩体的模拟第五节无限域单元及其应用第六节岩体工程中的弹塑性问题第七节无拉力分析及节理非线性分析第八节岩土工程三维非线性有限元程序系列第二页，共三十一页，2022年，8月28日

第一节概论一、有限元法在工程中的应用

有限元自50年代发展至今，以成为求解复杂的岩石力学及岩土工程问题的有力工具，并已为愈来愈多的工程科技人员所熟悉。在求解像弹塑性及流变,动力，非稳态渗流等时间相关性问题，以及温度场在求解像弹塑性及流变、动力、非稳定渗流等时间相关性问题，以及温度场、渗流、应力场的耦合问题等复杂的非线性问题中的效能已使它成为在岩石力学领域中应用最广泛的数值分析手段。特别是近十余年，在工程应用方面已有了较大的进展，并引起广大工程科技人员的兴趣。在岩土工程有关专业的大学生和研究生中，有限元已被列为独立的课程。在本章中将着重讨论有限元法在岩石力学中应用的有关问题。涉及有限元法基本原理，方法及求解技术方面的基础知识。第三页，共三十一页，2022年，8月28日有限元法基本方程就数学概念来说，有限元法是通过变分原理(或加权余量法)和分区插值的离散化处理把基本支配方程转化为线性代数方程，把求解待解域内的连续场函数转化为求解有限个离散点（节点）处的场函数值。显然，这种离散化的处理是一种近似，因面中有当单元划分得充分小时，才能保证满意的求解精度。由于单元充分小，在一个微小的单元内，未知场函数可以采用十分简单的代数多项式近似地表述。通常可取为如下的插值形式：

(4-1)第四页，共三十一页，2022年，8月28日下标m表示单元的节点数目，下标e表示单元的序号。有限元法即是以所有节点处的1值作为基本未知量。对于二维的工程物理问题，如热传导、渗流等，其基本支配方程为如下形式的准谐方程：

(4-2)边界条件为：①，边界在上②边界在上第五页，共三十一页，2022年，8月28日当=常量时，上式即退化为波松方程。由变分原理可知，满足方程式（4-2）及相应边界条件的场函数，应使如下的泛函数驻值，即：

(4-3)在对被考察域A进行离散化的情况下，此泛函驻值转化为多元函数的驻值：（4-4）

第六页，共三十一页，2022年，8月28日对于一个单元体的子域内，把方程式（4-1）带入方程式（4-3)并进行变分后即可得到单元的基本方程：

(4-5)

式中即单元特性矩阵：即为节点未知量组成的矢量：是由在上的边界条件给定的a（x,y）有关的矩阵，故仅对a≠0的二类边界上的单元才予考虑。对整个求解域A由变分驻值条件式（4-4）可得或写为

(4-6)第七页，共三十一页，2022年，8月28日

式中，[k]为由各单元特性矩阵及按节点号组集得到的总体特性矩阵：{u}为所有节点的待求值组成的矢量：{p}为与Q（x,y）有关的矢量及q(x,y)有关的矢量。解线性方程式（4-6）即可求得场函数在各单元节点处的值。对于弹性力学问题，可通过最小势能原理或虚功原理导出有限元法的基本方程。有限元法求解弹性力学问题通常以位移作为基本未知量。位移在直角坐标系中沿坐标轴X,Y,Z的分量分别表示为，故单元位移模式这时可写为

(4-7)式中，l为3×3阶单位矩阵；Nt称为单元的位移插值函数或形函数（对二维问题l为2×2阶，位移分量）；{}为由单元各节点位移分量形成的矢量称为移矢量，利用弹性力学的几何方程及物理方程可导出单元的应变及应力表达式：第八页，共三十一页，2022年，8月28日

（4-8）（4-9）应用最小势能原理或虚功原理可以推导出单元刚度矩阵的表达式

（4-10，a）

（4-10，b）

（4-10，c）第九页，共三十一页，2022年，8月28日

式中，{q}为分布体力；{p}为分布面力，通常为坐标的已知函数，最简单也最常见的情况是{q}及{p}均为常量，上式的积分则可简化。对于一般准谐方程描叙的工程物理问题及弹性力学问题，有限元法的公式表述及分析方法具有完全相同的格式，这是有限元法的主要特点。这种统一的格式有利于实现具有广泛通用性的计算机程序。二、岩体力学问题的特点

大多数岩土工程问题如结构-岩体相互作用，岩土边坡、地下工程等，都涉及无限域或半无限域。处理这类问题通常是在有限的区域进行离散化。为了使这种离散化不会产生大的误差，必须取足够大的计算范围，并应使假定的外边界条件尽可能接近真实状态。第十页，共三十一页，2022年，8月28日理论分析和计算实践表明，当由于结构或工程开挖释放荷载作用于岩体某一部位时，对周围岩体的应力及位移有明显影响的范围大约是开挖或结构物与岩体作用面的轮廓尺寸的2.5～3倍，在此范围之外，影响甚小，可忽略其影响。考虑到有限元离散误差和计算误差，为了保证必要的计算精度，计算范围应取不小于3～4倍，如图4-1所示。在这种情况下，外边界可采取两种方式处理：其一为在距离荷载作用足够远的外边界处位移为零；其二为假定外边界为受力边界（包括p=0）。这两种方式都同实际的无限域不完全一致，因而都有一定误差。这种误差随着计算区域的减小而增大，并且在靠近外边界处都远离外边界处的误差大，我们称此现象为“边界效应”。第十一页，共三十一页，2022年，8月28日三、初始地应力场与释放荷载由于自重及地质作用，自然状态的岩体处于一定的初始地应力状态，在结构荷载作用图4-1有限元离散化的计算范围第十二页，共三十一页，2022年，8月28日

下，岩体内的应力为荷载产生的应力与初始地应力之和。因此，初始地应力的大小直接影响到计算结果。正确决定岩体的初始应力状态是有限元分析中的一个重要问题，至今未能得到妥善的解决。由于初始地应力的存在，岩石的开挖将导致部分岩体卸载，成而使一定范围内的应力场发生变化。正确模拟这种开挖效果也是岩体力学问题的一个重要的特点，把由此所得的位移即作为由于工程开挖长生的岩体位移，由此所得的应力场同初始地应力场迭加即为开挖后的应力场，这种模拟开挖效果的方法可如图4-3所示。“开挖释放荷载”是按照已知的初始地应力由等效原则来确定的。

第十三页，共三十一页，2022年，8月28日第二节施工建造过程的模拟本节主要讨论如何正确模拟开挖及工程建造的过程一、开挖释放荷载对于具有已知初始应力场的岩体开挖工程，沿预计开挖线上各点的初始应力为已知，在进行离散化的情况下，可假沿开挖边界上两相邻节点之间的初始应力呈线性变化，则对于任一开挖边界点，开挖引起等效释放荷载为

(1)第十四页，共三十一页，2022年，8月28日

式中，上标i，i-1及i+1为沿开挖边界上的有限元网格的节点号如果坐标X，Y轴正沿主应力轴，则上式可简化为

(2)若原始应力场均为应力场，寄节点i、i-1、i+1等各点应力场相等，则上式可简化为

(3)第十五页，共三十一页，2022年，8月28日若X，Y轴同应力主轴重合，可进一步简化为

(4)考虑到存在一个初始应力场的情况，开挖后的实际应力场为初始应力场与开挖释放应力场之和。即以上是在已知开挖边界上各点初始应力的情况下计算各边界节点上的等效释放荷载，边界节点上的初始应力可按下述方法确定：第十六页，共三十一页，2022年，8月28日1、对于已知均匀初始应力场，有等效释放荷载可由(3)式求得

2、对于非均匀初始应力场，一般先用有限元法计算初始应力场，然后把要挖去的区域内各单元改变为“空单元”计算开挖结果当有限元计算直接给出各单元节点应力时，则开挖边界上各节点的应力也为已知。利用（2）式可以求出。第十七页，共三十一页，2022年，8月28日二、施工过程模拟为了模拟施工建造过程，在确定离散化网格时，必须要考虑各步施工的情况及结构特征。如下图为一地下洞室，开挖及支护分上、下两部分进行。施工步骤如下：

图4-7洞室断面及施工步骤第十八页，共三十一页，2022年，8月28日第三节

节理及不连续面的模拟一、平面问题节理单元——Goodman单元采用特殊的“节理单元”模拟地质不连续面是由R.EGoodman(1976)最先提出来的。一种无厚度4节点节理单元也常被称作Goodman单元，如下图所示，它是一段直接接触的平面，承受沿此界面的切向力及法向力。节理的应力-应变关系对局部坐标（s,n)可写为也可简化为第十九页，共三十一页，2022年，8月28日式中，分别为切向应力及法向应力；，为节理的切向及法向刚度；，为单元两侧对应点的相对切向及法向位移。

无厚度节理单元假定位移沿节理单元长度呈线性变化，则可导出节理单元对于局部坐标s-t的单元刚度矩阵为第二十页，共三十一页，2022年，8月28日对整体坐标的转换为式中，[T]为坐标变换矩阵，它是具有如下的分块形式的对角阵：

式中

上述节理单元的缺点是由于无厚度，计算中可能发生节理上下面相互“嵌入”的现象，必须对这种嵌入量做人为的限制，以免导致较大误差。此外，当节理的上下面发生相对转角位移时，也讲产生误差。Goodman曾于1976年对节理单元做了修正，考虑到相对转角的影响，假定以节理单元的中点作为计算点，则其应力-应变关系表示为

第二十一页，共三十一页，2022年，8月28日式中，为节理中点力矩；为相对转角，以节理的3—4节点的平面逆时针为正。仍假定位移沿长度线性变化，可导出节理单元刚度矩阵为二、变厚度节理单元下图所示的一种六节点变厚度的节理单元，具有更广泛的适用性，可用于等厚度、变厚度、平面及曲面的节理，这种单元实际上是细长四边形的等参单元，单元的位移函数及坐标变换关系可写为第二十二页，共三十一页，2022年，8月28日

这种较窄长的单元，在和方向采取不同阶次的插值函数(

方向为二次，方向为一次），即变厚度节理单元第二十三页，共三十一页，2022年，8月28日对于有一定厚度的夹层，可采用一般平面应变问题的本构关系，按一般平面等参单元处理，其单元刚度矩阵为式中det即为雅可比转换。应用高斯积分完成单元刚度计算，在方向可取三点积分，在方向可取两点积分。式中由几何方程导出，弹性矩阵可取一般各向同性或正交异性弹性矩阵，当采用正交异性弹性矩阵时要注意坐标至整体坐标之间的变换。当厚度很小时，可以作为等厚或无厚度的节理单元。此时各应力分量沿厚度方向为常量，则单元退化为上图的形式。插值函数可简化为第二十四页，共三十一页，2022年，8月28日三、三维问题节理单元

由二维节理单元可方便的推广至三维问题，常用三维节理单元如图所示单元位移插值函数可表示为

等厚度三维节理由单元第二十五页，共三十一页，2022年，8月28日对图a插值函数为图b插值函数为座标变换关系为第二十六页，共三十一页，2022年，8月28日考虑节理上、下面一对点的相对位移，对整体座标有式中为3×3阶单矩阵，为对整体座标的节点位移，通常需要考虑节理沿其法向及切向的相对位移及应力。为此，取节理中面任一点的法线方向作为轴建立局部动座标系则平面即为过单元该点的切平面，这时有