基于多孔介质土体分形特征的渗透系数研究

王 宇，谷艳昌 ，王士军 ，段祥宝

(1.江苏开放大学 建筑工程学院, 江苏 南京 210036； 2.南京水利科学研究院, 江苏 南京 210029； 3.水利部大坝安全管理中心, 江苏 南京 210029)

多孔介质土体由散体颗粒、胶结物质和不同阶次、不规则分布的孔隙及薄弱夹层构成，其颗粒与孔隙从原子尺度到晶粒尺寸范围均表现分形特征[1-4]。分形维数可以描述并量化具有分形特征的多孔介质土体的自相似性、复杂性及不规则性[5-8]。具体而言，质量分形维数综合反映了细颗粒含量及内在均匀程度，描述了土体颗粒分布特征；孔隙分形维数集中体现了孔隙通道大小及发育程度，反映了土体孔隙分布特征。多孔介质土体渗透破坏过程中，可动细颗粒从孔隙通道中不断流失，土体渗透系数、孔隙率、不均匀系数等物理力学参数发生变化，分形维数随之改变。研究多孔介质土体渗流分形特征，控制并预测渗透破坏形成及发展，对堤坝工程长效安全运行具有重要意义。

20世纪80年代，分形理论广泛应用于岩土工程，揭示了多孔介质土体分形特征，即土骨架和孔隙是分形的。流体在土骨架和孔隙形成的分形空间渗流时，也应当具有分形特征及其他一些特殊性质，各种分形模型在模拟多孔介质渗流问题上取得了较大成功[9-13]。近年来，国内外学者重点关注土体水力参数求取分形方法和渗流过程中分形维数变化规律，得出了一些有益的结论。徐永福等[14]建立了非饱和土孔隙分布的分形模型，推导了非饱和土的土-水特征曲线、渗透系数及扩散系数的表达式，提供了一种可供工程实际采用的非饱和土力学参数的简单方法。Tao等[15]利用土-水特征曲线分形模型给出了非饱和土的土-水特征曲线预测方法，揭示了不同初始孔隙比条件下分形维数近似不变的规律。段祥宝等[16]开展了天然土体及室内配制土样渗透变形试验，初步探索了土体孔隙率、不均匀系数、中值粒径、干密度与分形维数之间关系，指出土体质量分形维数可以鉴别渗透破坏形式。陶高梁等[17]开展了不同颗粒级配的砂土渗透破坏试验，分析了临界水力梯度与分形维数之间的关系，认为分形维数越大，临界水力梯度越小，解释了分形维数与渗透破坏的影响机制。

目前利用分形理论研究渗流集中在非饱和土持水能力（土水特征曲线），较少涉及多孔介质土体渗透能力，渗透系数及其影响因素研究尚显不足，渗透系数与分形维数之间非线性关系需要深入研究。本文基于多孔介质土体分形特征，推导渗透系数和孔隙率与分形维数之间分形关系解析式，探明多孔介质土体渗透系数影响因素，采用渗透破坏试验数据进一步论证渗透系数和孔隙率与分形维数之间关系，从而揭示多孔介质土体渗透破坏非线性特征，为土体渗透破坏形成机制及发展过程提供理论依据。

1 多孔介质土体分形特征

岩土体的变形、断裂、孔隙、密度等物理力学性质均应具有分形特征，这一思想已被大量试验结果所证实[18-21]。对于多孔介质土体而言，粒度分布是地质勘察的基本试验项，也是岩土力学性能评估的重要指标之一。Turcotte[22]提出粒子数目分布与尺寸之间具有分形关系，之后Tyler等[23]提出颗粒体积分布与尺寸也存在分形关系。鉴于土体颗粒数目和体积在实际工程测量中较为困难，在假设不同粒径土体颗粒具有相同密度前提下，可推导土体质量分布函数的分形表达式[7]。工程地质勘察中，土体颗粒级配曲线是以小于某一粒径d的土体质量百分含量M的分布特征来表示的，根据土体颗粒质量分布分形规律[16]，有：

式（1）两边取对数，有：

式中：直线斜率K=3-D，则质量分形维数D=3-K。式（2）表明土体颗粒级配曲线在双对数坐标轴系统中，小于某一粒径d的土体质量百分含量M具有线性分布规律。实际工程中可利用地勘资料，绘制土体颗粒级配曲线，在双对数坐标轴下绘出质量百分含量与粒径之间线性相关曲线，即土体粒度分布曲线，通过曲线斜率推算土体质量分形维数。土体质量分形维数定量表征土体颗粒质量累积特性，可以反映整个土层的物理力学性质。以某砂砾石为例，粒度分布曲线如图1所示，R2=0.99表明线性关系较好，直线斜率K=0.315 6，土体质量分形维数D=3-0.315 6=2.684 4。

图1 某砂砾石粒度分布曲线Fig.1 Sandy gravel size distribution curve

值得注意的是，由于多孔介质土体具有统计自相似特点，土体颗粒仅在一定的粒径范围内具备分形特征，即当x∈[x1,x2]，x1和x2为自相似性范围的上下限，式（2）分形规律成立。

2 不同试验土样统计分析

2.1 试验材料与方法

质量分形维数可以量化具有统计自相似特征的多孔介质土体颗粒分形结构，但由于拓扑结构存在空间差异性，统计自相似的多孔介质土体分形结构适用于一定的粒径范围。具有自相似性的粒径范围称为无标度区[12]，无标度区的确定是准确计算土体质量分形维数的前提。基于文献[24]中的37组渗透破坏试验成果，选取具有代表性的10种试验土样（见表1），进一步研究不同类型土体无标度区和质量分形维数。

表1 不同试验土样颗粒级配Tab.1 Particle size distribution of different experimental soils

试验土样来源于钱塘江中游金华段两岸，包括31组天然土样，6组人工配制土样。天然土样采用环刀法，将环刀刃口向下垂直压入原状土中，直至环刀筒充满土体，用削土刀切开环刀周围土壤，取出已装满土的环刀，削去并擦净环刀外面的土，尽量减少土体扰动；人工土样采用不同孔径振动筛配制。

2.2 无标度区的确定

分形几何基本理论中分形维数有容积维数、相似维数、关联维数、信息维数等类型，土体质量分形维数本质上属于容积维数。因此，确定多孔介质土体分形结构无标度区与计算容积维数的方法一致。计算容积维数方法有人工判定法、拟合误差法、强化系数法、相关系数检验法等。人工判定法优点在于计算简便、应用广泛，判定无标度区间精度较高，可避免无标度区出现较大“漂移”现象。

容积维数定义为N(r)=cr-D，其中r为分形结构量度的单位半径；N(r)为r量度分形结构所得的数量；c为分形结构不随r而改变的常数[1]。对容积维数定义式两边取对数，有 lnN(r)=lnc-Dlnr。因此，通过lnN(r)-lnr关系曲线可判别土体结构无标度区间。根据土体 lnN(r)-lnr关系曲线直线段范围可以确定多孔介质土体无标度区，不同试验土样无标度区范围如图2所示。图2中散点变化趋势出现2个拐点，拐点之间中部散点拟合后近似为一条直线，这一直线段所对应的量测尺度范围为无标度区间，即土体颗粒自相似性的尺度区域，其线性相关越好则土体分形结构的无标度区越准确。试验土样统计分析结果表明，不同类别试验土体无标度区差别很大，土体粒径越大，具有统计自相似性的无标度区范围越大，此时无标度区内细颗粒含量决定了土体大部分性质。圆砾无标度区粒径范围0.250～20.000 mm，中粗砂无标度区粒径范围0.250～10.000 mm，粉细砂无标度区粒径范围0.075～0.500 mm，黏土无标度区粒径范围0.002～0.005 mm。

图2 不同试验土样无标度区范围Fig.2 Scale-invariant space of different experimental soils

2.3 质量分形维数的计算

土体质量分形维数表征颗粒大小、数量及分布的均匀程度。以试验土样圆砾为例（土样编号为2，见表2），计算土体质量分形维数，粒度分布曲线如图3所示。由于土体颗粒具有非均质性、非连续性和不规则性，利用单一分形维数描述土体分形结构不够全面。若不分段统计，相关系数为0.67，质量分形维数D=3-0.865 2=2.134 8，可靠性较差；若分两段统计，各段拟合度较好，相关系数大于0.90，表明土体有两个质量分形维数，若对粒径小于0.10 mm进行统计，则质量分形维数D1=3-3.186 6＜0，表明土体均匀程度很高，此时已不具备分形规律；若对粒径大于0.10 mm的土体进行统计，则质量分形维数D2=3-0.265 4=2.734 6，相关系数为0.94，可靠性较高。该圆砾大于0.10 mm的颗粒含量约96%，对土体结构稳定性起主控作用，质量分形维数D=D2=2.734 6是合理的，此时质量分形维数可以表征土体大部分物理力学性质。此外，单一分形维数拟合精度低，无标度区范围偏大，质量分形维数偏小，平均偏离幅度12.12%；多段分形维数拟合曲线的相关系数在0.88～0.99，利用多段分形维数确定土体无标度区和计算质量分形维数更为精确，更符合工程实际。

图3 某圆砾的粒度分布曲线Fig.3 Round gravel size distribution curve

表2 不同试验土样质量分形维数与无标度区统计结果Tab.2 Statistical results of mass fractal dimension and scale-invariant space of different experimental soils

无标度区内土体细颗粒含量是决定土体质量分形维数的主要因素。当质量分形维数D=0时，表明土体由粒径相等的颗粒组成；按粒径分布确定的土体质量分形维数与其拓扑维数无直接关系，土体颗粒粒径愈小，无标度区细粒含量愈高，则质量分形维数就愈大，表明土体颗粒组成愈不均匀。需要说明的是，理论上确定型分形模型质量分形维数D∈（2，3），在实际应用中，对于相对均匀的砂性土、黏性土等细粒料，具有统计意义的分形模型质量分形维数D∈（1，2）也是有意义的。简言之，多孔介质土体质量分形维数表征颗粒质量累积特性，体现了土体颗粒级配、密实度、渗透性等物理力学性质，这对研究多孔介质土体渗透系数具有重要意义。

3 渗透系数与孔隙率分形公式理论推导

质量分形维数可以描述土体颗粒分布及质量累计的分形特征，但不能完全反映土体的渗透性能。多孔介质土体发生渗透破坏时（如管涌），土体可动细颗粒从孔隙通道不断流失，质量分形维数发生变化，在渗流侵蚀作用下，孔隙通道逐渐扩展，渗透性能也随之改变。因此，探究多孔介质土体渗透性能与质量分形维数之间分形关系尤为必要。渗透系数是反映土体渗透性能的重要指标，由于土体结构不规则性和复杂性，目前尚不能用统一的数学解析公式来表达渗透系数，工程上通常采用一些经验数值或经验公式，或根据室内土工试验、现场原位试验和数值试验反演来确定土体渗透系数。考虑到土体颗粒与孔隙从原子尺度到晶粒尺寸范围内均表现出分形特征，从测度与尺度幂律关系出发，分析多孔介质土体微观孔隙结构，推导渗透系数和孔隙率与分形维数之间解析关系，进一步探究多孔介质土体渗流分形特性。

假设多孔介质土体由一系列不等径平行毛管束构成，断面及空间分布见图4。单位截面上半径大于R的 毛细管总数N(R)与R呈幂律关系[12]，即：

图4 多孔介质毛管束模型Fig.4 Pipe bundle model of porous medium

式中：a为 分形系数；R为 孔径；D为分形维数。

式（3）两边求导，孔径R～R+dR之间毛管数为：

根据圆管中层流运动的普瓦索伊（Poiseuille）公式[25]，单个孔隙通道流量q(R)为：

式中：γ为水的重度；µ为黏滞系数；h1-h2为水头损失；L为 渗径长度；A为截面面积。

通过某界面的总流量Q可以看成该界面内每个孔隙通道的流量q(R)在最小孔径Rmin和最大孔径Rmax之间进行积分，总流量Q为：

将式（4）、式（5）代入式（6），整理得：

根据达西定律，有：

式中：k为渗透系数；v为渗流速度；J为水力梯度。

将式（7） 代入式（8），整理得：

根据孔隙率定义，有：

式中：n为 孔隙率；Vp为毛管束孔隙体积；V为总体积。

将式（4）代入式（10），整理得：

由式（9）和式（11）可知，多孔介质土体渗透系数和孔隙率是分形系数、分形维数及孔隙尺寸的函数。若多孔介质土体分形维数、分形系数和孔隙尺寸不同，则渗透系数和孔隙率也会不同。严格意义上，渗透系数与孔隙率之间很难找到一种解析关系，鉴于多孔介质土体多孔隙性和孔隙分布不规则性，其最大孔径远大于最小孔径，即：Rmin/Rmax→0，式（9）和式（11）可简化为：

进一步地，多孔介质土体渗透系数和孔隙率与分形维数之间解析关系表达式为：

式中：C、b为常数

由式（14）可知，渗透系数由分形系数、孔径大小、分形维数、流体黏滞系数共同决定；多孔介质土体渗透系数和孔隙率与分形维数之间解析关系为幂律函数，属于典型的非线性关系。

4 渗透系数与孔隙率分形关系试验验证

分析天然土样和人工土样渗透破坏试验数据[24]，进一步论证渗透系数、孔隙率与分形维数之间非线性关系。试验相关过程如下：①具体操作参照《土工试验方法标准》（GB/T 50123—2019），试验采取连续加压、逐级饱和的加水方式；②试验一律在室温下进行，土体几何尺寸严格按照土工试验方法标准执行，试验土体全部进行烘干处理，在充分饱和情况下开展试验；③试验过程中测量沿程测压管渗压水位及溢流量，观察溢水浑浊程度及土料表面变化情况，当出现水力坡降急剧下降或渗透系数和溢流量突增10倍以上时，试验结束。

试验结果表明：土体分形维数和孔隙率理论值与试验值变化趋势较为吻合，孔隙率随着分形维数的增大而减小（图5）。孔隙率试验值在理论值拟合光滑曲线附近呈现区域性振荡，最大相对误差为18%，可能与公式推导假定和试验系统误差有关。土体分形维数和渗透系数理论值与试验值变化趋势基本一致（图6）。渗透系数试验值比理论值偏大，最大相对误差为15%，可能与试验装置和试验条件有关。试验结果论证了多孔介质土体孔隙率和渗透系数与分形维数成反比关系。当分形维数D＞2.83时，孔隙率随着分形维数增大而减小的规律更为显著，而渗透系数随着分形维数增大而减小的规律不明显。究其原因，分形维数表征单位体积下颗粒逐级积累特性，填充在土骨架与孔隙之间可动细颗粒数量越多，土体孔隙率越小。假设一种理想情况，细颗粒不断填充土体骨架孔隙，土体结构形态随之发生变化，由松散的砂砾石逐渐过渡到中等密实的粉细砂，最后形成密实的黏土。当过渡到黏土阶段时，在颗粒之间吸着水和薄膜水形成的黏聚力作用下，孔隙率进一步减小，由于水气和气泡的存在，真正流过孔隙的水量减少，鉴于渗透系数表征流体通过骨架孔隙的难度程度，因而渗透系数减少相对缓慢。

图5 土体分形维数和孔隙率理论值与试验值对比Fig.5 Comparison between theoretical values and experimental results of soil fractal dimension and porosity

图6 土体分形维数和渗透系数理论值与试验值对比Fig.6 Comparison between theoretical values and experimental results of soil fractal dimension and permeability coefficient

多孔介质土体渗透系数、孔隙率与分形维数之间解析关系可以为实际工程渗透系数的选取提供参考依据。先根据地勘资料提供的颗粒组成，绘制土体颗粒级配曲线，在双对数坐标轴下绘出质量百分含量与粒径之间线性相关曲线，确定土体无标度区，计算土体分形维数；再根据式（14）和地勘资料里提供的孔隙率数值，求出土体渗透系数。与传统现场试验或室内试验测定渗透系数方法相比，采用分形方法求取渗透系数避免了取样扰动和测量误差，计算简便。鉴于多孔介质土体分形维数综合反映了颗粒累积特性和孔隙大小及分布特征，可以深入研究土体分形维数在渗透破坏过程中的演变规律，为岩土工程渗透破坏非线性力学行为提供新的认知。

5 结 语

（1）无标度区是准确计算多孔介质土体质量分形维数的前提，采用多段分形维数量化土体无标度区更符合工程实际。圆砾无标度区粒径范围为0.25～20.00 mm；中粗砂无标度区粒径范围为0.25～10.00 mm；粉细砂无标度区粒径范围为0.075～0.500 mm；黏土无标度区粒径范围为0.002～0.005 mm。

（2）无标度区内土体细颗粒含量是决定土体质量分形维数的主要因素。土体颗粒粒径愈小，无标度区颗粒比重愈大，细粒含量愈高，则土体质量分形维数愈大，表明土体颗粒组成愈不均匀。

（3）渗透系数与孔隙率之间幂律函数解析公式适用于层流状态下具有分形结构的多孔介质土体无标度区。渗透系数影响因子包括分形系数、孔径大小、分形维数及流体黏滞系数，其中分形维数是主要因素。

（4）土体孔隙率、渗透系数与分形维数成反比关系。分形维数越大，表明细颗粒填充在土体骨架孔隙数量越多，相应孔隙率越小。当分形维数大于2.83时，孔隙率随着分形维数增大而减小的规律显著，在颗粒之间吸着水和薄膜水形成的黏聚力影响下，渗透系数随着分形维数增大而减小的规律较不明显。

本文针对渗透系数、孔隙率与分形维数之间非线性关系进行了理论分析和试验验证，其适用范围及推广应用仍需进一步深入研究。