岩石力学的数值模拟(讲义)

第10章岩石力学的数值模拟

随着计算机软硬件技术的迅速发展，使岩石力学有了长足的进步，特别在岩石力学的数值计算和模拟方面发展尤为迅速，使得许多岩石力学解析方法难于解决的问题得以重新认识。正如钱学森在给中国力学学会“力学——迎接21世纪新的挑战”的一封信中对力学发展趋势总结的那样“今日力学是一门用计算机计算去回答一切宏观的实际科学技术问题，计算方法非常重要”。岩石力学和其他力学学科一样，需要数值计算方法并推动岩石力学的发展。

岩石介质不同于金属材料，在数值计算方面具有其独特的特点[205]：

（1）岩石介质是赋存于地壳中的各向异性天然介质。

（2）岩石介质被众多的节理、裂缝等弱面所切割而呈现高度的非均质性，而其物理、化学及力学性质具有随机性特点。

（3）岩石介质赋存时以受压为主，而且抗压强度远大于抗拉强度。

（4）岩石力学与工程问题在时空分布上较广，从本质上讲都是三维问题。

（5）岩石工程一般无法进行原型试验，而实验室测得的数据不能直接应用于工程设计和计算。

（6）岩石力学与工程具有数据有限问题。

数值计算方法经过几十年的发展，目前已形成许多种岩石力学计算方法，主要有有限元法、边界元法、有限差分法、离散元法、流形元法、拉格朗日元法、不连续变形法及无单元法等。它们各有优缺点，有限元的理论基础和应用比较成熟，在金属材料和构件的计算中应用十分成功，但它是以连续介质为基础，似乎与岩体的非连续性有一定差距，流形元等数值方法虽然考虑了岩体中节理效应，但其理论基础还不完全成熟。相信在不久的将来，肯定会出现完全适合于岩体材料和工程的数值计算方法[206~208]。

10.1 岩石力学的有限元分析[209~213]

有限元法（finite element method，FEM）是岩石力学数值计算方法中最为广泛应用的一种。自20世纪50年代发展至今，有限元已成功地求解了许多复杂的岩石力学与工程问题。被广大岩石力学研究与工程技术人员喻为解决岩石工程问题的有效工具。有限元法是根据变分原理求解数学物理方程的一种数值方法。有限元法把连续体离散成有限个单元，每个单元的场函数只包含有限个节点参量的简单场函数，这些有限个单元的场函数集合构成整个结构连续体场函数。根据能量方程和加权函数方程可建立有限个求解参数的方程组，求解这些离散方程组，就是有限元法的精髓所在。虽然求解时把连续函数转化为求解有限个离散点处的函数值，但只要单元划分得充分小时，足可以满足计算要求。

有限元法求解问题时一般遵循以下步骤：

（1）有限元计算模型的建立，包括模型单元的划分、确定边界条件。

（2）对单元体进行力学分析，包括求解节点位移、单元应变和单元应力。

（3）对计算模型进行分析。

（4）进行计算分析。

10.1.1 线弹性有限元法的基本方程

线弹性有限元是弹塑性有限元、损伤有限元、流变有限元等非线性有限元的基础。线弹性有限元假定岩石介质连续、均质、小变形和完全弹性。

有限元法求解弹性力学问题时通常以位移作为基本未知量，单元位移是以单元节点位移为基本未知量，选择合理的位移插值函数，将单元位移表达为节点坐标的连续函数，插值函数也可称为形函数。不同形状的单元具有不同的形函数。

图10-1为三种最常见单元形式，即三角形、四边形及四面体单元。它们的形函数分别

为：

图10-1 有限元的三种基本单元形式

（a ）三角形单元（b ）四边形单元（c ）四面体单元

三角形的形函数

),,( )(21m j i i y c x b a S

N i i i i =++= 式中，S 为三角形面积；m i i y x a =；k j i y y b -=；j k i x x c =。

四边形的形函数

)4,3,2,1( )1)(1(4

1=++=i N i i i ηηξξ 式中，位移量为i i u N u ∑=；i i v N v ∑=；i i x N x ∑=；i i y N y ∑=。

四面体的形函数

),,,( )(61m k j i i z d y c x b a V

N i i i i i =+++= k k

k m m

m j j

j i z y x z y x z y x a = k k m m j j i z y z y z y b 111-= k k

m m j j i z x z x z x c 111= k k m m j

j i y x y x y x d 111-=

式中，V 为四面体的体积。

单元在直角坐标轴中位移分量分别为u ，v ，w ，因此单元的位移矩阵为

}]{[],,[}{e e e N w v u u δ== （10.1）

式中，}{e u 为单元位移矩阵；[N]为形函数矩阵；}{e δ为单元节点位移列阵。

根据几何方程，对位移矩阵求偏导数，可以得到应变矩阵

}]{[}{e e B δε= （10.2）

式中，[B]为连续单元节点位移和单元应变的矩阵，也称为应变矩阵。对于三角形单元，[B]为常数矩阵，元素值取决于单元节点坐标差。

根据本构方程，可以得到单元节点位移与单元应力矩阵之间的关系

}]{][[}]{[}{e e e B D D δεσ== （10.3）

式中，[D]为弹性矩阵。应用虚功原理和最小势能原理可以推导出单元刚度矩阵的表达式

⎰=v

T e dv B D B K ]][[][}{ （10.4） 各单元的体积力和面力按照静力的等效原则移置到各单元的节点上，其等效节点力为

⎪⎩⎪⎨⎧==⎰⎰v T e

v T e ds Q N Q dv P N P }{][}{}{][}{ （10.5） 式中，{P e }为作用于单元体积力{P }的等效节点荷载；{Q e }为作用于单元面积力{Q }的等效节点荷载。设环绕某节点i 共有k 个单元，则i 节点上的外加荷载{R i }为

i k

i e si k i e vi i P Q R R ++=∑∑==1)(1)(}{}{}{ （10.6） 式中，P i 为作用于i 节点上集中力。

将各单元节点力与节点位移之间的关系叠加，形成以节点位移列阵为基本未知量的线性代数方程组

∑==n

e e K K 1][][ （10.7）

}{}]{[R K =δ （10.8）

求解上式有限个线性代数方程组，得到节点位移矩阵，根据相应的节点位移利用式（10.2）和（10.3）计算单元的应力应变值。

10.1.2 非线性问题的处理方法

从本质上讲，岩石力学问题都是非线性问题。这是因为一方面岩石材料的应力应变本构关系绝大多数呈现非线性特征，另一方面岩体的变形又大多是大变形。对于求解岩石力学的非线性问题，解析方法显得无能为力，而有限元可以求解岩石力学绝大部分的非线性问题。岩石力学的非线性问题可以分为三大类：①材料非线性，即岩石材料的本构关系为非线性，而变形的几何关系仍是线性的；②几何非线性，即岩石几何变形为非线性，而本构关系仍为线性；③两类非线性问题的组合，即岩石材料既是材料本构非线性，又是几何非线性。这三类非线性问题总体平衡方程组的共同特征都是非线性方程组，可表示为

}{}})]{[({σδδ=K （10.9）

式中，})][({δK 为刚度矩阵，它是位移}{δ的函数。求解这类非线性方程组一般有三种方法：迭代法、增量法和混合法。

迭代法又称为牛顿-拉斐逊法，对于一个变量的函数)(δf F =，如图10-2，它的迭代过程如下：设函数值F 由F 0增加至F a 通过切线法做第i 次近似值可由下式确定：

11)(---=∆i a i i F F K δ （10.10）

式中，K i-1为第i -1次迭代时的曲线切线斜率，那么最终的解为

i i i δδδ∆+=-1 （10.11）

牛顿-拉斐逊法的主要缺点是每次迭代过程中都要重新计算一下切线值，也就是刚度矩阵及其逆矩阵，因此花费机时较长，为了避免每次计算K i 值，每次迭代时都采用同一个初