王梦辉，涂福彬*，童俊，焦玉勇，郑翔
(1.中国地质大学(武汉)工程学院，武汉 430074；2.中建一局集团建设发展有限公司，北京 100102)
渗透性是裂隙岩体的一个重要特性，是分析地下油气运移与污染物迁移等渗透问题的基础。在岩体中由于裂隙的渗透性远远大于岩块的渗透性，因此裂隙是岩体中流体运移的主要通道，裂隙渗透性的研究足以确定裂隙岩体的渗透性。在实际岩体中，裂隙往往相互交错形成网状结构，裂隙岩体的渗透性研究存在一定的困难，对岩体中单一裂隙渗透规律的研究是了解复杂裂隙岩体渗透性的基础。
在岩体渗透性预测方面，Snow[1]通过光滑平行板裂隙水流实验提出了著名的立方定律，即认为通过裂隙的渗流量与裂隙宽度的三次方成正比。而自然界实际岩体的裂隙面往往是弯曲粗糙的，光滑平行板模型在表征裂隙表面时存在误差[2-3]。任志善等[4]通过对数值分析模型研究发现裂隙倾角、间距、隙宽和迹长对岩体的渗透特性均有明显影响。当流体流经粗糙裂隙时，由于节理表面粗糙度、接触面积的影响，流体渗流路径是曲折的[5]。朱红光等[6]研究表明，立方定律只适用于壁面粗糙度足够小的情况。为此，许多学者对立方定律进行了修正。Rose等[7]提出了迂曲度的概念，随后基于迂曲度修正的立方定律计算公式被用于预测含天然裂隙岩体的渗透性。夏才初等[8]采用节理表面分形维数和空腔组合形貌等参数来计算通过粗糙裂隙的渗流流量。肖维民等[9-12]通过归纳实验结果，提出了考虑曲折效应的粗糙节理渗流流量计算新公式。速宝玉等[13]对含天然裂隙的岩体渗流进行了研究。覃源等[14]采用实验手段研究了不同节理粗糙度对单裂隙渗流的影响。段慕白等[15]通过实验研究了裂隙粗糙度及隙宽对渗流的影响。李崴等[16]通过数值模拟方法研究了辐射流情况下裂隙岩体的渗透性。王玮等[17]采用基于正、负相关指标的复合因子模型评估裂隙岩体的渗透系数。
采用数值方法求解不可压缩黏性流体稳态流动的Navier-Stokes方程是近年来精细确定裂隙岩体渗透性的主要手段。Nguyen等[18]采用快速傅里叶变换方法求解细观表征体元内不可压缩黏性流体稳态流动的Navier-Stokes方程，并通过平均流速和压强梯度之比确定多孔介质的渗透系数。王志良等[19]、张戈等[20]采用格子Boltzmann方法模拟岩体裂隙内的流体流动。上述方法只模拟了一定长度裂隙内的流体流动，忽略了裂隙延伸长度对渗流的影响。钟振等[21]采用数值模拟手段发现粗糙单裂隙的渗透性存在明显的尺寸效应。Zhang等[22]通过微裂隙渗流实验也发现裂隙的渗透性存在尺寸效应。
为此，有必要施加周期性边界条件(periodic boundary conditions,PBC)来消除裂隙延伸长度所引起的误差。首先采用不可压缩单元离散不可压缩流体的稳态纳维-斯托克斯方程，并施加周期性边界条件构造裂隙岩体渗透性的有限元分析方法。接着采用该方法分析了不同节理面粗糙度系数(joint roughness coefficient,JRC)裂隙的渗透性。
通过求解Navier-Stokes方程，可以获得流体在裂隙内的流动速度[18]。对于渗流问题，流体流动往往处于恒定层流状态。假定流体不可压缩，忽略体积力，取长度为L的表征体元(representative volume element,RVE)，如图1所示。
图1 无限长平行板周期性结构示意图Fig.1 Periodic structure diagram of infinite length parallel plate
不可压缩黏性流体稳态流动的Navier-Stokes方程为
(1)
∇·u=0
(2)
对表征体元施加周期性边界条件，则有
u(x+L,y)=u(x,y)
(3)
(4)
式中：x为平行于无限长平行板的方向；y为垂直于无限长平行板的方向。
采用四节点等参数单元对速度和周期性分布的压强进行插值，并得离散方程为
(5)
式(5)中：K为刚度矩阵；G为耦合矩阵；F为外力等效节点力向量。
(6)
(7)
(8)
式(5)的数值稳定形式为
(9)
式(9)中：π为保证数值求解稳定而引入的伪变量;L为周期性压强对应的刚度矩阵;Q为周期性压强和稳定变量的耦合刚度矩阵;M为稳定变量对应的刚度矩阵。
式(9)可以简化为
(10)
(11)
(12)
(13)
基于式(10)～式(13)，编写了ABAQUS子单元程序实现裂隙渗流的有限元求解，并通过Python脚本施加周期性边界条件。
采用不同JRC值来描述裂隙表面粗糙度，探究裂隙面粗糙度系数、裂隙宽度以及裂隙两表面相互错动对渗流的影响。国际岩石力学学会推荐给出的10级JRC值对应的粗糙度剖面曲线如图2[23]所示，JRC值为0表示裂隙面为光滑平面。
0～10 cm为裂隙的长度图2 不同JRC值范围的粗糙度剖面[23]Fig.2 Roughness profiles with different JRC values[23]
地质构造的作用极易引起裂隙两表面的相互错动，而错动将影响裂隙的宽度，进而改变裂隙内部流体的流速，影响渗透性。JCR值为18～20时，两表面相互错动的裂隙模型如图3所示。
图3 两表面相互错动的裂隙模型Fig.3 Fracture model with dislocation of two surfaces
Rose等[7]提出了迂曲度的概念来表征渗流路径的曲折特性，迂曲度τ为沿着裂隙走向由一端走至末尾的实际长度Le与裂隙两端直线距离L之比，可表示为
(14)
将立方定律公式采用迂曲度修正之后得
(15)
现从无限长的光滑平行板中选取表征体元，如图1所示。通过施加周期性边界条件来计算无限长的光滑裂隙间的渗流。并将有限元计算的结果与泊肃叶流动的解析解进行对比。泊肃叶流动的流速分布为
(16)
选取参数如表1所示。将有限元方法计算得到的结果与泊肃叶流动的解析解进行对比，如图4所示，有限元计算结果与解析解相吻合。可见，所提出的裂隙渗流有限元计算方法是可行的。
表1 模型参数Table 1 Parameters of the model
图4 有限元方法与解析解对比Fig.4 Comparison between result of finite element method and analytical solution
对JRC值10～12隙宽1 mm的一条裂隙在水力坡度为0.01的情况下进行有限元分析，比较有无周期性边界条件计算结果的差异，如图5所示。不施加周期性边界条件裂隙左右两侧细观上的流速不同，这与沿走向延伸较长的裂隙中渗流不符合，施加周期性边界条件对延伸较长的裂隙非常重要。
图5 有无周期性边界条件结果对比Fig.5 Comparison of results with and without periodic boundary conditions
通过模拟隙宽3 mm的不同JRC值裂隙内流体流动，探究相同水力坡度下节理面粗糙程度对通过裂隙的流体平均流速和最大流速的影响。不同JRC值裂隙内渗流的最大流速和平均流速如图6所示。随着裂隙表面粗糙程度的增大，渗流的平均流速整体上呈逐渐减小的趋势。细观上渗流的最大流速并不是随着节理面粗糙程度的增大而逐渐增大，而是与节理的局部起伏角相关，这与覃源等[14]的研究结果相似。
图6 不同JRC值裂隙内流体流动情况(隙宽3 mm)Fig.6 Fluid flow in fracture with different JRC values (crack width 3 mm)
对不同水力坡度下隙宽3 mm的不同JRC值的裂隙进行有限元分析，计算不同水力坡度下通过不同JRC值裂隙的单宽流量值，结果如图7所示。通过不同JRC值裂隙的单宽渗流量与水力坡度成正比，符合达西定律。
计算隙宽3 mm的不同JRC值裂隙在水力坡度为0.01时的渗流量，将有限元计算结果与引入迂曲度修正后的立方定律公式计算所得渗流量进行比较，如图8所示。在节理面粗糙程度较低时，单宽流量的迂曲度修正值与有限元计算结果较为接近，但是当节理面粗糙程度较大时采用迂曲度修正的立方定律公式高估了单宽流量的值，在JRC值为18～20时约高估了19.4%。通过裂隙的单宽流量整体上随着节理面粗糙程度的增加而降低。
图8 数值计算结果与迂曲度修正值对比(隙宽3 mm)Fig.8 Comparison between the numerical results and the modified value of tortuosity (crack width 3 mm)
取水力坡度为0.01，对不同JRC值的裂隙取不同的隙宽分别进行有限元计算，得到通过裂隙的流体平均流速如图9所示。裂隙宽度为1 mm时，和光滑平行板裂隙内流体的平均流速相比，JRC值为18～20的裂隙内流体的平均流速下降了27.66%，当裂隙宽度为3 mm和5 mm时，这一值分别下降了21.38%和19.71%。由此可见，随着裂隙隙宽的增大，裂隙表面的粗糙程度对渗流的影响减小。当裂隙宽度由1 mm分别变化到3 mm和5 mm时，光滑平行板裂隙内的流体平均流速分别增加了9倍和25倍，而JRC值为18～20的裂隙内流体的平均流速分别增加了9.64倍和27.74倍。由此可见，JRC值越大，裂隙隙宽变化对渗流的影响越显著。这是由于粗糙度即裂隙表面的起伏作用范围有限，随着裂隙隙宽的逐渐增大，粗糙度对渗流的影响逐渐减小。
图9 不同宽度裂隙内流体平均流速对比Fig.9 Comparison of average velocity of fluid in different width fractures
不同JRC值及隙宽下裂隙内流体最大流速对比如图10所示。相同隙宽的裂隙，流体最大流速并没有随着裂隙面JRC值的变化呈现出规律性的变化趋势。流体最大流速与裂隙面的局部粗糙度以及粗糙作用范围有关，在裂隙宽度相同的情况下，裂隙面局部曲折程度越大流体最大流速越大。最大流速反映了裂隙内流体流动速度的不均匀性，裂隙宽度越小节理面粗糙程度对最大流速的影响越大。当裂隙宽度为5 mm时，各JRC值的裂隙中流体最大流速差别已较小，这说明节理面JRC值对裂隙内渗流的影响与节理面最大起伏和裂隙宽度的比值有关，当这一比值较小时，节理面JRC值对裂隙内渗流的影响也较小。
图10 不同宽度裂隙内流体最大流速对比Fig.10 Comparison of relative maximum velocity of fluid in different width fractures
选取JRC值为18～20、初始隙宽5 mm的一条裂隙，水力坡度值取为0.01，计算裂隙两表面相对错动不同距离时流体的渗流速度。裂隙两表面相互错动5 mm时的细观流速分布如图11所示。在隙宽最小处，流体流速最大，这是由于通过该裂隙的流体是稳定流，各时刻通过各横截面的流量相等，在隙宽最窄处流体的渗流断面面积最小，因而流速最大。在裂隙中靠近裂隙面的位置流体流速最小，裂隙中央流体流速最大。
图11 裂隙两表面相对错动5 mm时的流速分布Fig.11 Distribution of fluid velocity with relative displacement of 5 mm between two surfaces of fracture
裂隙两表面不同错动距离下渗流速度如图12所示。裂隙内渗流的平均流速随着裂隙两表面的相对错动距离的增加而逐渐降低，错动距离增加，裂隙最小宽度变窄，渗流量降低。通过裂隙的最大渗流速度并未随着裂隙两表面相对错动距离的增大呈现规律性的变化。
图12 裂隙两表面相对错动状态下的渗流速度(JRC为18～20)Fig.12 Seepage velocity at relative dislocation of fracture surfaces (JRC is 18～20)
选取JRC值为18～20、初始隙宽5 mm的一条裂隙，计算裂隙两表面相互错动不同距离时不同水力坡度下的单宽流量值，计算结果如图13所示。发现在裂隙两表面相互错动后，通过裂隙的单宽渗流量仍然与水力坡度呈线性关系，即符合达西定律。裂隙两表面相互错动距离越大，裂隙的渗透系数越小。裂隙两表面相互错动后，裂隙形状变得复杂，流体的渗流路径也变得更加曲折，因此裂隙的渗透系数减小。
图13 裂隙两表面不同错位情况下单宽流量随水力坡度的变化Fig.13 Variation of unit width seepage flow with hydraulic gradient under different dislocation of fracture surfaces
使用ABAQUS自定义不可压缩单元来求解不可压缩黏性流体稳态流动的Navier-Stokes方程，引入周期性边界条件求得压强梯度作用下表征体元的渗流速度，数值预测JRC值对渗流的影响，得出以下结论。
(1)仅选用一小段裂隙作为表征体元而忽略裂隙延伸长度计算得到的裂隙渗透性有一定的误差，施加周期性边界条件能有效地消除这种误差。
(2)节理面粗糙度系数JRC值对裂隙的渗透性有影响，随着裂隙宽度的增大，影响减小。当裂隙宽度和水力坡度一定时，裂隙的渗透性整体趋势上随着JRC值的增大而减小。