在采矿工程中, 经常需要对无限或半无限介质中的地下工程进行围岩应力和稳定性分析。实践证明, 有限元法解决诸如此类问题是非常有效的, 但是针对某一具体工程, 有限元法最致命的缺点是人为给定某种计算边界条件, 并需对分析的工程问题取足够大的区域进行离散, 以确保计算精度。对此监凯维奇(Zieckwicz)提出用边界积分方程与有限元耦合法来解决无限或半无限域问题, 这对提高计算精度, 减少数据存贮是很有益的。然而对于孔洞布置复杂, 尤其是对具有弹性以外的应力—应变行为的复合介质问题(三种以上介质), 使用有限元与边界元耦合法所表现的不足之处在于其本构关系的理论推导繁琐, 计算程序编制复杂, 对于某些非线性问题甚至无法实现。

一   无限元刚度矩阵的导出及耦合程序结构

无限单元是Bettess与Ungless于1977年先后提出来的, 其核心思想是利用位移的半解析解将坐标形函数与位移衰减函数相乘积来构造位移形函数, 通过坐标变换的关系将无限域变换成有限域, 从而实现距工程问题无限远处的位移衰减为零的真实边界条件。由此可见, 在无限元中衰减函数的选取是至关重要的, 衰减函数的种类很多, 但其共同的特点是要能反映无限远处位移衰减为零的变化规律, 同时确保无限单元刚度矩阵的广义积分收敛。本文针对静载问题选用

指数型衰减函数, 并以平面问题为例, 有限元取四节点等参元, 无限元取六节点非等参元。

1.无限单元刚度矩阵。无限单元如图 1所示。

图  1  六节点无限元

下载: 全尺寸图片 幻灯片

图 1中, ①、②为边界点, ③、④为衰减控制点, ⑤、⑥为无穷远点, 这样坐标xoy与ζoη之间的关系是:

其中Ni是坐标形函数，且为

由于无限单元域内位移与节点位移之间的关系为:

其中N_i为位移形函数, 取位移衰减函数

则位移形函数

依照弹性理论, 应变与位移之间的关系是:

将(5)式展开, 得无限单位的几何矩阵

其中|J|是雅可比行列式。

依照弹性理论, 利用泛函的驻值条件, 导出无限单元的刚度矩阵是

利用代换

则无限单元的刚度矩阵为

其中﹝D﹞是岩体的弹性矩阵。

2.节理无限元刚度矩阵。如图 2所示。

图  2  节理无限元

下载: 全尺寸图片 幻灯片

它是无厚度六节点节理无限元。节理无限元其模型概念是节理沿着节理方向伸向无限远处, 其理论推导是以节理上下盘的一对对应点的相对位移的法向分量ε_Д和切向分量ε_Б作为节理无限元的法向应变和切向应变, 并由此导出节理无限元的刚度矩阵, 其过程类似于普通单元, 限于篇幅, 直接给出其刚度矩阵。

其中

kn、ks分别是节理的法向刚度和切向刚度。

3.耦合程序结构。有限元与无限元耦合法可以方便的处理复合介质问题和非线性问题, 它克服了边界元法处理非均匀复合介质(三种以上介质)所存在的种种困难, 表现出良好的模拟性能，耦合程序结构见图 3。

图  3  耦合程序结构框图

下载: 全尺寸图片 幻灯片

二   边界条件的研究

为了研究有限元的人为零位移边界耦与合法的无限元边界对计算结果的影响, 我们以一条地处无限域的园形巷道为例, 计算使用的单元网格见图 4。通过改变网格边界的位置即取计算区域为1.5r, 2r, 2.5r, 4r, 6r, 其中区域向外扩展时, 内层单元划分的比例不变, 把在各种计算半径下有限元人为零位移边界与无限元边界的计算结果加以系统研究, 并以巷道表面节点的最大位移与解析解比较, 求其误差率, 同时以误差率的大小来定量描述边界条件对计算结果的影响。

图  4  研究边界影响的单元网络

下载: 全尺寸图片 幻灯片

误差率

式中:u⁰_max、是有限元人为零位移边界或无限元边界计算所得的最大垂直位移, u_max是解析解求得的最大位移。

对于无限域, 半径r=200厘米的巷道表面最大位移为:

式中:σ_v、为初应力的垂直分量；k₀为水平应力与垂直应力的比值；E是弹性模量；μ是泊松比, 用(10)式求得巷道表面的最大位移u_max=0.05厘米。

根据采矿工程的特性, 如果认为地下工程结构表面位移或围岩中的应力(数值分析值)与解析解之间的误差小于10%, 则这种计算半径和网格划分是合理的这一结论^﹝3﹞, 那么按照图 4所示的网格, 求得巷道表面位移的误差率与计算半径之间的关系见表 1。从表 1与图 5可以得出, 如果要使位移误差率小于10%, 采用有限元的零位移边界需要取6倍的巷道半径作为计算区域进行离散, 而采用无限元计算边界时仅需要取2.5倍的巷道半径作为计算区域, 这时的位移误差率为9.4%, 而应力误差却小于5%, 由此可见, 在相同的计算半径和同样的网格划分条件下, 边界条件对计算结果的影响很大。为了减少数据准备的工作量, 节省机时, 降低计算费用, 采用无限元计算边界是很有意义的。

表  1  误差率与计算半径之间的关系

下载: 导出CSV 

| 显示表格

图  5  误差率与计算半径的关系曲线

下载: 全尺寸图片 幻灯片

三   衰减系数的研究

仍以图 4的园形巷道为例, 取2.5倍的园形巷道半径进行有限元离散。这时用有限元人为零位移边界算得巷道表面最大位移u_max=0.0314厘米, 用解析法求得巷道表面的最大位移u_max=0.05厘米, 而用有限元与无限元耦合法算得的最大位移值却随衰减系数L值的不同而变化(见表 2)。从图 6可以看到, 不论L取何值, 耦合解都界于有限元解与解析解之间, 当L=0.5时, 耦合解的u_max=0.0318厘米, 它与有限元解0.0314厘米很接近, 当L不断增大时, u_max也不断增大, 并不断逼近解析值, 当L大于5时, 其值又疏远于解析解, 由此可见, 当L在0.5~1.5时, 耦合法的计算值与解析解的差值变化很大, 而当L大于2时, 这种差值的变化趋于缓和, 当L=4—5时, 计算精度最好。

表  2  巷道表面最大位移值与衰减系数的L关系  单位:厘米

下载: 导出CSV 

| 显示表格

图  6  巷道周边位移与衰减系数L的变化曲线

……analytical solusions; —infinite element solusions; —·—finite element solusions

下载: 全尺寸图片 幻灯片

四   衰减控制点的研究

仍以园形巷道为例, 这时无限元的衰减系数L=4.5, 以m=DI/DF的比值来研究衰减控制点对计算精度的影响。其中DI是衰减控制点到有限元边界的径向长度, DF是有限元最外层单元的径向长度, 并以2.5倍的圆形巷道进行有限元离散。通过计算可知, 巷道最大位移值随m值的变化而变化(见表 3)。从巷道周边最大位移值与m值的变化曲线(图 7)可见, 不论m取何值, 有限元与无限元耦合解所得的巷道周边最大位移值都界于有限元解与解析解之间。由此可见, 无限元的衰减控制点只要按照有限元最外层单元相同的比例布置, 就可以满足工程对计算精度的要求, 而当m=3时, 耦合法的计算精度最好。

图  7  巷道周边最大位移与m值的变化曲线

……analytical solusions; —infinite element solusions; —·—finite element solusions

下载: 全尺寸图片 幻灯片

表  3  巷道周边最大位移与m值的关系

下载: 导出CSV 

| 显示表格

五   非均匀应力场的初应力问题

在岩土工程的数值计算中, 首先需要确定初始应力场。目前的一般做法是在计算边界施加荷载或位移, 并同时考虑介质的自重。对于非均匀初应力场, 如果在计算区域的边界采用无限元, 则上述求初应力的方法就不适用了。首先是因为无限元主要是针对无限或半无限域由于岩体局部开挖所产生的围岩位移在距开挖岩体无限远衰减为零这一事实提出来的, 因此在无限元计算边界上不允许再附加任何位移和应力条件, 其次无限或半无限岩体在原岩应力场作用下产生的原岩位移不存在着无限远处原岩位移为零这一事实。按照原岩位移这种性质所选取的任何位移形函数都将造成无限单元刚度矩阵的广义积分发散。目前国内外已有许多文献把非均匀初应力问题看成是无限元急待解决的一个课题。

众所周知, 在地下工程问题的数值计算中, 人们所关心的是在初应力场作用下, 由于开挖所引起的一系列力学效应, 对开挖前原岩本身的位移是不感兴趣的。通常的计算方法是先求出原岩应力场和原岩位移, 然后在计算开挖效应时, 首先保留原岩应力, 并把原岩位移充零, 以计算出由于开挖而在岩体中产生的次生应力和位移。基于这一事实, 我们可以采用回避无限元目前无法处理非均匀应力场的初应力问题, 采用有限元零位移边界计算初应力值。为此在编制耦合程序时只需实现人为零位移边界计算初应力场, 同时把计算区域外围的岩体自重转化成作用在边界上的荷载加以考虑, 当计算开挖效应时, 采用无限元计算边界, 这时不再考虑无限域岩体的自重。从而可以简便、顺利地解决非均匀应力场的初应力问题。计算结果表明, 这种方法是非常有效的。

六   结论

1.由于引入无限元边界, 它实现了无限远处位移为零的真实边界, 从而使有限元的模拟性能更好、更完善。

2.无限元概念清楚, 简便易行, 使有限元与无限元耦合法在程序编制上十分便利, 它克服了边界元法处理复合介质问题所存在的各种困难。

3.在采矿工程中, 人们只对开挖岩体附近很小一个区域的应力和位移感兴趣, 这样在使用耦合法时, 只需在开挖面附近布置很少几层单元, 就能获得满意的计算精度, 从而大大减少了计算准备的工作量, 降低了计算费用。

4.无限元的核心是选取衰减函数, 通过对指数型衰减系数的研究, 可以根据工程特性选取衰减系数, 从而使衰减函数的选取规范化。

5.无限元衰减控制点对计算精度的影响很大, 通过对m值的研究, 只要使m=1即衰减控制点与相邻的有限元按照相同的比例布置, 就能保证计算精度满足工程要求。

6.采用有限元的人为零位移边界确定初使应力场, 用无限元处理开挖时的有限元计算边界, 并将计算区域外围岩体的自重转化成边界载荷, 这是使用有限元与无限元耦合法处理非均匀应力场初应力问题的一条十分有效的途径。