《岩石力学的数值模拟(讲义)(共26页).doc》由会员分享，可在线阅读，更多相关《岩石力学的数值模拟(讲义)(共26页).doc（26页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上第10章 岩石力学的数值模拟随着计算机软硬件技术的迅速发展，使岩石力学有了长足的进步，特别在岩石力学的数值计算和模拟方面发展尤为迅速，使得许多岩石力学解析方法难于解决的问题得以重新认识。正如钱学森在给中国力学学会“力学迎接21世纪新的挑战”的一封信中对力学发展趋势总结的那样“今日力学是一门用计算机计算去回答一切宏观的实际科学技术问题，计算方法非常重要”。岩石力学和其他力学学科一样，需要数值计算方法并推动岩石力学的发展。岩石介质不同于金属材料，在数值计算方面具有其独特的特点205：（1）岩石介质是赋存于地壳中的各向异性天然介质。（2）岩石介质被众多的节理、裂缝等弱面所切
2、割而呈现高度的非均质性，而其物理、化学及力学性质具有随机性特点。（3）岩石介质赋存时以受压为主，而且抗压强度远大于抗拉强度。（4）岩石力学与工程问题在时空分布上较广，从本质上讲都是三维问题。（5）岩石工程一般无法进行原型试验，而实验室测得的数据不能直接应用于工程设计和计算。（6）岩石力学与工程具有数据有限问题。数值计算方法经过几十年的发展，目前已形成许多种岩石力学计算方法，主要有有限元法、边界元法、有限差分法、离散元法、流形元法、拉格朗日元法、不连续变形法及无单元法等。它们各有优缺点，有限元的理论基础和应用比较成熟，在金属材料和构件的计算中应用十分成功，但它是以连续介质为基础，似乎与岩体的非连
3、续性有一定差距，流形元等数值方法虽然考虑了岩体中节理效应，但其理论基础还不完全成熟。相信在不久的将来，肯定会出现完全适合于岩体材料和工程的数值计算方法206208。10.1 岩石力学的有限元分析209213有限元法（finite element method，FEM）是岩石力学数值计算方法中最为广泛应用的一种。自20世纪50年代发展至今，有限元已成功地求解了许多复杂的岩石力学与工程问题。被广大岩石力学研究与工程技术人员喻为解决岩石工程问题的有效工具。有限元法是根据变分原理求解数学物理方程的一种数值方法。有限元法把连续体离散成有限个单元，每个单元的场函数只包含有限个节点参量的简单场函数，这些有限
4、个单元的场函数集合构成整个结构连续体场函数。根据能量方程和加权函数方程可建立有限个求解参数的方程组，求解这些离散方程组，就是有限元法的精髓所在。虽然求解时把连续函数转化为求解有限个离散点处的函数值，但只要单元划分得充分小时，足可以满足计算要求。有限元法求解问题时一般遵循以下步骤：（1）有限元计算模型的建立，包括模型单元的划分、确定边界条件。（2）对单元体进行力学分析，包括求解节点位移、单元应变和单元应力。（3）对计算模型进行分析。（4）进行计算分析。10.1.1 线弹性有限元法的基本方程线弹性有限元是弹塑性有限元、损伤有限元、流变有限元等非线性有限元的基础。线弹性有限元假定岩石介质连续、均质、
5、小变形和完全弹性。有限元法求解弹性力学问题时通常以位移作为基本未知量，单元位移是以单元节点位移为基本未知量，选择合理的位移插值函数，将单元位移表达为节点坐标的连续函数，插值函数也可称为形函数。不同形状的单元具有不同的形函数。图10-1为三种最常见单元形式，即三角形、四边形及四面体单元。它们的形函数分别为：图10-1 有限元的三种基本单元形式（a）三角形单元（b）四边形单元（c）四面体单元三角形的形函数式中，S为三角形面积；。四边形的形函数式中，位移量为；。四面体的形函数 式中，V为四面体的体积。单元在直角坐标轴中位移分量分别为u，v，因此单元的位移矩阵为 （10.1）式中，为单元位移矩阵；N为
6、形函数矩阵；为单元节点位移列阵。根据几何方程，对位移矩阵求偏导数，可以得到应变矩阵 （10.2）式中，B为连续单元节点位移和单元应变的矩阵，也称为应变矩阵。对于三角形单元，B为常数矩阵，元素值取决于单元节点坐标差。根据本构方程，可以得到单元节点位移与单元应力矩阵之间的关系 （10.3）式中，D为弹性矩阵。应用虚功原理和最小势能原理可以推导出单元刚度矩阵的表达式 （10.4）各单元的体积力和面力按照静力的等效原则移置到各单元的节点上，其等效节点力为 （10.5）式中，Pe为作用于单元体积力P的等效节点荷载；Qe为作用于单元面积力Q的等效节点荷载。设环绕某节点i共有k个单元，则i节点上的外加荷载R
7、i为 （10.6）式中，Pi为作用于i节点上集中力。将各单元节点力与节点位移之间的关系叠加，形成以节点位移列阵为基本未知量的线性代数方程组 （10.7） （10.8）求解上式有限个线性代数方程组，得到节点位移矩阵，根据相应的节点位移利用式（10.2）和（10.3）计算单元的应力应变值。10.1.2 非线性问题的处理方法从本质上讲，岩石力学问题都是非线性问题。这是因为一方面岩石材料的应力应变本构关系绝大多数呈现非线性特征，另一方面岩体的变形又大多是大变形。对于求解岩石力学的非线性问题，解析方法显得无能为力，而有限元可以求解岩石力学绝大部分的非线性问题。岩石力学的非线性问题可以分为三大类：材料非线
8、性，即岩石材料的本构关系为非线性，而变形的几何关系仍是线性的；几何非线性，即岩石几何变形为非线性，而本构关系仍为线性；两类非线性问题的组合，即岩石材料既是材料本构非线性，又是几何非线性。这三类非线性问题总体平衡方程组的共同特征都是非线性方程组，可表示为 （10.9）式中，为刚度矩阵，它是位移的函数。求解这类非线性方程组一般有三种方法：迭代法、增量法和混合法。迭代法又称为牛顿-拉斐逊法，对于一个变量的函数，如图10-2，它的迭代过程如下：设函数值F由F0增加至Fa通过切线法做第i次近似值可由下式确定： （10.10）式中，Ki-1为第i-1次迭代时的曲线切线斜率，那么最终的解为 （10.11）牛
9、顿-拉斐逊法的主要缺点是每次迭代过程中都要重新计算一下切线值，也就是刚度矩阵及其逆矩阵，因此花费机时较长，为了避免每次计算Ki值，每次迭代时都采用同一个初始的K0，如图10-3，此方法称为修正的牛顿-拉斐逊法，具体的计算公式如下 （10.12） 图10-2 牛顿-拉斐逊迭代法 图10-3 修正的牛顿-拉斐逊法两种迭代法比较而言，修正的牛顿-拉斐逊法的迭代次数多于牛顿-拉斐逊法，但它省去了每次迭代时重新计算赐度矩阵Ki。计算时间的多少，不仅取决于迭代次数或收敛速度，还取决于每次迭代所花费的时间，在一般情况下常刚度法在总体计算时间上比较合理，对于非线性很强的材料，常刚度法并不适用。增量法也是把非线
10、性问题线性化的一种处理方法。如图10-4把总荷载F划分成若干个增量，逐级施加荷载增量进行求解，有限元计算公式为 （10.13）那么总位移和总荷载为 （10.14）增量法的误差较大，最终误差为各级增量时误差之和，为了减小误差，有时对一级增量后，加上一个修正正项加以误差矫正，计算公式为 （10.15）图10-4 增量法由于增量法和迭代法各有其优缺点， 目前有限元常常采用增量法和迭代法的结合，即混合法，一般将非线性关系分成若干个荷载增量段，在每一个增量段内用迭代法求解逼近非线性解，最终累加求得各级荷载作用下的最终应力与应变。10.1.3 岩石弹塑性有限元分析12岩石弹塑性本构关系是岩石主要非线性问题
11、之一。岩石的弹塑性是指岩石材料的应力应变关系在屈服之前呈线性关系，当应力达到屈服应力时，应力应变关系就变为非线性。由于弹塑性模型中应变不仅依赖于受载的应力状态，而且与加载路径有关，因此一般弹塑性本构关系不能用应力应变全量关系准确描述，只能用能反映加载路径有关的应力应变增量关系描述。在岩石非线性弹塑性本构关系有限元分析中一般采用初应力法和初应变法求解非线性平衡方程组。初应力法是将荷载以微小增量形式逐级加在模型上，每加一级荷载增量，就会产生相应的位移增量、变形增量和应力增量。对于具有初应力的弹塑性应力应变本构关系可以写成 （10.16）式中，为初应力；为塑性矩阵，它与加载前的应力水平有关，而与应力
12、增量无关。初应力法是通过对的处理将应力修正到正确的水平上，初应力不仅与加载增量前应力水平有关，还与本级所加荷载增量引起的变形增量有关，如图10-5。增量形式的平衡方程为 （10.17）图10-5 初应力法式中，K0为线弹性计算中的总刚度矩阵；为矫正荷载项，它由下式决定： （10.18）由于随位移变化而变化，因此计算时必须进行迭代求解。初应力法求解按照下述计算步骤实现：（1）把全部荷载划分成若干个增量，在每一级增量段内，按照增量弹塑性平衡方程进行求解。（2）计算各单元的应力增量及当前的应用 （10.19）式中，下标i表示第i级荷载增量；j表示第j次迭代。（3）根据岩石的屈服准则，由各单元应力判断
13、单元是否屈服，对于塑性单元，计算应力修正项并修正应力 （10.20）（4）塑性单元通过修正项计算等效节点力，所有塑性单元的等效节点力叠加构成总的修正荷载矢量 （10.21）（5）在修正荷载作用下进行下次迭代运算，此时基本方程为 （10.22）重复进行（2）（5）步计算直至所有的塑性单元达到收敛精度要求。再进行下一步的荷载增量计算。（6）重新施加下一级荷载增量，重复计算（1）（5）步，直至计算完毕。通过累加各级荷载作用的计算结果，求得总位移和总应力。一般初应力法的收敛速度比较缓慢，因此通常采用常刚度和变刚度法相结合的方法加速收敛。初应变法按照弹塑性增量理论，总的应变量可表示为 （10.23）式中
14、，为弹性应变增量；为塑性应变增量，类似初应力法可以把塑性应变增量看做初应变。因为在计算过程中的计算格式和弹性系统中的初应变一致。图10-6 初应变法从图10-6中可以看出，荷载增量dF所对应的位移为，按线弹性计算为，两者的位移增量之差称为初始位移，它是由材料塑性引起的附加位移。与模型系统初始位移对应的单元位移为，那么单元的初应变为 （10.24）10.1.4 岩石节理单元的有限元方法岩体中常存在大量节理，而节理的变形性质和岩体力学性质有十分明显差别。从某种程度上讲，节理的存在决定着岩体的力学性质。因此分析节理的变形性质，对于研究岩体工程的变形情况至关重要。在进行有限元分析时，这种节理一般采用专
15、门的节理单元来处理。节理单元首先由Goodman提出19，214，Goodman认为节理不是一个实体，它只在若干点处接触，因此节理单元也应满足这一特点，图10-7为无厚度节理单元，节点1与4，2与3具有相同的坐标，沿节理面的应力分别为法向应力和剪切应力。把节理单元两侧对应的位移定义为应变，相对的法向位移称为法向应变，相对的切向位移称为切向应变，那么节理单元的本构关系为 （10.25）图10-7 无厚度节理单元式中，为单元的刚度矩阵，对于节理单元长度上任何一点s处的应变可定义为界面两侧相应位移之差 （10.26） （10.27）上式可用矩阵形式表示为根据一般化公式导出节理单元对应于局部坐标的刚度
16、矩阵 （10.28）式中，分别为法向和切向刚度。对于整体坐标的刚度矩阵通过坐标变换矩阵得到，具体如下 （10.29）式中，T为变换矩阵，它具人分块形式的对角阵 （10.30）式中，。这种无厚度节理单元适用于模拟直接接触的界面，而对于有一定厚度的弱面或夹层不适宜。10.2 岩石力学的边界元分析边界元法（boundary element method, BEM）是和有限元法同步发展的一种数值计算方法。与有限元相比有以下特点：边界元法把一个均质区域看做一个大单元，只把它的边界离散化，区域内不划分单元，场变量处处连续；对于无限区域，场变量自动满足无穷边界条件及自然表面状态。对于半无限区域，场变量也自动
17、满足无穷远边界条件及自然表面状态。有限元法是全区域离散化，而边界元是把基本方程转化为边界积分方程，只对边界离散化建立相应的方程组进行求解。这样边界元使三维问题降为二维问题求解，使二维问题转化为一维问题求解，当物体的表面积和体积之比较小时，边界元的划分单元数要比有限元少数倍和十几倍，这样也使待解的方程数目、处理和存储的数据量降低同样的倍数，大大节省了机时和算题费用。当仅需求解物体内部几个点的应力时，有限元仍不得不划分整个区域，才能确定这几个点的应力值。而边界元当知道边界的应力解时，就可以根据需要去求物体内部个别点的解。在应力梯度较高处，有限元法的剖分密度常常受到限制，而边界元可以方便地确定应力梯
18、度的分布。但随着计算机硬件的飞速发展，边界元的优势逐渐显得不很明显。边界元法矩阵方程中系数阵的元素结构比有限元法刚度阵中的元素复杂。有限元刚度阵属带状稀疏阵，而边界元法的系数阵为满阵，因此对于面积和体积之比较大的薄壁结构而言，边界元的优越性就不明显。边界元比较适合求解无限区域和半无限区域问题，如深埋巷道是一个典型的例子。但边界元在计算非均质介质问题、非线性问题以及模拟工程开挖过程等方面不如有限元方便有效。有限元与边界元划分单元的比较如图10-8。图10-8 有限元法与边界元法比较（a）力学模型和边界条件（b）有限元单元划分（c）边界元单元划分求解边界方程组主要有Massonet和Rizze分别
19、提出的直接和间接两种数值解法。直接法是直接建立关于边界的积分方程，通过离散化求解边界未知参数，并进而求解计算区域内场函数值。间接法是设立一个在求解区域内含有若干系数的满足基本方程组的解，代入边界条件求出系数，进而求得边界及区域内各点的场函数值。10.2.1 边界元分析原理215217在无限的弹性体内作用一单位力引起周围的应力和位移的解析解是边界元求解弹性体问题的基础，如图10-9在平面无限体内i点作用一x方向的单位力，其基本的弹性解在1848年由Kelvin解出 （10.31） （10.32）式中，为i点l方向单位力引起j点k方向的应力和位移；为Kronecker符号，当l=k时，当时，；r为
20、i，j两点之间的距离（矢径）；，为j点作用面外法向对于k和l的方向余弦；为矢径r对外法向n的方向倒数；。当不考虑体力时，根据功的互等定理和以上的Kelvin基本解，可以建立弹性体边界上i，j两点之间的应力和位移关系（如图10-10） （10.33） 图10-9 开尔文基本解条件 图10-10 边界积分方程的建立式中，的大小取决于边界情况，当边界为光滑曲线时，等于0.5；当边界曲线不光滑时，的值根据刚体位移原则求解。当边界作用分布力时，j点绕行一周，按叠加原理积分得 （10.34）式（10.34）就是边界元中的边界积分方程。10.2.2 边界单元与边界的离散边界元是通过离散边界来求解边界积分方程
21、。对于一个确定区域的边界进行离散，划分成有限个线段，称为边界单元。根据单元的节点数目和节点模式，边界单元可分为常量单元、线性单元、二次单元等。常量单元及线性单元均为直线段，常量单元以单元的中点为节点，每个单元有一个节点，场变量在单元内是一个常数。而线性单元的场变量在单元内呈线性变化。单元的模式与有限元相似，可以表示成插值函数形式，即， （10.35）式中，m为单元的节点数目；插值函数由拉格朗日插值公式给出；为节点i处的值。如图10-11几个常见的插值函数如下：常量单元， 线性单元， ， 二次单元， ， ， 图10-11 常见的几种边界单元（a）常量单元（b）线性单元（c）二次单元对于常量的边界
22、元，边界积分方程可简化为 （10.36）对于平面问题，上式有2n个方程，可写成矩阵形式 （10.37）式中，分别为常量边界单元中点的位移列阵和应力列阵；G，H为系数矩阵。边界是元法的边界条件一般都为混合边界条件，即在一部分边界上位移和应力已知，而另一部分未知。为了方程求解，把所有的已知量移到等式的右边，未知量移到等式右边，经过变换后式（10.37）可简化为 （10.38）式中，X为包含所有位移和应力未知量列阵；F为已知量列阵；A为系数矩阵。求解此方程可求得全部边界节点上未知量，由此进而求得计算域内任意一点的位移和应力值。如图10-12，根据开尔文基本解和Betti的互等定理，可以得到计算域内任
23、意一点的位移和边界点外力功的关系式 （10.39）图10-12 计算域内i点与边界j点作用力与位移的关系图边界离散后，对于常量单元，上式可改写为 （10.40）欲求在i点l方向的位移分量，可建立边界积分方程 （10.41）式中，分别为i点l方向作用单位力于j点k方向引起的应力位移开尔文基本解。欲求计算域内i点的应力可由几何方程和Lame方程求得 （10.42）将边界积分方程代入上式可得 （10.43）对于二维问题，上式中系数D，S分别为 （10.44） （10.45）式中，。10.2.3 边界元与有限元耦合如前所述，边界元和有限元在计算时各有优缺点，为了发挥各自的优点，提高求解精度和解题效率，
24、A.Wexler于1972年提出了边界元和有限元耦合求解的思路和方法。以后，J.R.Osias和M.Margulies等对边界元和有限元的耦合求解进行了深入的研究。为了讨论方便，把有限元的基本方程改写为如下形式 （10.46）式中，F为边界力的等效节点力；D为体力的等效节点力。由于边界力可以由面力P与面积之积获得，那么有限元方程可改写为 （10.47）一般边界元方程在考虑体力的情况下也可写为 （10.48）从式（10.47）和（10.48）可以看出，有限元方程和边界元方程具有类似的形式，因此从解题方式上提供了建立两者耦合的基本条件。把一个计算域人为地分成两个子域，对两个不同的子域分别采用边界元
25、和有限元进行求解。如图10-13有限元的边界为S2，边界元的边界为S1，两者的共同边界为S3。对两个区域分别建立有限元方程和边界元方程，再利用两者边界S3上的平衡和协调条件建立耦合方程组。耦合方程组有两种方式建立：把边界元方程转化为等价的有限元方程；把有限元方程转化为等价的边界元方程。一般在耦合求解时前者采用得较多。图10-13 边界元和有限元耦合的区域划分边界元方程转化为有限元方程时，首先对边界元方程进行改造，对式（10.48）中G进行求逆可得 （10.49）因为边界元法中P为分布力，要转化为有限元中的等效节点力，必须在上式两侧左乘一个M矩阵，那么可得到与有限元形式一致的方程 （10.50）
26、式中； （10.51）需要注意的是有限元法的刚度矩阵为对称矩阵，而由边界元方程转化来的为非对称，为了完全耦合，必须通过最小二乘法使其对称化。对称化后的刚度系数为 （10.52）一个建筑结构在地基上沉降问题比较适合耦合求解，如图10-14，把建筑结构部分划分为有限元，地基部分为边界元。边界元的离散化有以下两种形式，由地面延伸到到无限远处，边界元的划分须近似截取一个有限区，或利用半无限平面问题的基本解引入到无限边界元。本算例采用耦合成有限元方程进行求解。计算结果和有限元结果比较如表10-1。可以看出耦合法和有限元计算结果十分一致，但计算时间大为减少。表10-1 结构顶部垂直位移计算结果比较（mm）
27、有限元解-339-97135361600耦合法解-350-105135370617 图10-14 有限元和边界元耦合求解的算例1210.3 岩石力学有限差分方法（FLAC）有限差分方法（finite difference method，FDM）是从一般的物理现象出发建立相应的微分方程，经离散后得到差分方程，再进行求解的方法。这种方法可能是最古老的求解方程组的数值方法之一。差分方程在计算机出现以前用一般的手摇计算器也可以求解。随着计算机的不断发展和其他计算方法的兴起，有限差分法曾一度受到冷遇，但20世纪80年代末由美国ITASCA公司开发的FLAC（fast Lagrangian analysi
28、s of continua）程序广泛采用差分方法进行求解，在岩土工程数值计算中得到了广泛的应用。10.3.1 FLAC程序的简介218FLAC是为岩土工程应用而开发的连续介质显式有限差分计算机软件，主要适用于模拟计算岩土类工程地质材料的力学行为，特别是岩土材料达到屈服极限后产生的塑性流动。材料通过单元和区域表示，根据研究对象的形状构成相应的网络结构。每个单元在外载和边界约束条件作用下，按照约定的线性和非线性应力，应变关系产生力学响应。FLAC软件建立在拉格朗日算法基础上，特别适用于模拟材料的大变形和扭曲转动。FLAC程序设有多种本构模型，可模拟地质材料的高度非线性（包括应变软化和硬化）、不可逆
29、剪切破坏和压密、黏弹性（蠕变）、孔隙介质的流固耦合、热力耦合以及动力学行为等。另外，程序还设有边界单元，可以模拟断层、节理和摩擦边界的滑动、张开和闭合行为。支护结构，如衬砌、锚杆、可缩性支架或板壳等与围岩的相互作用也可以在FLAC程序中进行模拟。同时用户可根据需要在FLAC中创建自己的本构模型，进行各种特殊的修正和补充。FLAC程序具有强大的后处理功能，用户可以直接在屏幕上绘制或以文件形式创建和输出多种形式的图形，使用者还可以根据需要，将若干个变量合并在同一幅图形中进行研究分析。由于FLAC程序主要为地质工程应用而开发的岩石力学数值计算程序，包括反映地质材料力学效应的特殊要求。FLAC设计有7
30、种材料本构模型：（1）各向同性弹性材料模型。（2）横观各向同性弹性材料模型。（3）Coulomb-Mohr弹塑性材料模型。（4）应变软化、硬化塑性材料模型。（5）双屈服塑性材料模型。（6）节理材料模型。（7）空单元模型，可用来模拟地下开挖和煤层开采。FLAC程序采用的是快速拉格朗日方法，基于显式差分来获得模型的全部运动方程（包括内变量）的时间步长解。程序将计算模型划分为若干个不同形状的三维单元，单元之间用节点相互连接。对某一个节点施加荷载之后，该节点的运动方程可以写成时间步长的有限差分形式。在某一个微小的时间内，作用于该点的荷载只对周围的若干节点（相邻节点）有影响。根据单元节点的速度变化和时间
31、，程序可求出单元之间的相对位移，进而可以求出单元应变；根据单元材料的本构方程可求出单元应力。随着时间的推移，这一过程将扩展到整个计算范围，直到边界。这样程序可以追踪模型从渐进破坏直至整体破坏的全过程。FLAC程序将计算单元之间的不平衡力，并将此不平衡力重新加到各节点上，再进行下一步的迭代运算，直到不平衡力足够小或者各节点的位移趋于平衡为止。图10-15为FLAC程序的计算循环示意图。FLAC和有限元相比具有以下一些特点：（1）由Maxti和Cundall于1982年提出的混合离散法适用于塑性破坏荷载和塑性流动的模拟，这个方法比有限元中普遍采用的归约积分法更为合理。图10-15 FLAC程序中的
32、计算循环示意图（2）FLAC虽然具有动力方程，但模拟系统实质上是静止的。这使得FLAC不需要数值上卸载就可以遵循物理的失稳过程。（3）FLAC采用显式表达方式，对于任意非线性应力应变关系计算所用的时间和线性关系基本一致，而且它并不需要存储任何矩阵，这就意味着计算机一般的内存就可以计算大量的单元，而且大变形计算所花费的计算和小变形基本一样，因为它没有刚度矩阵。当然FLAC也有不足之处：用FLAC计算线性问题比同样的有限元要慢，FLAC对非线性、大变形问题及岩土物理失稳的计算更有效；FLAC的计算时间由模拟系统最长固有周期和最短固有周期之比确定。对于一些特定的模型，如弹性模量和单元尺寸变异较大的模
33、型，计算效率非常低。10.3.2 有限差分拉格朗日法的基本原理和方程在有限差分方法中，控制方程组的每一个变量都可以直接用离散点的场变量代数形式直接表达，其特点是直接求解基本方程和相应的定解条件近似解。具体步骤是首先将求解域划分成网络，然后在网络节点上用差分方程近似表示微分方程，当采用较多节点时，近似解的精度可以得到很大的提高。而有限元是将连续求解区域离散成一组有限个相互连接的单元组合体，利用每一个单元内假设的近似函数来分片求解区域上待求的未知场函数。单元内的近似函数通常由未知场函数及其导数在单元各节点的数值和插值函数来表达。这样有限元分析中的未知函数及其导数在各个节点上的数值就成为新的未知量，
34、一经求出这些未知量，就可以通过插值函数计算出各单元内场函数的近似值，进而得到整个求解域上的近似解。有限差分和有限元两种方法都有一组代数方程需要求解，虽然那些方程是用不同的方法推导而得，但其最终的形式相同，一般有限元程序把单元矩阵组成一个整体的总刚度矩阵，而有限差分方法并不需要，因为在每一计算步都产生一个有限差分方程。FLAC采用拉格朗日分析方法，由于它不需要形成整体刚度矩阵，大变形计算时在每一个计算步都很容易修正坐标。位移增量施加到坐标上以致网格随着材料的移动和变形，这就是所谓的拉格朗日法，与此对应的欧拉法，其材料移动和变形是相对于固定的网格。10.3.2.1 连续介质的场方程对于一个二维的连
35、续介质而言，其运动方程可表示如下 （10.53）式中，为介质密度；为坐标；为重力加速度。在FLAC程序中，单元应变率是计算各主要参数的纽带，由于非线性应力应变本构关系不具有惟一性，一般用增量形式表示，各向同性最简单的弹性本构关系张量差分形式如下 （10.54）式中，G，K分别为剪切和体积模量；为时间步。10.3.2.2 差分方程三角形单元的差分方程一般由高斯离散定律推导而得 （10.55）式中，为沿封闭表面边界的积分；f为标量、矢量和张量；xi为位置矢量；为面积A上的积分；ni为表面s上的单位法向量。把面积A上f的平均值定义为 （10.56）把上式代入式（10.55），就可以得到 （10.57
36、）对于一个三角形单元，其有限差分形式可变为 （10.58）式中，为三角形的边长；为每条边上的平均值。由引可以把应变率用节点速度来表示 （10.59）式中 （10.60）式中，标记（a）和（b）为三角形边的两个相邻节点，把应变率代入本构关系式（10.54），就可以获得应力张量。一旦应力计算出来，相应施加到各节点上的力也就确定。如图10-16，每个三角形节点由两个相邻边上应力合成得到，因此 （10.61）对于包含两个三角形的四边形单元，由一个三角形形成的节点，其节点力是三角形两边应力合成之和，有两个三角形形成的节点，其节点力是两者的平均值。划分网络中节点力是其周围所有单元对该节点的矢量和，它包含单
37、元体重力和作用在节点上的外加荷载。如果单元体处于平衡或稳定状态流动（塑性流动），那么该节点力将为零，否则该节点有一个加速度 （10.62） 图10-16 三角形单元节点力的计算模型（a）具有速度矢量的三角形单元（b）节点力矢量式中，上标表示相应变量计算的时间。对于大变形问题，上式再积分一次以确定网格节点的新坐标 （10.63）有关FLAC计算实例请见文献219，220。10.4 岩石力学的离散元分析有限元法、边界元法和有限差分法都是基于连续介质力学的数值计算方法，它们都要求计算对象满足变形的连续性条件。但工程岩体往往被节理和结构面所切割形成明显的节理岩体，具有明显的不连续性，用连续介质力学的数
38、值计算方法难于处理。针对不连续岩体的变形和运动的求解，Cundall于1971年提出了一种新的数值计算方法离散单元法（distinct element method，DEM），并用汇编语言编制了计算程序。1978年Cundall又将最初的汇编语言全部翻译成FORTRAN文本，形成离散元的基本程序。到1985年，他们完成了目前广泛采用的离散元数值分析程序UDEC（universal distinct element code）。离散元法由王泳嘉教授等人于20世纪80年代中期介绍到我国，发展十分迅速。目前在矿业、铁道及水利等行业已被广泛应用221-223。10.4.1 离散元法基本原理有限元法的理
39、论基础是基于最小势能原理，边界元法的理论基础是Betti互等定理，而离散元法的理论基础则为最简单、最基本的牛顿第二运动定律。它与其他数值方法不同的是离散元法首先将计算区域划分成有限个独立的多边形块体单元，单元与单元之间具有一定的初始接触状态，单元在外力或自重的作用下转动和移动，最终达到平衡状态或匀速运动。离散元法划分的单元，依据力学性质可以分为刚性体和变形体，依据形状可分为多边形和圆形。由于圆形变形体计算比较复杂，本节只介绍二维刚性体多边形计算模型。离散单元并不像有限单元必须符合三个基本方程，平衡方程、变形协调方程和本构方程，只要符合平衡方程即可，对于变形体还须符合本构方程。离散单元之间的接触
40、一般有三种方式：角角接触、边角接触或边边接触，如图10-17。对于图中块体B，其相邻块体分别对块体B通过接触点作用一组力，（=15），加上块体B的重力，它们对块体B的重心产生合力F和合力矩M，即 （10.64）图10-17 块体的集合及作用于个别块体上的力式中，yi为块体受力作用点坐标；xc，yc为块体重心坐标。如果上式中的合力和合力矩不等于零，那么块体B将按照牛顿第二定律的规律运动。令块体B的质量为m，转动惯量为I，块体绕其重心转动角度为，那么块体的运动方程为 （10.65）式中，u为块体位移；g为重力加速度。对上式可采用差分方法进行求解，对x方向的运动方程采用向前分格式进行数值积分，这样可
41、以得到岩块质心沿x方向的速度和位移 （10.66）式中，t0为初始时间；为计算时步。对于y方向的运动和转动，都有类似的算式。虽然离散单元可被视为不变形的刚体，但在单元接触的边界仍存在变形，如图10-18设块体之间的法向力Fn，两块体之间的“叠合”位移为un，两者之间成正比关系 （10.67）式中，Kn为接触面的法向刚度系数。 图10-18 离散单元间的作用力（a）块体之间压缩（b）块体之间剪切因为块体所受的剪切力与块体运动和加载历史无关，所以对于剪切力和剪切位移关系宜用增量形式表示 （10.68）式中，Ks为接触面的剪切刚度系数；为两块体之间的相对位移。 上面两式所表示的力和位移是弹性关系，它
42、的极限值符合Coulomb-Mohr准则，即 （10.69）式中，C和分别是接触面的黏结力和内摩擦角。当块体接触面上的剪切力大于极限值时，块体之间就会产生滑动。而当块体受到张力相互分离时，块体产生张拉破坏，作用在块体表面上的法向力和剪切力随之消失。块体的运动是一个不可逆过程，为了避免岩块在运动过程中发生振动，在离散元法计算中采用加阻尼的方法来耗散系统在震动过程中的动能。在静力平衡问题中，加阻尼来吸收系统的动能，以使系统达到稳定状态。阻尼模型一般可采用具有集中质量的Voigt模型，如图10-19，模型阻尼系统的自由振动微分方程为 （10.70）图10-19 具有集中质量的弹性阻尼系统式中，c为阻尼系数。根据振动理论，临界阻尼系数为。在实际计算时，阻尼系数取值应取略小于临界值，以使岩块运动没有回弹。引入阻尼并考虑岩块位移后，离散单元法的基本运动方程为 （10.71）式中，m为单元质量；c是黏性阻尼系数；k是为弹性刚度系数；f是单元所受的外荷载；u为位移；t是时间。上式可用动态松弛法进行求解。图10-20表示动态松弛法求解力和位移的计算循环。计算循环是以时步进行差分计算，在每一时步内进行一次迭代，根据前一次迭代所得到的块体位置，求出接触力作为下一次迭代的出发点，如