王 庆，姚 俊，谭文禄

(中电建生态环境集团有限公司，广东 深圳 518102)

在边坡[1-2]、坝基[3]、隧洞[4]等岩石工程中，经常会遇到地下水与岩石(体)之间的耦合问题，其作用机理是：一方面，当岩石(体)变形时会引起孔隙率、节理厚度等变化，从而引起渗流场的变化；另一方面渗流场变化会引起孔隙水压力和有效应力的变化，从而引起岩石(体)应力应变场的变化。对应力渗流耦合最早的研究是Biot的三维固结理论[5]，但该理论中渗透系数在变形过程中为常数。有关学者的研究表明[6-9]，渗透系数在岩石的变形破坏过程中，会发生较大的变化，可能会影响工程的安全性，法国Malpasset拱坝失事正是没有考虑到渗透系数变化而导致的[10]，因此有必要研究岩石变形过程中渗透系数的变化规律。

目前应力渗流耦合的研究方法主要分为宏观唯象和细观2类方法。宏观唯象方法又可以分为直接法和间接法。直接法是指通过室内室外实验对实验数据进行拟合得到经验公式的方法[11-12]，这些公式只能反映岩石变形过程中渗透系数的变化的某一阶段，不能体现应力峰值后的变化规律。为了克服经验公式的这一缺点，学者们开始研究全应力-应变过程中的渗透系数变化规律[13-16]，但由于公式的参数太多，形式受岩石性质、实验条件等影响太大，而且不能解释渗透系数变化的物理力学机理，因此不利于在工程中应用。间接法则是通过引入中间变量，间接建立应力(应变)与渗透系数之间的关系，最常用的中间变量是孔隙率[17]。间接法虽然能减少模型参数数量，但岩石应力应变全过程不同阶段的渗透系数变化的机理不同，如峰前渗透系数急剧增大的主要原因是微裂纹扩展，而不是孔隙率变化，因此用一个中间变量很难描述全过程的渗透系数变化规律。细观方法是通过描述岩石在变形过程中细观结构(孔隙、微裂纹等)的变化，从而引起渗透系数变化的方法[18-20]。细观方法可以解释岩石变形过程中的渗透系数变化机理，但计算过程中需要记录每一个微裂纹的倾向、倾角、开度等参数，因此需要借助计算机程序实现。

为了提出既能解释渗透系数的变化机理，又方便应用的模型，许多学者开始从等效连续介质的角度将岩石细观结构的变化等效为损伤演化[21-23]。将损伤作为连接应力(应变)和渗透系数之间的桥梁，建立岩石、混凝土、水泥砂浆等类岩石脆性材料的弹性损伤-渗流耦合模型，这个耦合过程中有2个重要的准则，一个是损伤的演化准则，另一个则是损伤过程中的渗透系数变化准则，本文在传统的弹性损伤模型和渗透性变化模型基础上，提出改进的弹性损伤-渗流耦合模型。

1 改进的弹性损伤模型

传统的岩石弹塑性损伤模型，虽然有可以考虑塑性和可以区分拉压破坏的优点，但是该本构关系也有一些缺点。首先是输入参数过多，弹塑性损伤本构关系中拉压曲线是分开的，而且要输入塑性强化的应力～应变曲线以及损伤因子随着开裂应变的变化曲线，这些曲线的参数过多；另外一点是使用弹塑性损伤本构关系计算，结果很难收敛，而且收敛性也不能控制。鉴于这两个缺点，这里在Mazars模型的基础上，改进其收敛性和网格依赖性，提出一种改进的弹性损伤模型。Mazars模型是一种描述混凝土损伤的本构关系，在该模型中，应力应变关系如下：

(1)

(2)

(3)

(4)

其中，ABS()是取绝对值函数，如果f为正，结果为f；如果为负，结果为0。

在Mazars模型中，等效应变作为内变量来判断损伤的初始化和演化，但是这里是将有效应力分解为拉压部分，按照这种思想，将拉压有效应力的等效应力作为损伤初始化准则，形式如下：

(5)

(6)

(7)

(8)

(9)

式中：

(10)

(11)

(12)

(13)

式(6)、(8)是受拉损伤的初始化准则，其形式和最大拉应力准则类似，式(7)、(9)是受压损伤的初始化准则，由于受压破坏中含有压缩和剪切破坏的组合，其形式类似于Drucker-Prager屈服准则，K是Drucker-Prager屈服准则的锥角，R0是双轴压缩强度和单轴压缩强度比，f±是拉伸和压缩强度。

Mazars模型中的损伤演化准则形式如下：

(14)

(15)

材料的软化阶段会导致边值问题的不适定，式(14)、(15)中的损伤演化准则对网格的依赖性很大，为了减小网格依赖性，这里采用断裂能的手段[25]来修正损伤演化准则。常数A+和B-可以写成如下断裂能的形式：

(16)

(17)

(18)

式中G±——拉压断裂能；lc——微裂纹的特征长度；Ve——单元体积；E0——初始弹性模量。

考虑拉损伤和压损伤的总的损伤因子定义如下：

d=1-(1-d+)(1-d-)

(19)

以上的改进策略可以减小网格依赖性，但是材料软化在隐式分析中，如ABAQUS Standard, 还存在收敛性问题，这里采用正则化准则来处理收敛性问题。正则化准则的主要目的是使软化材料的切线模量在增量步足够小的时候保持正定。根据式(1)中的应力-应变关系，可以得到Jacobian矩阵的形式如下:

(20)

为了提高收敛性，ABAQUS子程序UMAT中使用一种损伤因子的黏性正则化准则，此时，不直接使用式(14)、(15)来计算损伤因子，而是使用如下形式的正则化损伤因子：

(21)

(22)

(23)

(24)

从式(23)和式(24)中可以看出：

(25)

因此，使用正则化损伤因子计算得到的Jacobian矩阵形式如下：

(26)

需要注意黏性系数η的选择，过小的黏性系数导致收敛所需要的增量步很小；过大的黏性系数可能会过度延迟刚度的弱化，导致软化阶段的应力-应变曲线偏离实际过多。

2 损伤过程中渗透系数的变化

岩石的压缩过程中渗透性试验中反应出渗透系数的变化通常具有几个阶段。根据Zhang[26]的试验数据描绘的典型渗透系数变化的几个阶段见图1。A是应力-应变曲线中损伤开始的点，对应于渗透系数开始快速增大的点，B点是应力-应变曲线中的峰值应力点，对应于渗透系数曲线的转折点，B点之后渗透系数增长速度有所降低，然后达到渗透系数最大的点C，对应于应力-应变曲线的转折点，C点之后应力下降速度减缓。OA阶段渗透系数变化不大，但是有时候在开始的时候渗透系数有所减小，原因是在压应力作用下孔隙变化和初始微裂纹的闭合。AC阶段渗透系数的增大主要是因为微裂纹数目的增多以及微裂纹的张开，C点之后渗透系数的减小主要是因为试件破坏之后在围压的作用下宏观裂纹闭合造成的。

图1 脆性岩石的典型应力-应变曲线和渗透系数-应变曲线

目前描述脆性岩石压缩过程中渗透系数的变化的经验模型和理论模型只能描述渗透系数变化的某个阶段，通常是增大的阶段，并不能同时描述渗透系数变化的整个过程。假设岩石在压缩过程中的渗透系数的变化主要由孔隙闭合引起的渗透系数减小和局部损伤引起的渗透系数增大这2个因素造成的。描述孔隙闭合引起的渗透系数减小现象的理论有Kozeny-Carman模型，形式如下：

(27)

式中Kp——孔隙闭合引起的渗透系数变化；c——一个与细观结构有关的常数；n——孔隙度。

假设岩石颗粒不可压缩，Vp是孔隙的体积；Vs是岩石颗粒的体积；e0是初始孔隙比；Δe是孔隙比的变化，则体积应变为：

(28)

考虑到孔隙比和孔隙率的关系e=n/(1-n)，代入式(28)可以得到孔隙率和体积应变的关系：

(29)

式中n0——初始孔隙率。

将式(29)代入式(27)中，得到渗透系数随着体积应变变化关系为：

(30)

(31)

在描述损伤引起的渗透系数增大方面，本文将从局部损伤的角度推导出渗透系数的变化，可以看出在精度要求不高的时候，理论公式和经验公式的形式类似。在细观结构中存在着初始裂纹和后来发展的微裂纹，细观结构中的渗流主要受到这些微裂纹的影响，微裂纹中的渗流满足平行板定理，其渗透系数为：

(32)

式中u——微裂纹的张开度；ξ——微裂纹的粗糙度。

虽然式(32)比较简单，但是需要知道微裂纹的几何产状，而且本文的计算方法是基于连续介质的，所以按照微裂纹的思想，用损伤带代替微裂纹来描述渗流通道的改变，见图2。一个面积为S的截面受到单位水头梯度的作用，材料没有损伤部分的渗透系数为K0，材料的局部损伤带的宽度为λlc，在具体计算过程中可以看作是单元的特征长度。等效的思想是让损伤带和微裂纹的流量相等，具体如下：

(33)

a) 微裂纹 b) 损伤带

由式(33)可以得到损伤带的渗透系数为：

(34)

ξ代表了微观渗流通道的弯曲程度。最重要的是微裂纹开度u的计算，本文通过改进的弹性损伤模型来得到裂纹开度。在扩展有限元中，裂纹的张开是通过位移跳跃来表示的，而在渐进损伤理论中，裂纹的张开是通过损伤带的法向累积位移表示的，所以这里只考虑拉损伤造成的渗透系数变化，在不考虑卸载造成的裂纹闭合的情况下，裂纹的张开可以表示为：

(35)

假设损伤带内的应变是均匀的，式(35)的积分结果为:

(36)

要得到裂纹张开度和损伤的关系，需要求得式(14)的反函数的解，这里反函数没有显式解，需要通过迭代等数值方法得到，假设反函数的解已经得到，形式如下:

(37)

将式(37)代入到式(36)中得到微裂纹的开度：

(38)

将式(38)代入到式(34)中就可以得到损伤带的渗透系数，但因为式(37)没有显式解，所以需要简化。之前提到过理论推导的损伤带渗透系数的公式在计算精度要求不高的时候，和指数型经验公式有类似的形式，通过这一关系可以将损伤带渗透系数简化为指数型公式。将式(14)的指数项展开成泰勒级数，并取一次项，得到反函数的简化解为：

(39)

至此，由孔隙闭合引起的渗透系数减小和局部损伤引起的渗透系数增大都考虑到了。所以，能同时考虑渗透系数变化不同阶段的模型可以通过结合Kozeny-Carman公式和指数型公式得到，形式如下：

K=Kp(1+αdβ)

(40)

式中，α、β是与微裂纹几何产状有关的常数。将式(30)代入式(40)中得到损伤过程中渗透系数变化公式为：

(1+αdβ)

(41)

3 数值算例

3.1 模型边界和参数

选取水泥砂浆材料作为研究对象，具体成分是砂子、水泥和水，其水灰比为0.4，砂子重量占40%，水泥浆和砂子的密度分别为2 000、2 650 kg/m3，在不考虑气孔的情况下，细观结构中基质和砂子的体积比分别为67%、33%。砂子采用的粒径小于0.5 mm的河砂，假设砂子的粒径满足在区间[0.3,0.5]之间的均匀分布。有限元模型大小为4×4×4(mm3)，位移边界条件为在模型下边界施加固定约束；由于需要研究峰值后的应力应变规律，故采用位移边界来代替荷载边界，在模型上边界匀速施加0.04 mm的位移，对应的最大竖向应变为0.01；渗流边界为在模型上边界施加1 kPa恒定不变的孔隙水压力，下边界为自由透水边界，其他边界为不透水边界，见图3。基质和颗粒的力学和渗流参数见表1、2。

图3 水泥砂浆的细观模型

表1 细观成分的力学参数

表2 细观成分的渗流参数

3.2 结果分析

由于水泥砂浆材料的细观模型中含有基质(水泥浆)和颗粒(砂子)2种组分，分析水泥砂浆的应力应变和渗流特性必须综合基质和颗粒的应力场和渗流场，因此结果分析中将采用等效宏观应力和等效宏观渗透系数来表示水泥砂浆的特性。等效宏观应力是细观模型应力场的体积平均值：

(42)

等效宏观渗透系数的处理方式与等效宏观应力有所不同，首先根据细观模型的流速场计算等效宏观流速：

(43)

然后根据达西定律求解等效宏观渗透系数：

(44)

单轴压缩应力-应变曲线和渗透系数变化曲线见图4。可以看出宏观应力-应变曲线和室内实验数据吻合得很好，室内实验中的试件由于初始裂纹的存在，导致在弹性阶段初有一个压密阶段，峰后阶段改进的损伤模型软化曲线比较平滑。由于渗透实验数据的缺乏，所以这里不进行渗透系数曲线的对比，但可以看出渗透系数的变化曲线与Zhang[26]的试验数据一致，可以分为3个阶段，但应力峰值后的渗透系数减小阶段比Zhang的试验数据平缓，其原因可能是本文中的模型不考虑卸载后细观裂纹的闭合效应。相应于应力应变曲线上A和B点的损伤和流速云，见图5、6。在峰值应力点A，损伤在基质和颗粒接触的附近区域开始发展。在残余应力点B，损伤区域逐渐发展到充满基质。由于颗粒不透水，所以渗流只发生在基质中，颗粒中没有流速场。

a)应力-应变曲线

b)渗透系数-应变曲线

a)A点

b)B点

a)A点

b)B点

4 结论

类岩石材料的应力-损伤-渗流耦合模型需要研究变形过程中的损伤演化和渗透系数随损伤变化这两个规律。本文提出一种改进的弹性损伤模型；同时将孔隙压密阶段和裂纹扩展阶段的渗透系数变化过程统一起来，提出改进的全应力应变过程中渗透系数-损伤演化模型。最后通过在水泥砂浆材料的细观模型上实施单轴压缩数值实验，得出以下结论。

a) 数值实验的应力应变曲线和室内实验结果吻合，表明改进的弹性损伤模型用来模拟类岩石脆性材料的有效性。

b) 数值结果可以明显识别渗透系数在应力应变全过程中的不同阶段，在应力达到峰值前由于孔隙不断压密，渗透系数不断减小；在应力峰值附近，应变为0.004时，渗透系数突然增大，是因为此时开始出现损伤区域；在应力峰值附近，应变为0.007时，渗透系数达到最大值。

c) 由于缺少渗透系数的室内实验数据，本文的渗透系数-损伤演化模型的参数还需进一步通过实验确定，但渗透系数曲线趋势基本能反映出Zhang[24]的试验数据规律。