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

第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 值，每次迭代时都采用同一个初
始的K 0，如图10-3，此方法称为修正的牛顿-拉斐逊法，具体的计算公式如下
10)(--=∆i a i F F K δ （10.12）
图10-2 牛顿-拉斐逊迭代法 图10-3 修正的牛顿-拉斐逊法
两种迭代法比较而言，修正的牛顿-拉斐逊法的迭代次数多于牛顿-拉斐逊法，但它省去了每次迭代时重新计算赐度矩阵K i 。

计算时间的多少，不仅取决于迭代次数或收敛速度，还取决于每次迭代所花费的时间，在一般情况下常刚度法在总体计算时间上比较合理，对于非线性很强的材料，常刚度法并不适用。

增量法也是把非线性问题线性化的一种处理方法。

如图10-4把总荷载F 划分成若干个增量F ∆，逐级施加荷载增量进行求解，有限元计算公式为
}{}]{[1i i i F K ∆=∆-δ （10.13）
那么总位移和总荷载为
∑=∆+=n
i i 10}{}{}{δδδ
∑=∆+=n
i i F F F 10}{}{}{ （10.14）
增量法的误差较大，最终误差为各级增量时误差之和，为了减小误差，有时对一级增量后，加上一个修正正项加以误差矫正，计算公式为
}{}{}]{[11--'∆+∆=∆i i i i F F K δ （10.15）
图10-4 增量法
由于增量法和迭代法各有其优缺点， 目前有限元常常采用增量法和迭代法的结合，即混合法，一般将非线性关系分成若干个荷载增量段，在每一个增量段内用迭代法求解逼近非线性解，最终累加求得各级荷载作用下的最终应力与应变。

10.1.3 岩石弹塑性有限元分析[12]
岩石弹塑性本构关系是岩石主要非线性问题之一。

岩石的弹塑性是指岩石材料的应力应变关系在屈服之前呈线性关系，当应力达到屈服应力时，应力应变关系就变为非线性。

由于弹塑性模型中应变不仅依赖于受载的应力状态，而且与加载路径有关，因此一般弹塑性本构关系不能用应力应变全量关系准确描述，只能用能反映加载路径有关的应力应变增量关系描述。

在岩石非线性弹塑性本构关系有限元分析中一般采用初应力法和初应变法求解非线性平衡方程组。

初应力法是将荷载以微小增量形式逐级加在模型上，每加一级荷载增量}{i dF ，就会产生相应的位移增量}{i d δ、变形增量}{i d ε和应力增量}{i d σ。

对于具有初应力的弹塑性应力应变本构关系可以写成
}{}]{[}{0σεσd d D d p += （10.16）
式中，}]{[}{εσd D d p -=为初应力；][p D 为塑性矩阵，它与加载前的应力水平有关，而与应力增量无关。

初应力法是通过对}{0σd 的处理将应力修正到正确的水平上，初应力}{0σd 不仅与加载增量前应力水平有关，还与本级所加荷载增量引起的变形增量}{εd 有关，如图10-5。

增量形式的平衡方程为
}{}{}]{[0F d dF d K '+=δ （10.17）
图10-5 初应力法
式中，[K 0]为线弹性计算中的总刚度矩阵；}{F d '为矫正荷载项，它由下式决定：
∑⎰='dV d D B F d p T }{][][}{ε （10.18）
由于}{F d '随位移变化而变化，因此计算时必须进行迭代求解。

初应力法求解按照下述计算步骤实现：
（1）把全部荷载划分成若干个增量，在每一级增量段内，按照增量弹塑性平衡方程进行求解。

（2）计算各单元的应力增量及当前的应用
⎪⎩⎪⎨⎧+===-j i j i j
i j i j i j i j i d d D d d B d }{}{}{}]{[}{}]{[}{1σσσεσδε （10.19） 式中，下标i 表示第i 级荷载增量；j 表示第j 次迭代。

（3）根据岩石的屈服准则，由各单元应力判断单元是否屈服，对于塑性单元，计算应力修正项并修正应力
⎩⎨⎧-='=j p j i j
i j i p j p d d D d }{}{}{}]{[}{σσσεσ （10.20） （4）塑性单元通过修正项j p d }{σ计算等效节点力，所有塑性单元的等效节点力叠加构成总的修正荷载矢量
∑='dV d B F d j p T j i }{][}{σ （10.21）
（5）在修正荷载作用下进行下次迭代运算，此时基本方程为
j i j i dF d K }{}]{[=δ （10.22）
重复进行（2）～（5）步计算直至所有的塑性单元达到收敛精度要求。

再进行下一步的荷载增量计算。

（6）重新施加下一级荷载增量}{1+i dF ，重复计算（1）～（5）步，直至计算完毕。

通过累加各级荷载作用的计算结果，求得总位移}{u 和总应力}{σ。

一般初应力法的收敛速度比较缓慢，因此通常采用常刚度和变刚度法相结合的方法加速收敛。

初应变法按照弹塑性增量理论，总的应变量可表示为
p e d d d }{}{}{εεε+= （10.23）
式中，e d }{ε为弹性应变增量；p d }{ε为塑性应变增量，类似初应力法可以把塑性应变增量p d }{ε看做初应变。

因为在计算过程中p d }{ε的计算格式和弹性系统中的初应变}{0εd 一致。

图10-6 初应变法
从图10-6中可以看出，荷载增量dF 所对应的位移为B du ，按线弹性计算为A du ，两者的位移增量之差称为初始位移，它是由材料塑性引起的附加位移。

与模型系统初始位移对应的单元位移为0δd ，那么单元的初应变为
}]{[}{00δεd B d = （10.24）
10.1.4 岩石节理单元的有限元方法
岩体中常存在大量节理，而节理的变形性质和岩体力学性质有十分明显差别。

从某种程度上讲，节理的存在决定着岩体的力学性质。

因此分析节理的变形性质，对于研究岩体工程的变形情况至关重要。

在进行有限元分析时，这种节理一般采用专门的节理单元来处理。

节理单元首先由Goodman 提出[19，214]，Goodman 认为节理不是一个实体，它只在若干
点处接触，因此节理单元也应满足这一特点，图10-7为无厚度节理单元，节点1与4，2与3具有相同的坐标，沿节理面的应力分别为法向应力n σ和剪切应力s τ。

把节理单元两侧对应的位移定义为应变，相对的法向位移v ∆称为法向应变n ε，相对的切向位移称为切向应变s ε，那么节理单元的本构关系为
}{][}{εσ'=D （10.25）
图10-7 无厚度节理单元
式中，]['D 为单元的刚度矩阵，对于节理单元长度上任何一点s 处的应变可定义为界面两侧相应位移之差
)()(12314v v L s v v L s n -+
-⎪⎭⎫ ⎝⎛-=ε （10.26） )()(12314u u L s u u L s s -+
-⎪⎭⎫ ⎝⎛
-=ε （10.27） 上式可用矩阵形式表示为
}]{[}{δε'=B
根据一般化公式导出节理单元对应于局部坐标的刚度矩阵
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------='n n n n s s s s n n n n s s s s n n n n s s s s n n n n s s s s e K K K K K K K K K K K K K K K K K K K K K K K K K K K K K K K K K 200020020002020200
0020200202000020202000200
20002][ （10.28） 式中，n K ，s K 分别为法向和切向刚度。

对于整体坐标的刚度矩阵通过坐标变换矩阵得到，具体如下
]][[][][T K T K e
T e '= （10.29） 式中，[T ]为变换矩阵，它具人分块形式的对角阵
]][],[],[],[[][0000T T T T T = （10.30）
式中，⎥⎦
⎤⎢⎣⎡-=θθθθcos sin sin cos ][0T 。

这种无厚度节理单元适用于模拟直接接触的界面，而对于有一定厚度的弱面或夹层不适宜。

10.2 岩石力学的边界元分析
边界元法（boundary element method, BEM ）是和有限元法同步发展的一种数值计算方法。

与有限元相比有以下特点：①边界元法把一个均质区域看做一个大单元，只把它的边界离散化，区域内不划分单元，场变量处处连续；②对于无限区域，场变量自动满足无穷边界条件及自然表面状态。

对于半无限区域，场变量也自动满足无穷远边界条件及自然表面状态。

有限元法是全区域离散化，而边界元是把基本方程转化为边界积分方程，只对边界离散化建立相应的方程组进行求解。

这样边界元使三维问题降为二维问题求解，使二维问题转化为一维问题求解，当物体的表面积和体积之比较小时，边界元的划分单元数要比有限元少数倍和十几倍，这样也使待解的方程数目、处理和存储的数据量降低同样的倍数，大大节省了机时和算题费用。

当仅需求解物体内部几个点的应力时，有限元仍不得不划分整个区域，才能确定这几个点的应力值。

而边界元当知道边界的应力解时，就可以根据需要去求物体内部个别点的解。

在应力梯度较高处，有限元法的剖分密度常常受到限制，而边界元可以方便地确定应力梯度的分布。

但随着计算机硬件的飞速发展，边界元的优势逐渐显得不很明显。

边界元法矩阵方程中系数阵的元素结构比有限元法刚度阵中的元素复杂。

有限元刚度阵属带状稀疏阵，而边界元法的系数阵为满阵，因此对于面积和体积之比较大的薄壁结构而言，边界元的优越性就不明显。

边界元比较适合求解无限区域和半无限区域问题，如深埋巷道是一个典型的例子。

但边界元在计算非均质介质问题、非线性问题以及模拟工程开挖过程等方面不如有限元方便有效。

有限元与边界元划分单元的比较如图10-8。

图10-8 有限元法与边界元法比较
（a ）力学模型和边界条件（b ）有限元单元划分（c ）边界元单元划分
求解边界方程组主要有Massonet 和Rizze 分别提出的直接和间接两种数值解法。

直接法是直接建立关于边界的积分方程，通过离散化求解边界未知参数，并进而求解计算区域内场函数值。

间接法是设立一个在求解区域内含有若干系数的满足基本方程组的解，代入边界条件求出系数，进而求得边界及区域内各点的场函数值。

10.2.1 边界元分析原理[215~217]
在无限的弹性体内作用一单位力引起周围的应力和位移的解析解是边界元求解弹性体问题的基础，如图10-9在平面无限体内i 点作用一x 方向的单位力，其基本的弹性解在1848年由Kelvin 解出
⎪⎭
⎪⎬⎫⎪⎩⎪⎨⎧⎪⎪⎭⎫
⎝⎛∂∂-∂∂--⎥
⎦⎤⎢⎣⎡∂∂∂∂+∆-∂∂--=l k k l k l lk lk n x r n x r v x r x r v n r r v )21()21()1(41*πσ （10.31） ⎥⎦
⎤⎢⎣⎡∂∂∂∂+∆--=k l lk lk x r x r r v v G u 1ln )43()1(81*π （10.32） 式中，*lk σ，*lk u 为i 点l 方向单位力引起j 点k 方向的应力和位移；lk ∆为Kronecker 符号，当l =k 时，1=∆lk ，当k l ≠时，0=∆lk ；r 为i ，j 两点之间的距离（矢径）；k n ，l n 为j 点作用面外法向对于k 和l 的方向余弦；n
r ∂∂为矢径r 对外法向n 的方向倒数；r x x r l l =∂∂；r
x x r k k =∂∂。

当不考虑体力时，根据功的互等定理和以上的Kelvin 基本解，可以建立弹性体边界上i ，j 两点之间的应力和位移关系（如图10-10）
**u u u c i i σσ=+ （10.33）
图10-9 开尔文基本解条件 图10-10 边界积分方程的建立
式中，i c 的大小取决于边界情况，当边界为光滑曲线时，i
c 等于0.5；当边界曲线不光滑时，i c 的值根据刚体位移原则求解。

当边界Γ作用分布力时，j 点绕行一周，按叠加原理积分得
Γ=Γ+⎰⎰Γ
Γd u d u u c i i **σσ （10.34） 式（10.34）就是边界元中的边界积分方程。

10.2.2 边界单元与边界的离散
边界元是通过离散边界来求解边界积分方程。

对于一个确定区域的边界Γ进行离散，划分成有限个线段，称为边界单元。

根据单元的节点数目和节点模式，边界单元可分为常量单元、线性单元、二次单元等。

常量单元及线性单元均为直线段，常量单元以单元的中点为节点，每个单元有一个节点，场变量在单元内是一个常数。

而线性单元的场变量在单元内呈线性变化。

单元的模式与有限元相似，可以表示成插值函数形式，即
∑==m i i i u N u 1， ∑==m
i i i N 1σσ （10.35）
式中，m 为单元的节点数目；插值函数i N 由拉格朗日插值公式给出；i
u 为节点i 处的值。

如图10-11几个常见的插值函数如下：
常量单元 1=m ， 1=i N
线性单元
2=m ， )1(211ξ-=
N ， )1(2
12ξ+=N 二次单元 3=m ， ξξ)1(211-=N ， )1)(1(212ξξ+-=N ， ξξ)1(2
13+=N
图10-11 常见的几种边界单元
（a ）常量单元（b ）线性单元（c ）二次单元
对于常量的边界元，边界积分方程可简化为
⎰∑⎰∑∆Γ∆ΓΓ=Γ+d u d u u c n
n i i *1*1σσ （10.36）
对于平面问题，上式有2n 个方程，可写成矩阵形式
}]{[}]{[σG u H = （10.37）
式中，}{u ，}{σ分别为常量边界单元中点的位移列阵和应力列阵；[G ]，[H ]为系数矩阵。

边界是元法的边界条件一般都为混合边界条件，即在一部分边界上位移和应力已知，而另一部分未知。

为了方程求解，把所有的已知量移到等式的右边，未知量移到等式右边，经过变换后式（10.37）可简化为
}{}]{[F X A = （10.38）
式中，{X }为包含所有位移和应力未知量列阵；{F }为已知量列阵；[A ]为系数矩阵。

求解此方程可求得全部边界节点上未知量，由此进而求得计算域内任意一点的位移和应力值。

如图10-12，根据开尔文基本解和Betti 的互等定理，可以得到计算域内任意一点的位移和边界点外力功的关系式
⎰⎰Γ
ΓΓ=Γ+d u ud u i σσ** （10.39）
图10-12 计算域内i 点与边界j 点作用力与位移的关系图
边界离散后，对于常量单元，上式可改写为
⎰∑⎰∑ΓΓΓ=Γ+d u d u u n
n i *1*1σσ （10.40）
欲求在i 点l 方向的位移分量，可建立边界积分方程
⎰⎰Γ
ΓΓ=Γ+d u d u u k lk k lk i σσ** （10.41） 式中，*lk σ，*lk u 分别为i 点l 方向作用单位力于j 点k 方向引起的应力位移开尔文基本解。

欲求计算域内i 点的应力可由几何方程和Lame 方程求得
⎪⎪⎭
⎫ ⎝
⎛∂∂+∂∂+∂∂∆-=i j j i i i j i ij x u x u G x u v Gv 212σ （10.42） 将边界积分方程代入上式可得 ⎰⎰ΓΓΓ-Γ=d u S d D k kij k kij ij σσ （10.43）
对于二维问题，上式中系数D ，S 分别为
},,,2],,,)21{[()1(41
k j i k ij i kj j ki kij r r r r r r v r
v D +∆+∆-∆--=π （10.44）}
)41()
,,2)(21(),,,,(2],,,4),(,)21[(2{)1(422ij k lk i lk j j i k k i j k j i k j i j ij k ij kij n v n n r r n v r r n r r n v r r r r v r v n r r v G S ∆--∆+∆+-+++-∆+∆-∂∂-=π （10.45） 式中，i i x r r ∂∂=,，j
j x r r ∂∂=,，k k x r r ∂∂=,。

10.2.3 边界元与有限元耦合
如前所述，边界元和有限元在计算时各有优缺点，为了发挥各自的优点，提高求解精度和解题效率，A.Wexler 于1972年提出了边界元和有限元耦合求解的思路和方法。

以后，J.R.Osias 和M.Margulies 等对边界元和有限元的耦合求解进行了深入的研究。

为了讨论方便，把有限元的基本方程改写为如下形式
}{}{}]{[D F U K += （10.46）
式中，F 为边界力的等效节点力；D 为体力的等效节点力。

由于边界力可以由面力P 与面积之积获得，那么有限元方程可改写为
}{}]{[}]{[D P M U K += （10.47）
一般边界元方程在考虑体力的情况下也可写为
}{}]{[}]{[B P G U H += （10.48）
从式（10.47）和（10.48）可以看出，有限元方程和边界元方程具有类似的形式，因此从解题方式上提供了建立两者耦合的基本条件。

把一个计算域人为地分成两个子域，对两个不同的子域分别采用边界元和有限元进行求解。

如图10-13有限元的边界为S 2，边界元的边界为S 1，两者的共同边界为S 3。

对两个区域分别建立有限元方程和边界元方程，再利用两者边界S 3上的平衡和协调条件建立耦合方程组。

耦合方程组有两种方式建立：①把边界元方程转化为等价的有限元方程；②把有限元方程转化为等价的边界元方程。

一般在耦合求解时前者采用得较多。

图10-13 边界元和有限元耦合的区域划分
边界元方程转化为有限元方程时，首先对边界元方程进行改造，对式（10.48）中G 进行求逆可得
}{}){}]{([][1P B u H G =-- （10.49）
因为边界元法中{P }为分布力，要转化为有限元中的等效节点力，必须在上式两侧左乘一个
[M ]矩阵，那么可得到与有限元形式一致的方程
}{}{}]{[D F U K '+'=' （10.50）
式中
][][]][[1K H G M '=-；}{}{]][[1D B G M '=-；}{}]{[F P M '= （10.51）
需要注意的是有限元法的刚度矩阵为对称矩阵，而由边界元方程转化来的][K '为非对称，为了完全耦合，必须通过最小二乘法使其对称化。

对称化后的刚度系数为
)(21ji ij ij k k k '+'= （10.52） 一个建筑结构在地基上沉降问题比较适合耦合求解，如图10-14，把建筑结构部分划分为有限元，地基部分为边界元。

边界元的离散化有以下两种形式，由地面延伸到到无限远处，边界元的划分须近似截取一个有限区，或利用半无限平面问题的基本解引入到无限边界元。

本算例采用耦合成有限元方程进行求解。

计算结果和有限元结果比较如表10-1。

可以看出耦合法和有限元计算结果十分一致，但计算时间大为减少。

有限元解
-339 -97 135 361 600 耦合法解 -350 -105 135 370 617
图10-14 有限元和边界元耦合求解的算例[12]
10.3 岩石力学有限差分方法（FLAC ）
有限差分方法（finite difference method ，FDM ）是从一般的物理现象出发建立相应的微分方程，经离散后得到差分方程，再进行求解的方法。

这种方法可能是最古老的求解方程组的数值方法之一。

差分方程在计算机出现以前用一般的手摇计算器也可以求解。

随着计算机的不断发展和其他计算方法的兴起，有限差分法曾一度受到冷遇，但20世纪80年代末由美国ITASCA 公司开发的FLAC （fast Lagrangian analysis of continua ）程序广泛采用差分方法进行求解，在岩土工程数值计算中得到了广泛的应用。

10.3.1 FLAC 程序的简介[218]
FLAC 是为岩土工程应用而开发的连续介质显式有限差分计算机软件，主要适用于模拟计算岩土类工程地质材料的力学行为，特别是岩土材料达到屈服极限后产生的塑性流动。

材料通过单元和区域表示，根据研究对象的形状构成相应的网络结构。

每个单元在外载和边界
约束条件作用下，按照约定的线性和非线性应力，应变关系产生力学响应。

FLAC软件建立在拉格朗日算法基础上，特别适用于模拟材料的大变形和扭曲转动。

FLAC程序设有多种本构模型，可模拟地质材料的高度非线性（包括应变软化和硬化）、不可逆剪切破坏和压密、黏弹性（蠕变）、孔隙介质的流固耦合、热力耦合以及动力学行为等。

另外，程序还设有边界单元，可以模拟断层、节理和摩擦边界的滑动、张开和闭合行为。

支护结构，如衬砌、锚杆、可缩性支架或板壳等与围岩的相互作用也可以在FLAC程序中进行模拟。

同时用户可根据需要在FLAC中创建自己的本构模型，进行各种特殊的修正和补充。

FLAC程序具有强大的后处理功能，用户可以直接在屏幕上绘制或以文件形式创建和输出多种形式的图形，使用者还可以根据需要，将若干个变量合并在同一幅图形中进行研究分析。

由于FLAC程序主要为地质工程应用而开发的岩石力学数值计算程序，包括反映地质材料力学效应的特殊要求。

FLAC设计有7种材料本构模型：
（1）各向同性弹性材料模型。

（2）横观各向同性弹性材料模型。

（3）Coulomb-Mohr弹塑性材料模型。

（4）应变软化、硬化塑性材料模型。

（5）双屈服塑性材料模型。

（6）节理材料模型。

（7）空单元模型，可用来模拟地下开挖和煤层开采。

FLAC程序采用的是快速拉格朗日方法，基于显式差分来获得模型的全部运动方程（包括内变量）的时间步长解。

程序将计算模型划分为若干个不同形状的三维单元，单元之间用节点相互连接。

对某一个节点施加荷载之后，该节点的运动方程可以写成时间步长的有限差分形式。

在某一个微小的时间内，作用于该点的荷载只对周围的若干节点（相邻节点）有影响。

根据单元节点的速度变化和时间，程序可求出单元之间的相对位移，进而可以求出单元应变；根据单元材料的本构方程可求出单元应力。

随着时间的推移，这一过程将扩展到整个计算范围，直到边界。

这样程序可以追踪模型从渐进破坏直至整体破坏的全过程。

FLAC程序将计算单元之间的不平衡力，并将此不平衡力重新加到各节点上，再进行下一步的迭代运算，直到不平衡力足够小或者各节点的位移趋于平衡为止。

图10-15为FLAC 程序的计算循环示意图。

FLAC和有限元相比具有以下一些特点：
（1）由Maxti和Cundall于1982年提出的混合离散法适用于塑性破坏荷载和塑性流动的模拟，这个方法比有限元中普遍采用的归约积分法更为合理。

图10-15 FLAC程序中的计算循环示意图
（2）FLAC虽然具有动力方程，但模拟系统实质上是静止的。

这使得FLAC不需要数值上卸载就可以遵循物理的失稳过程。