流体裂隙中的渗流计算方法：基于域模拟和压力分布

UDEC在裂隙中的流动是用域来分析的,如图5.2所示,域的标号为①～⑤,假设各种等压的流体充满了域,域与域之间通过接触发生作用,其中A～F即为接触点,其中节理为域①、③、④,两个节理的交点用域②来表示,域⑤为空洞。

在利用UDEC进行计算时,可以用网格对块体进行划分,采用三角形单元网络来对块体内的应力与位移进行计算。对于可变形块体来说,节点既可以存在于顶点上,也可以存在于块体的边缘上。如图5.2所示,点D就是处在边缘上的点,它将节理划分成了域③和域④,同时流体在裂隙中的流动就是通过这两个域来计算的。若对网格进行进一步的划分,那么在块体的边缘上将会产生更多的节点,从而节理也将被分割成更多细小的域。从这里也可以看出,利用域来分析流体在裂隙中的流动,其数值精度很大程度上取决于可变形块体的网络尺寸。在实际分析研究中,可以根据实际需求,然后确定相应的精度值,再进而确定网络尺寸。

图5.2 通过域模拟流体在裂隙中的流动

当不考虑重力时,假定流体压力在流动域中的分布是呈均匀分布的,当不考虑重力时,流体压力则按照线性分布的静水压力来计算。地下水的流动是由相邻域之间的压力差来决定的,而流动域中流体压力的大小则是由流动域汇总的中心压力差的大小来决定的。根据块体的接触条件不同,裂隙岩体中流体的流速计算方法有两种。

1)点接触

点接触(即角-边缘接触,如在图5.2中的接触点F,或角-角接触),液体由压力为p1的区域流到压力为p2的区域的流动速度表示为

式中:Kc为点接触渗透系数;Δp=p2-p1+ρwg(y2-y1);ρw为流体密度;y2、y1为域的中心的Y坐标。

2)边-边接触

在边-边接触的情况下,可以定义接触长度,如图5.2中的lD和lE分别表示触点D和E的长度,长度定义为距离的一半,到最近的左侧接触,再加上距离最近的右侧接触距离的一半。在这一个平面可以使用断裂流动的立方定律(威瑟斯庞等,1980)。然后,流速由式(5.5)得出。

式(5.5)可用于点接触,只要给这些触点分配最小长度。式(5.5)表示即使在两个域压力为零时,也可以在接触处发生流动。在这种情况下,重力可能导致流体从未完全饱和的区域迁移。然而,有两个因素需要考虑:

(1)表面渗透率应随着饱和度的降低而减小,对于零饱和度,渗透率应为0;

(2)流动性不能从零饱和度域中引出流体。

为了解决因素(1),流量(表观渗透率)可以由式(5.5)乘以一个因子fs来计算,一个关于饱和度S的函数为

式(5.9)是经验公式,它的性质是如果S=0,则fs=0;如果S=1,则fs=1(即渗透率为完全饱和不变,饱和度为0)。此外,式(5.10)的导数在S=0和S=1时为0,这在物理上是合理的。在式(5.10)中,S值被认为是流入发生时的饱和度;因此,流入不能在完全不饱和的域发生。水力隙宽公式一般为

式中:a0是法向应力为0处的裂隙隙宽;un是节理的法向开度(张开为正)。

水力隙宽最小值amin,假设为孔隙低于机械闭合不影响接触的渗透率;水力隙宽最大值amax假设为效率在显式计算中的变化。图5.3描述了UDEC中水力开度a与法向应力σn之间的关系。(www.xing528.com)

图5.3 UDEC中水力隙宽a与法向应力σ之间的关系

在UDEC中,每个时步的力学计算决定了系统的几何更新,几何更新使接触产生了新的隙宽,使区域产生了新的体积。通过上述公式可以计算出接触点上的流量。然后,考虑到流入该区域的净流量,并且由于围岩可能进一步发生运动,导致压力重分布。这个区域的压力就变为

式中:P0为前一个时步的孔隙水压力;Q为从围岩接触到的域的流量之和;Kw是流体的体积模量;ΔV=V-V0,Vm=(V+V0)/2,其中V和V0分别为新的域体积和旧的域体积。

如果按式(5.12)计算的新区域的压力为负,则压力设置为0,流出量降低了饱和度S,为

式中:S为在前面的时间域饱和度。只要S<1,则令压力值为0;应用式(5.12)这种处理方法确保了流体质量守恒,冗余的域体积或用于改变压力或用于改变饱和度。

根据上面的阐述可以看出,裂隙岩体渗流场与应力场的离散元求解过程及步骤如下:

(1)首先,把上一时步的结果作为初始条件,再利用相应的边界条件,这样就得到了本时步域与域之间的压力差;

(2)然后,根据相邻域之间的压力差计算出系统内部各个接触点的渗流量,当孔隙水压力为负值的情况时,渗透系数按照非饱和渗透系数公式进行折减,除此之外的其他情形均按照饱和渗透系数进行计算;

(3)再按照裂隙网络中的地下水的流动效应和裂隙网络自身变形来更新孔隙水压力;

(4)然后,再把求得的孔隙水压力和其他作用力(如接触力、弹性力及荷载力等)加起来,再对其进行离散元的迭代计算,最后得到新的位移和应力的分布情况。

根据新得到的块体之间的相对位移从而得到新的裂隙水力隙宽、域体积,及渗透系数等水力特征。然后重复进行上述的(1)～(4)步骤,一直到迭代收敛为止,得到裂隙岩体稳定的应力场和渗流场。此模型的一大优点就是可以使裂隙岩体中的渗流效应可以直接参与到离散元的显式平衡迭代之中来,而不必去求解大量的渗流方程组;同时,还可以考虑岩体在大变形情况下的渗透计算问题。

对于显式算法,为了保证其数值的稳定性,对流体时步有如下限制,即

式中:V为域体积;∑Ki为域周围所有接触的渗透系数的总和。

然后,再计算得出域的流体时步Δt所要求的值,取其最小值来作为最终的计算时步。

在进行流体的瞬态分析时,对于时步的要求则显得更加严格,这样一来往往会使计算耗时比较长,特别是对于裂隙隙宽远大于域的范围时的情形。除此之外,流体充满节理也会增加节理的表观节理刚度因而可能会通过减小计算时步以保证其数值的稳定性。