基于颗粒流理论的层状非均质介质等效渗透系数研究

基于颗粒流理论的层状非均质介质等效渗透系数研究
陈小茜;曾斌;王春晖;郑束宁;彭丁茂
【摘要】It is difficult to obtain the hydraulic conductivity in the work of evaluating the geological hazards that are caused by seepage flow,such as landslides and collapses,especially under the conditions of complex geological structures.This paper simulates the scenarios in which the streamline is intersected with the layered heterogeneous media at any angle using the Particle Flow Code and Darcy's Law,considering that the groundwater streamline is parallel,orthogonal or oblique to the layered heterogeneous media.The simulated hydraulic conductivity is in good agreement with the theoretical hydraulic conductivity.In conclusion,we note that the method of the particle flow simulation is more practical in evaluating the equivalent hydraulic conductivity under complicated geological conditions.The results provide a theoretical basis and reference for the seepage destruction problems in engineering construction.%在评价由于渗流引起的滑坡、塌陷等地质灾害时,由于地层结构的复杂性,导致其渗透系数难以准确获取.文章根据达西定律,利用基于颗粒流理论的数值模拟方法,对地下水流线与层状非均质介质层面成任意角度的情景进行模拟,分别选择地下水流线斜交、垂直和平行层状非均质介质层面等三种工况,定量计算含水层的等效渗透系数,模拟得到的渗透系数和基于等效渗透理论的计算值基本吻合.结果表明:颗粒流模拟方法在求取复杂地质条件下的等效渗透系数时,具有更强的实用性,可为复杂地质条件下工程建设中因渗流导致的灾害等问题的研究提供理论依据和参考.
【期刊名称】《水文地质工程地质》
【年(卷),期】2018(045)002
【总页数】6页(P7-12)
【关键词】等效渗透系数;颗粒流数值模拟;渗流;层状非均质介质
【作者】陈小茜;曾斌;王春晖;郑束宁;彭丁茂
【作者单位】中国地质大学(武汉)环境学院,湖北武汉430074;中国地质大学(武汉)环境学院,湖北武汉430074;杭州交通投资建设管理有限公司,浙江杭州 310005;浙江省交通规划设计研究院,浙江杭州 310006;浙江省交通规划设计研究院,浙江杭州310006
【正文语种】中文
【中图分类】P641.2
渗透破坏是导致多种地质灾害发生的重要原因之一，如由隧道渗水引起的地表沉降、雨水入渗导致的滑坡、土坝渗漏引起的坝体失稳等[1～3]。

为了减少渗流引起的灾害，复杂地层结构中渗流问题的研究则变得十分重要。

渗透系数是渗流问题中一个非常重要的参数，它反应了地层的透水性能[4～5]。

对于裂隙介质，前人多用等效多孔介质、双重介质等方法求取地层渗透系数[2～3]。

对于非均质孔隙介质，现阶段多采用现场试验或经验法求取渗透系数。

但现场试验费用高、耗时长，且获取的渗透系数通常是局部的、近似的；经验法求取的参数则偏差较大，难以满足计算要求的精度。

另外，等效渗透理论也常用于计算渗透系数[2]。

对于相对简单的地质
条件，基于等效渗透系数理论的计算方法简洁、精度较高；但对于复杂的地质条件，
其计算方法则很难获取精确的结果。

颗粒流程序(PFC，Particle Flow code)通过离散元法模拟规则颗粒介质之间的运
动及其相互作用，适用于固体力学大变形等问题[6～8]。

颗粒流程序不但可以表征宏观物质的物理特性，而且可以反映其他方法无法实现的细观特性。

前人已用颗粒流程序在工程地质和地质灾害等领域进行过诸多探索，包括对坝基渗流、砂土渗透、堤防渗透变形等工况的模拟研究[9～13]。

但尚无利用颗粒流程序模拟非均质层状
含水介质的等效渗透系数等方面的研究。

本文采用颗粒流离散元法探讨地层等效渗透系数的求取，从渗流过程的细观角度出发，分析颗粒的运动规律及应力分布特征，重点讨论介质结构对渗透系数的影响；并将模拟结果与经典理论公式计算结果对比，讨论数值模拟方法的可行性与有效性。

1 颗粒流数值模拟方法
PFC中用圆形颗粒代替土颗粒，引入“域”和“管道”的概念来表征和模拟真实
的流体。

“域”是一系列封闭的颗粒链，链上的每个链接都是一个连接接触。

“域”代替水压进行模拟，将压力等效成体力施加在颗粒上，在模拟中不断更新。

假想固体中流体的通道在颗粒接触处，并与颗粒接触处相切，称该通道为“管道”。

利用“管道”和“域”的连接模拟颗粒与流体的耦合作用。

PFC中流固耦合计算的流动方程、压力方程和求解条件如下[5～7]：
(1)流动方程
流体管道相当于一个平行通道，长度为L′、孔径为a，在垂直平面方向上为单位厚度，管道内的流速q(单位时间内的体积)为：
(1)
式中：K′——传导系数/(cm·s-1)；
P2-P1——相邻域的压力差/ Pa。

(2)压力方程
周围管道流入每个域的流量和为∑q，在单位时间步长Δt下，流体压力增量ΔP(流入为正)为：
(2)
式中：Kf——流体的体积模数/ kPa；
Vd——域的表观体积/mm3。

(3)耦合方式
在处理耦合过程时，采用以下方法：①通过接触的张开与闭合或接触力的变化实现通道间隙的变化；②通过改变研究区内域的力学特征来改变压力；③区域孔隙压力对其内部颗粒有推移作用。

(4)求解方法
应用显示求解方法，将流量方程应用于所有的管道，并将压力方程应用于所有的域之间交替求解。

假设某个域内存在扰动压力ΔPp，由于扰动流域里的流量可以从式(1)计算得：
(3)
由水流流入引起的响应压力变化ΔPr为：
(4)
式中：N——域所连接的管道数；
a——管道孔径/mm；
R——域周围颗粒的平均半径/mm。

(5)求解条件
保证模型运行稳定的条件是水流入引起的压力变化必须小于扰动压力，当两者相等时可求出临界时间步长：
(5)
(6)渗透系数K
(6)
式中：A——流体通过的横截面积/m2;
Δh——上下游水头差/m;
ΔL——流程长度/m。

取长度L=10 mm范围内所有截面的流量总和，在Δt时间内：
(7)
则渗透系数K的表达式为：
(8)
2 渗流模型构建
2.1 模型结构的设计
利用PFC软件生成圆形颗粒，模拟实际情况中的细土颗粒，渗流模型的尺寸为20 mm×20 mm。

经计算分析，当模型尺寸较小、其他细观参数相同时，由颗粒分布不同引起渗透系数的变化为同一数量级内的微小变化[11]，可直接忽视。

故本文中颗粒半径采用Rmax到Rmin的均匀分布，Rmax=0.58 mm，Rmin=0.35 mm。

对上、下边界的颗粒进行固定，标记为不排水边界(如图1(a)中上、下侧的绿色颗粒)。

通过PFC生成“域”和“管道”所组成的网络。

模型结构示意图见图1，图
1(a)中圆形颗粒为土颗粒，图1(b)中黑色圆点为“域”，“域”之间的浅绿色线段为“管道”，图1(c)中蓝色线段为力链。

模型左侧水压固定为2×105 Pa，右侧水压固定为0。

为减少边界对渗流模型的影响，对文中所有模型的左、右侧分别预留1.5 mm进行加压，则实际的渗流途径
为17 mm。

因无法测量某一过水断面的流量，故选取模型L=10 mm处的中间范围测量总流量。

图1 模型结构示意图Fig.1 Schematic diagram of the model structure
2.2 模拟参数的选取
用PFC模拟地层1、地层2，两地层中的细观参数——体积模量、域的管道数量
等保持一致。

为使模型符合实际渗流情况，经大量数据调试可得最佳数值，参数取值如表1所示。

表1 数值模拟基本参数Table 1 Basic parameters of the numerical simulation 体积模量Kf/kPa域的管道数量N横截面积A/m2水头差Δh/m渗径
ΔL/m10620 02200 017
PFC中的传导系数是指流体在渗流通道内的流动速度，是流体细观运动参数。

白
若虚等[11]利用PFC探讨渗透系数的细观影响因素，认为渗透系数与传导系数K′
呈线性相关。

本文通过改变传导系数模拟地层1、地层2的渗透系数。

结果见表2。

表2 不同传导系数下的流量和渗透系数Table 2 Flow and hydraulic conductivity with different conductivities地层分类地层1地层2传导系数
K′/(cm·s-1)1×10-81×10-9水流流量/mm33 16×1063 11×105渗透系数
K/(cm·s-1)1 32×10-51 32×10-6
3 模拟合理性验证
3.1 理论等效渗透系数公式
自然界中，地下水流线通过具有不同渗透系数的地层边界时，会像光线通过一种介
质进入另一种介质一样发生折射，这种现象是为了保持通过每个过水断面的流量相等。

水流从渗透性较好的地层进入渗透性较差的地层，流线会变得更加稀疏，等水线间隔变小。

对于层状非均质介质，往往假想为均质介质，该介质的水力梯度及含水层厚度和原含水层相等，Kb作为假想介质——层状非均质介质的渗透系数[13]。

在岩土体渗流问题求解中，往往应用流线和地层成特殊角度下(0°、90°)的等效渗
透系数公式较多。

对于流线和地层成任意角度的等效渗透公式，何勇等[14]提出公式(9)，并采用IGW软件进行了相关实例的分析。

通过IGW软件计算证实了此公
式在非均质含水层渗流计算中的可行性及准确性。

(9)
式中：M1M2——地层1、2的含水层厚度/m；
α——等效流线和层面的角度；
α1、α2——流线和层面1、2的角度。

图2 流线与含水层层面斜交示意图Fig.2 Schematic diagram of the scenario with the streamline oblique to the aquifer
流线与地层1和地层2斜交，地层渗透系数见表2，水流流经地层1，穿过两地层界面时发生折射，角度由α1变为α2。

地层1、2厚度相同，M1=M2(图2)。

(10)
(11)
cos2α1=0.94
(12)
cos2α2=1
(13)
由式(9)得Kb=2.39×10-6 cm/s。

3.2 模拟验证结果及其分析
由PFC软件模拟计算可得：1区的流量Q1=4.27×105 mm3/s，2区的流量
Q2=4.27×105 mm3/s，将Kb=2.39×10-6 cm/s(传导系数为1.81×10-6 m/s)
代
入PFC软件计算得出模拟流量Q3=4.55×105 mm3/s，此时Q1=Q2≈Q3，说明利用PFC软件模拟层状非均质介质的等效渗透系数的方法是合理的。

4 特殊工况等效渗透系数渗流模拟
对地下水流线斜交、垂直和平行层状非均质介质层面等三种工况进行模拟。

图3(a)中流线和地层1的夹角为73°，M1=M2=4.8 mm。

图3(b)中，流线垂直于地层
层面，M1=M2=10 mm。

图3(c)中，流线平行于地层层面，M1=M2=8.5 mm。

地层1的渗透系数K1=1.32×10-5 cm/s，地层2的渗透系数K2=1.32×10-6
cm/s。

模型左侧边界加压2×105 Pa，模型右侧边界水压为0。

红色箭头为流速矢量场，此模型边界条件为定水头边界，虽然边界水压分布对水流方向稍有影响，但流速矢量场大体上垂直等水头边界，符合地下水动力学理论。

保持地层1、2的渗透系数不变，对厚度M1、M2取不同值进行模拟。

当循环达
到一定步数时，渗流稳定，得到渗透系数模拟值Ka，其与理论渗透系数Kb具有
较高的拟合度，符合等效渗透公式。

M1、M2、Ka、Kb计算结果见表3。

图3 三种工况下的模拟结果示意图Fig.3 Schematic diagram of the simulation results in three scenarios
模型运行初期，水压从左到右分布不均，渗透系数越大，水压越大；随着模型运行时间的增加，渗流模型趋于稳定，水压分布逐渐均匀，且呈线性分布(图4)。

由此
可知，在水力梯度相同的情况下，水压在土层之间的传递时间和传导系数成反比。

三种工况下模拟值Ka和理论值Kb拟合方程的斜率均接近为1，截距小，总体上PFC计算结果与公式计算结果基本一致，Ka和Kb的关系如图5所示。

由误差分
析可知(表3)，模拟值与理论值存在一定的差异，原因在于理论公式计算中直接忽
略了地层边界对计算结果产生的影响。

在水流和地层层面斜交的工况下，保持其他参数不变，改变模型大小，发现非均匀介质地层等效渗透系数计算结果的可靠性与模型尺度有关：对于小尺度模型，可直接用PFC软件求取等效渗透系数；但当模型尺度较大时，地层的非均质程度更高，建模难度较大，建议选取不同位置的地层分别建模计算[15]。

也有学者建议保持地层结构不变，对模型等比例缩放，但研究结果仍需通过实践检验。

因此在大尺度范围内，如何应用PFC软件精确求取等效渗透系数，尚需深入研究。

表3 不同地层厚度下的模拟渗透系数和理论渗透系数的结果对比Table 3 Comparison between the simulated and theoretical hydraulic conductivities with different rock thicknesses地层1厚度M1/mm地层2厚
度M2/mm渗透系数K1/(cm·s-1)渗透系数K2/(cm·s-1)理论等效渗透系数
Kb/(cm·s-1)模拟渗透系数Ka/(cm·s-1)误差8 02 01 32×10-51 32×10-64
63×10-64 35×10-6-6 04%7 03 01 32×10-51 32×10-63 53×10-63 36×10-6-4 81%6 04 01 32×10-51 32×10-62 85×10-62 75×10-6-3 51%流线斜交于地
层5 05 01 32×10-51 32×10-62 39×10-62 26×10-6-5 44%4 06 01 32×10-
51 32×10-62 05×10-61 99×10-6-3 00%3 07 01 32×10-51 32×10-61
80×10-61 72×10-6-4 44%2 08 01 32×10-51 32×10-61 61×10-61 55×10-6-3 73%2 517 51 32×10-51 32×10-62 81×10-62 54×10-6-9 61%5 015 01
32×10-51 32×10-64 29×10-64 24×10-6-1 17%7 512 51 32×10-51 32×10-65 78×10-65 68×10-6-1 73%流线平行于地层10 010 01 32×10-51 32×10-67 26×10-67 18×10-6-1 10%12 57 51 32×10-51 32×10-68 75×10-68 83×10-
60 90%15 05 01 32×10-51 32×10-610 20×10-61 05×10-62 15%17 52 51 32×10-51 32×10-611 70×10-61 20×10-52 56%1 016 01 32×10-51 32×10-61 39×10-61 42×10-6-2 33%3 513 51 32×10-51 32×10-61 62×10-61
64×10-61 23%6 011 01 32×10-51 32×10-61 93×10-61 92×10-60 51%流线垂直于地层8 58 51 32×10-51 32×10-62 40×10-62 45×10-62 08%11 06 01 32×10-51 32×10-63 16×10-63 20×10-61 26%13 53 51 32×10-51 32×10-64 63×10-64 49×10-6-3 02%16 01 01 32×10-51 32×10-68 63×10-67
99×10-6-7 41%
图4 渗流与地层平行情况下不同时段水压分布图Fig.4 Schematic diagram of pressure distribution at different time intervals in the scenario with the streamline parallel to the layer
图5 模拟结果Ka与理论渗透系数Kb关系图Fig.5 Relationship between the simulation results Ka and theoretical coefficient of permeability Kb
5 结论
(1)采用PFC计算非均质层状介质的等效渗透系数得到的结果与使用理论公式计算得到的结果基本一致，Kb与Ka关系式的斜率约为1，最大误差为-9.61%；表明基于颗粒流理论的数值模拟方法适用于小尺度条件下等效渗流模型的计算及等效渗透系数的求取。

(2)采用PFC数值模拟方法，可得到渗流模型内各质点、任意时段的水压、流速矢量场，且水力梯度相同时，水压在土层之间的传递时间和传导系数成反比。

这一结果可为研究管涌、岩溶塌陷等的临界水力梯度、临界水压提供参考依据。

(3)对复杂地质条件下的非均质含水层的渗透系数的求取，应采用适用性更强的PFC进行模拟计算。

参考文献：
[1] 刘印,张冬梅,黄宏伟.盾构隧道局部长期渗水对隧道变形及地表沉降的影响分析[J].岩土力学,2013,34(1): 290-298.[LIU Y,ZHANG D M,HUANG H W. Influence of long-term partial drainage of shield tunnel on tunnel deformation and surface settlement[J].Rock and Soil Mechanics,2013,34(1): 290-298.(in Chinese)]
[2] 张奇华,李玉婕,袁东,等. 地下水封洞库水幕孔注水试验及岩体等效渗透参数分析[J]. 岩土力学, 2015, 36(9):2648-2658. [ZHANG Q H,LI Y J,YUAN D, et
al.Water injection test about water curtain borehole for underground water-sealed cavern and analysis of rock equivalent permeability parameter[J].Rock and Soil Mechanics, 2015, 36(9):2648-2658. (in Chinese)] [3] 简文星,丁志雷,李世金.巴东县城白土坡裂隙岩体等效渗透系数数值分析[J].水文地质工程地质,2015,42(3):65-68.[JIAN W X, DING Z L, LI S J. Numerical analysis of equivalent coefficient of permeability of the Baitupo slope in Badong County [J].Hydrogeology & Engineering Geology,2015,42(3):65-68. (in Chinese)]
[4] 张人权,梁杏,靳孟贵,等.水文地质学基础[M].北京:地质出版社. 2010. [ZHANG R Q,LIANG X,JIN M G, et al. Fundaments of hydrogeology
[M].Beijing:Geological Publishing House,2010.(in Chinese)]
[5] JACOB Bear. Hydraulics of groundwater[M]. New York: Dover Publications, 2007.
[6] Itasca Consulting Group Inc.PFC2D particle flow code in 2 dimensions-optional features [M].Minneapolies: Itasca Consulting Group Inc,2005. [7] Itasca Consulting Group. Verification Problems and Example Applications[M]. Minneapolies: Itasca Consulting Group Inc,2005.
[8] Cundall P A. PFC2D User ’s Manual (Version2. 0) [R]. Minnesota: Itasca Consulting Group Inc,1999.
[9] 李识博,王常明,王刚城,等.松散堆积物坝基渗透淤堵试验及颗粒流模拟[J].水利学报,2012,43(10):1163-1170.[LI S B, WANG C M, WANG G C, et al. Experimental study on seepage and silting and seepage of granular pile foundation [J]. Journal of Hydraulic Engineering, 2012, 43(10): 1163-1170. (in Chinese)]
[10] 张刚,周健,姚志雄.堤坝管涌的室内试验与颗粒流细观模拟研究[J].水文地质工程地质,2007,34(6): 83-86. [ZHANG G,ZHOU J, YAO Z X. Study on mesomechanical simulation of piping with model tests and PFC2D[J]. Hydrogeology & Engineering Geology, 2007,34(6): 83-86. (in Chinese)] [11] 白若虚,程雪松,郑刚.关于土渗透系数颗粒流细观参数研究 [J].低温建筑技术,2012(1):64-67. [BAI R X, CHENG X S, ZHENG G. Study on the mesoscopic parameters of soil Hydraulic conductivity [J]. Low Temperature Architecture Technology, 2012(1): 64-67. (in Chinese)]
[12] 倪小东,朱春明,王媛. 基于三维离散-连续耦合方法的堤防工程渗透变形数值模拟方法[J].土木工程学报, 2015, 48(增刊1): 159-165. [NI X D, ZHU C M, WANG Y. Analyzing the behavior of seepage deformation by three-dimensional coupled continuum-discontinuum methods [J]. China civil Engineering Journal,2015,48(Sup1): 159-165. (in Chinese)]
[13] 薛禹群.地下水动力学[M].北京:地质出版社,1997.[XUE Y Q. Groundwater dynamics[M]. Beijing: Geological Publishing House, 1997. (in Chinese)] [14] 何勇,袁明,陈旭,等.等效渗透系数在非均质含水层渗流计算中的应用[J]. 电网与清洁能源, 2010, 26(4): 73-76. [HE Y, YUAN M, CHEN X, et al. Application of
equivalent penetration coefficient method in seepage flow calculation of non-homogeneous aquifer [J]. Power System and Clean Energy, 2010,
26(4): 73-76. (in Chinese)]
[15] 张弛,仵彦卿,覃荣高.渗透系数升尺度对非均质含水层溶质迁移影响研究[J].水文地质工程地质,2014,41(5): 20-25.[ZHANG C,WU Y Q, TAN Y G. Research on effect of hydraulic conductivity upscaling on groundwater solute transport in heterogeneous aquifer [J]. Hydrogeology & Engineering Geology, 2014,41(5):20-25. (in Chinese)]。