有限变形应变与裂隙网络渗透系数耦合模型

　　摘要　根据有限变形力学理论，编制了弹性均质岩层在采矿影响下位移和应变的有限元计算程序，并根据岩层位移可改变裂隙宽度，进而引起裂隙 渗透系数变化的思路，实现了有限变形应变与裂隙网络渗透系数的耦合分析，并就一个采矿引起顶板岩层裂隙渗透系数变化的模型进行了分析计算。结果表明，纵向裂隙的渗透系数增大约4倍，而在不考虑顶板岩层离层的情况下，横向裂隙的渗透系数减少约4倍。可见顶板岩层的变形对裂隙网络的渗透系数有显著的影响。
　　关键词　有限变形　裂隙介质　渗透率　耦合　模型
　　中国图书资料分类法分类号　P641.2

THE COUPLING MODEL OF LARGE DEFORMATION STRAIN
　　AND PERMEABILITY IN FISSURE NET

　　Li Yulin¹　Chen Zhida²　Hu Weiyue¹　Li Jingsheng¹
　　(Xi'an Branch,CCRI)¹　(Beijing Graduate School¹，China University of Mining and Technology)²

　　Abstract　On the basis of deformation mechanics,the paper promotes FEM program for the displacement and strain calculation of elastic,homogenous rock layers under mining,and realizes the coupling analysis of large deformation strain and fissure net permeability according to relations between displacement of rock layer and breadth of fissure.There is a model analysis and calculation for the change of fissure permeability of roof rock layers under mining.The results show that vertical fissure permeability increases about 4 times,and horizontal one decreases about 4 times under the condition of non-considering roof rock layers break.That is to say,the deformation in roof rock layer have a great influence on the fissure permeability.
　　Keywords　finite deformation;fisstured medium;permeability;coupling;models

1　引言
　　应力(应变)对岩体裂隙网络渗透系数的影响已越来越引起大家的注意，在水利水电坝基裂隙渗透性与扬压力的分析和计算及核废料处理领域已开展了许多有益的工作，但在煤炭开采研究领域相关工作较少。目前应力(应变)与裂隙网络渗透系数的分布参数耦合模型有以下几种类型；
　　a.　裂隙损伤系统分析模型^［1］；
　　b.　分别考虑岩块和裂隙的位移对岩体总应变贡献的应变与裂隙网络渗透系数耦合模型^［2］；
　　c.　不考虑岩块变形，用裂隙单元的分析方法直接求解裂隙位移实现应变与裂隙网络渗透系数的耦合模型^［3］。
　　上述模型中，作者认为第二种模型具有更好的适用性。值得指出的是，上述耦合模型中涉及岩体应力、应变及位移计算时，都采用的是小变形力学理论，这种耦合模型在岩体变形很小时是适用的，而采矿引起的顶板变形是典型的大变形现象。因而本文借鉴了Ouyang Z和Elsworth D^［2］的应变与裂隙网络渗透系数的耦合方法，建立了新的基于有限变形应变与裂隙网络渗透系数的耦合模型，并就一个采矿模型进行了分析和计算。

2　有限变形位移(应变)场数值计算的控制方程
　　有限变形问题的增量变分方程：

式中　——物性矩阵；
　　　——有限变形应变速率的物理分量；
　　——变形体中一点速度对实时拖带坐标的协变导数的物理分量；
　　　——减去渗透压力后作用于边界S_p上的等效外力速率；
　　　——重力速率；
　　　——求解域。
增量变分方程是数值解法的基础。
　　用实时物化拖带坐标系法推导获得的用相关矩阵表达的增量变分方程如下：

式中　D——弹性矩阵；
　　ΔU^l_k——节点k的位移增量；
　　　N——形函数；
　　ΔP_i——外力增量；
　　Δq_i——重力增量；
　　　B——形函数对拖带坐标物理分量的导数构成的形状矩阵；
　　　G——速度对初始拖带坐标系的协变导数与其对实时拖带坐标系协变导数的转换矩阵。

式中　Uⁱ｜_j——位移对初始拖带坐标系的协变导数，i,j=1,2；
　　　　xⁱ——拖带坐标的物理分量，i=1，2。
　　将研究区剖分为4节点等参数单元，用2×2点高斯积分求积，权系数W₁=1,W₂=1,所以有：

式中　　M——单元总数；
　　　　N₁——边界力作用单元数；
　　　　J——雅可比矩阵；
　　　　t——变形体厚度。

上式计算的结果为一组代数方程，表示为：

［K］［ΔU］=［R］，

式中　［K］——结构刚度矩阵；
　　　［ΔU］——节点位移增量列向量；
　　　［R］——扣除裂隙水压力后的有效外载荷向量。
解此方程组即可获得有限变形体中各节点的位移值，进而求得各单元的应变。需要指出的是，由于［K］是节点位移的函数，因而方程组是非线性的，这里选用Newton-Ramson迭代法求解。

3　裂隙网络渗透系数与应变的关系
　　设b为裂隙宽度，S为裂隙间距，μ_k为流体粘滞系数，g为重力加速度。对于应力作用状态下的裂隙岩体，变形后裂隙在x方向和y方向的位移分别为ΔU_jx，ΔU_jy，则有：

k_x，k_y表示裂隙网络在x和y方向的渗透系数主值。
　　岩体变形是由构成岩体的岩石变形和裂隙网络变形叠加产生的，因此　　

ΔU=ΔU_j＋ΔU_s，

ΔU为岩体中一点的位移，ΔU_j为裂隙产生的位移，ΔU_s为岩石产生的位移。
　　利用有限单元法求出的位移是岩体中一点的总位移，如果能够将裂隙岩体的变形合理地分配到裂隙和岩块上，也就是说，找到一个比例系数，进而求得裂隙变形的位移，就可以求得岩体变形后的渗透系数值。Ouyang Z和Elsworth D^［2］提出一个岩体模量折合参数的概念：

式中　E——构成岩体的岩石弹性模量；
　　　E_m——岩体的弹性模量；
　　　R_m——岩体模量折合系数。
经推导可以获得应变与渗透系数耦合的方程式^［2］：

当R_m=0时，说明岩体刚度很小，岩体的变形都集中在裂隙上，岩石本身的变形与裂隙相比非常小；当R_m=1时，表明岩体的弹性模量与岩石的弹性模量相等，岩体的变形均匀分配在裂隙与岩石中；当0＜R_m＜1时，说明岩体变形由裂隙变形和岩石变形共同组成，随R_m值不同，按不同比例分配，R_m取值如表1所示。

表1　岩体参数分类