《岩石力学的数值模拟》由会员分享，可在线阅读，更多相关《岩石力学的数值模拟（32页珍藏版）》请在金锄头文库上搜索。
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）四面体单元三角形的形函数N. = 2s （气 + b产 + % y）（i = i, j, m）式中，S为三角形面积；a = xy ； b. = y . - y ； c. = xx_。i i m i j k
6、 i k j四边形的形函数N. = 4（1+草）（1+n.n）（i = i,2,3,4）式中，位移量为u = %气；v =却匕；x = 吒气；四面体的形函数1N=:(a + b x + c y + d z)(i =i, j, k, m)i6V iiiix. y.z1yzjjjb=-1jja =x yzyzimmmimmx yz1yzk J kkkk1 xz1xyjjjJc =1 xzd=-1xyimmimm1 xz1xykkkJ k式中，V为四面体的体积。单元在直角坐标轴中位移分量分别为u，v，W，因此单元的位移矩阵为（10.1）, = u, v, w = N %式中，ue为单元位移矩阵；N为
7、形函数矩阵；5 e为单元节点位移列阵。根据几何方程，对位移矩阵求偏导数，可以得到应变矩阵 = B5（10.2）式中，B为连续单元节点位移和单元应变的矩阵，也称为应变矩阵。对于三角形单元，B 为常数矩阵，元素值取决于单元节点坐标差。根据本构方程，可以得到单元节点位移与单元应力矩阵之间的关系七 = D 气 = D B 5 （10.3）式中，D为弹性矩阵。应用虚功原理和最小势能原理可以推导出单元刚度矩阵的表达式Ke = f BTDBdv（10.4）各单元的体积力和面力按照静力的等效原则移置到各单元的节点上，其等效节点力为Pe = i N T Pdv Q = kNTQds（10.5）e v式中，Pe为
8、作用于单元体积力P的等效节点荷载；QJ为作用于单元面积力Q的等效 节点荷载:设环绕某节点，共有k个单元，则，节点上的外加荷载RJ为一 一 一W 一 一（10.6）R. = LR(e) + Z Q(e) + P i=1i=1式中，Pj为作用于i节点上集中力。将各单元节点力与节点位移之间的关系叠加，形成以节点位移列阵为基本未知量的线性 代数方程组（10.7）K = E Ke e=1K5 = R（10.8）求解上式有限个线性代数方程组，得到节点位移矩阵，根据相应的节点位移利用式（10.2） 和（10.3）计算单元的应力应变值。10.1.2非线性问题的处理方法从本质上讲，岩石力学问题都是非线性问题。这
9、是因为一方面岩石材料的应力应变本构 关系绝大多数呈现非线性特征，另一方面岩体的变形又大多是大变形。对于求解岩石力学的 非线性问题，解析方法显得无能为力，而有限元可以求解岩石力学绝大部分的非线性问题。 岩石力学的非线性问题可以分为三大类：材料非线性，即岩石材料的本构关系为非线性， 而变形的几何关系仍是线性的；几何非线性，即岩石几何变形为非线性，而本构关系仍为 线性；两类非线性问题的组合，即岩石材料既是材料本构非线性，又是几何非线性。这三 类非线性问题总体平衡方程组的共同特征都是非线性方程组，可表示为K （5 ）5 = 。（10.9）式中，K（5 ）为刚度矩阵，它是位移5 的函数。求解这类非线性方
10、程组一般有三种方 法：迭代法、增量法和混合法。迭代法又称为牛顿-拉斐逊法，对于一个变量的函数F = f （5 ），如图10-2，它的迭代过程如下：设函数值F由F0增加至Fa通过切线法做第i次近似值可由下式确定：K._1（A5 i） = FF._1（10.10）式中，K.1为第i-1次迭代时的曲线切线斜率，那么最终的解为5 i =5 i_1 +A5 i（10.11）牛顿-拉斐逊法的主要缺点是每次迭代过程中都要重新计算一下切线值，也就是刚度矩 阵及其逆矩阵，因此花费机时较长，为了避免每次计算K.值，每次迭代时都采用同一个初始的犬0,如图10-3，此方法称为修正的牛顿-拉斐逊法，具体的计算公式如下K
11、 0 i）a-Fi-1图10-2牛顿-拉斐逊迭代法图10-3修正的牛顿-拉斐逊法（10.12）两种迭代法比较而言，修正的牛顿-拉斐逊法的迭代次数多于牛顿-拉斐逊法，但它省去 了每次迭代时重新计算赐度矩阵K。计算时间的多少，不仅取决于迭代次数或收敛速度， 还取决于每次迭代所花费的时间，在一般情况下常刚度法在总体计算时间上比较合理，对于 非线性很强的材料，常刚度法并不适用。增量法也是把非线性问题线性化的一种处理方法。如图10-4把总荷载F划分成若干个 增量AF，逐级施加荷载增量进行求解，有限元计算公式为K,-1A8 J = AF.（10.13）那么总位移和总荷载为8 = 8 0 + E A8.i=
12、1F = + AF.（10.14）i=1增量法的误差较大，最终误差为各级增量时误差之和，为了减小误差，有时对一级增量 后，加上一个修正正项加以误差矫正，计算公式为Kt-1A8 . = AF. + -（10.15）十枷击占图10-4增量法由于增量法和迭代法各有其优缺点，目前有限元常常采用增量法和迭代法的结合，即 混合法，一般将非线性关系分成若干个荷载增量段，在每一个增量段内用迭代法求解逼近非 线性解，最终累加求得各级荷载作用下的最终应力与应变。10.1.3岩石弹塑性有限元分析12岩石弹塑性本构关系是岩石主要非线性问题之一。岩石的弹塑性是指岩石材料的应力应 变关系在屈服之前呈线性关系，当应力达到屈
13、服应力时，应力应变关系就变为非线性。由于 弹塑性模型中应变不仅依赖于受载的应力状态，而且与加载路径有关，因此一般弹塑性本构 关系不能用应力应变全量关系准确描述，只能用能反映加载路径有关的应力应变增量关系描 述。在岩石非线性弹塑性本构关系有限元分析中一般采用初应力法和初应变法求解非线性 平衡方程组。初应力法是将荷载以微小增量形式逐级加在模型上，每加一级荷载增量亦.，/就会产生相应的位移增量*、变形增量决和应力增量d。, 。对于具有初应力的弹 塑性应力应变本构关系可以写成d。 = Dp de + d。0（ 10.16）式中，d。 = -Dp&为初应力；D,为塑性矩阵，它与加载前的应力水平有关，而
14、与应力增量无关。初应力法是通过对d。的处理将应力修正到正确的水平上，初应力d。不仅与加 载增量前应力水平有关，还与本级所加荷载增量引起的变形增量de 有关，如图10-5。增 量形式的平衡方程为K0d5 = dF + dF（10.17）图10-5初应力法式中，K0为线弹性计算中的总刚度矩阵；dF为矫正荷载项，它由下式决定：dF = Ef BTD,de dV（10.18）由于dF随位移变化而变化，因此计算时必须进行迭代求解。初应力法求解按照下述计算 步骤实现：（1）把全部荷载划分成若干个增量，在每一级增量段内，按照增量弹塑性平衡方程进 行求解。（2）计算各单元的应力增量及当前的应用(10.19)dj = Bdd. j db = D d 如 i j b. = b. . 1 + db .1 I J I J-1I J式中，下标l表示第i级荷载增量；J表示第J次迭代。（3）根据
