2杨忠翰田东弘

(1.东北林业大学 土木工程学院， 哈尔滨 150040；2.济宁市勘测院, 山东 济宁 272000 )

饱和土作为一种两相介质，是由土颗粒及孔隙中的流体混合而成的。因此在研究饱和土体的应力-应变关系时，必须要考虑土颗粒与流体间的相互作用关系。在受到外界应力作用时，土颗粒骨架的压缩会使饱和土体中孔隙水压力升高，导致土体的有效应力降低，并使得土体的抗剪强度下降，从而造成土体的液化。对于孔隙水压力的计算，土颗粒间孔隙中的流体流动及流量变化对孔隙水压力的影响往往被忽略，在计算孔隙水压力时，计算值往往小于实际值，导致人们对于土体的有效应力及抗剪强度计算存在误差。本文基于孔隙数目-尺寸分形模型[1-2]，推导出孔隙度与分形维数之间的数学关系，并进一步基于饱和土中孔隙水流动方程及水压力方程，推导出孔隙间流量变化导致的孔隙水压力增量与分形维数及流量变化的定量关系，并进一步导出饱和土体中任一点孔隙水压力计算公式。

1 岩土材料的孔隙度与分形维数的演化关系

孔隙-固体分形模型形成于20世纪90年代，早在1989年，Neimark[3]发展了自相似多尺度逾渗系统，用于表征孔-固相界面为分形的无序离散的多孔介质；1994年，综合分形孔隙和固体的数量-尺寸分布，Perrier[4]独立地提出了一个土壤结构的多尺度模型；基于Neimark给出的定义，Perrier等[5]进一步系统地建立了孔隙-固体分形模型。该模型可用于表征具有表面分形、分形孔隙尺寸分布、分形固体质量及分形孔隙质量特性的多孔介质。因为岩土类材料从微观角度讲，是一种由大小各不相同的土颗粒堆积而成的一种多孔介质，此类材料的构造过程与分形十分相似，因此，可以基于分形模型对岩土类材料的孔隙度进行定量分析。

1.1 分形维数的计算

基于孔隙-固体分形模型给出的定义，在一个自相似系统中，分形物体的累计测量数目与单元尺寸之间通常存在幂律关系[6]，即

N(ri)=k(ri)-D,i=0,1,2,...,。

(1)

则分形维数D为

(2)

式中：N(ri)是长度为ri的单元个数；k为单位长度的初始元个数；i为迭代次数；D为分形维数。式(1)为基于分形理论的孔隙和固体的累计数目-尺寸分布的普适关系，与自相似系统中初始元维数及构造方法无关。

1.2 孔隙度的计算

1992年，Tyler和Wheatcraft曾指出土体材料中可能同时存在孔隙数目-尺寸和固体数目-尺寸分形幂律关系[7]，但当时还没有解释此现象的理论。而Perrier等[5]所提出的孔隙-固体分形模型则提供了一种明确的理论模型。郑瑛等[8]基于此模型对多孔材料的孔隙结构进行了分形描述。本文基于孔隙-固体分形模型的构造过程，对土体材料孔隙度与分形维数之间的关系进行如下推导。

Vs=(n3-m)i(R/ni)3。

(3)

则孔隙体积可表示为

Vp=R3-(n3-m)i(R/ni)3=

(4)

则孔隙度可表示为

(5)

由式(2)可得出立方体分形维数为

(6)

式中：n的直观含义为初始元样本内某一边长方向上的样本边长与最高阶次孔隙平均直径之比，即样本的测量尺度(初始元等份数)；m表示最高阶次孔隙的个数[3]。

由此可见，分形维数D可以反映样本内孔隙的大小及其分布空间的关系。

由式(6)可得，m=n3-nD，代入式(5)可得孔隙度与分形维数的数学关系式，即

φ=1-(nD-3)i。

(7)

该估计是先验过程的期望分布F0(x)和经验分布估计Fn(x|x1,x2,···,xn)的加权平均.由于经验分布函数是阶梯函数,为得到一个光滑的分布估计,用核估计代替经验分布函数Fn(x|x1,x2,···,xn),则总体的密度估计为

(8)

由式(8)可以看出，在测量尺度一定的情况下，通过分形维数即可求得土体多孔介质的孔隙度。基于此公式，可进一步利用分形维数对土体孔隙水压力进行定量描述。

2 饱和土孔隙水压力计算模型

2.1 测量球概念的提出

对于完整的岩土材料而言，孔隙水存在于其内部的土颗粒间的孔隙当中。想要通过模型反映土体内部孔隙水压力的大小及分布情况，准确地表征土体内部孔隙的分布情况是十分重要的。因此，在考虑孔隙水压力时，Okada等[9]提出了测量球的概念。测量球是一个虚拟球体空间，其球心与土颗粒球心重合，半径为土颗粒半径2倍。如图1所示，测量球是一个由若干土颗粒及颗粒间流体组成的两相体系，其优点在于可以将饱和土样的孔隙度与孔隙水压力等物理及力学参数作离散处理。根据Okada等的研究，测量球半径为对应的颗粒半径2倍时，空间内部包含土体颗粒数目适中，能够获得准确的体积应变。该模型假定固体颗粒为刚体，而孔隙中流体的体积可以发生改变，且流体的体积变化可以产生对应的孔隙水压力。

图1 Okada的测量球示意图Fig.1 Okada’s measuring sphere

2.2 孔隙水压力增量计算

饱和土体为固液两相介质，土颗粒间的孔隙全部充满流体，因此对于饱和土中任意位置测量球，其孔隙度φ可以定义为

(9)

式中：φ表示测量球孔隙度；Vv表示测量球孔隙体积；Vw表示流体体积；V表示测量球体积。将式(7)代入式(9)得

Vw=φV=[1-(nD-3)i]V。

(10)

图2表示相邻两测量球的位置关系，假设测量球内的土颗粒与孔隙为均匀分布，故两测量球的过流面积可表示为两测量球的接触面积与平均孔隙度乘积，由Darcy定律可知，单位时间内流过两测量球交界面的渗透流量为

(11)

图2 相邻两测量球位置关系Fig.2 Position relation between two adjacent measuring spheres

则A，B两测量球之间Δt时间内的流量可表示为

(12)

式中：pA，pB表示两测量球的孔隙水压力；ρ表示20 ℃水的密度；g表示重力加速度。由此可知，孔隙度可以确定相邻两测量球的过流面积，进而确定孔隙水流的交换速度。

本文基于相邻两测量球中心距离与测量球半径的数值关系，分别计算出不同情况下两测量球的接触面积，即过流面积。相邻两测量球的位置关系可分为相交或相离2类。其中相交情况又可分为两测量球部分相交或半径较小测量球处于半径较大测量球之内。当两测量球相交时，接触面积可通过两测量球A,B的相对位置进行计算。不妨假设测量球A与测量球B的半径关系为RA

其中Rt=0.5(RA+RB+LAB)。当两测量球相离，即LAB≥RA+RB时，S=0；当半径较小测量球处于半径较大测量球之内，即LAB≤RB时，S=πRB2。

在Δt时间内，一个测量球与其周围相邻测量球发生孔隙水流量交换，假设在流量交换过程中，测量球孔隙度不发生改变，则交换过程中，测量球的流量总变化可表示为∑ΔQ。根据孔隙水压力方程可得流量交换导致的孔隙水压力增量，即

(14)

式中Ew为孔隙水压缩模量。将式(10)代入式(14)可得

(15)

式(15)可视为饱和状态下孔隙水压力增量与土体分形维数之间的函数关系。实际上在饱和土体内部发生孔隙水流量交换时，土体孔隙度及渗透系数会发生微小的变化。此处假定孔隙度不变且渗透系数为常数是作简化处理。

2.3 测量球内土颗粒所受水压力计算

土体在饱和状态下，其内部完全由土颗粒及孔隙间流体组成。其中对于完全浸没在流体中的土颗粒而言，所受浮力远小于孔隙水压力，故浮力可忽略不计。假设测量球内各孔隙均为相互连通状态，即测量球内部孔隙水压力处处相同。如图3所示，只考虑单一测量球内部孔隙水压力的情况下，完全浸没在测量球内部流体中的土颗粒所受流体压力是平衡的；而边界处部分浸没的土颗粒所受流体压力为非平衡的。

图3 测量球内不同位置孔隙水压力示意图Fig.3 Schematic diagram of pore water pressure at different locations in the measuring sphere

因此，对于测量球内部任意一点的孔隙水压力，只需考虑边界处位置即可。对于边界区域的宽度，应满足测量球内最大颗粒位于边界区域时，该颗粒部分浸没于孔隙流体中而非完全浸没。因此可根据最大颗粒半径限定边界区域宽度，即0

式中：x,y,z分别表示测量球球心的三维坐标；xi,yi,zi分别表示测量球内任一土颗粒中心的三维坐标。

确定测量球边界条件后，将孔隙水压力在其作用面上积分，即可得到测量球不平衡区域内任一点孔隙水压力，即

(17)

式中：ri为边界土颗粒i与测量球接触区域所构成圆的半径。式(18)即孔隙水在饱和土孔隙间发生流量交换时，假设土体孔隙度不发生改变的前提下，推导的基于分形理论的饱和土中任意一点孔隙水压力计算模型。

3 数值算例

选取一个饱和的均质粉质黏土边坡作为算例，以验证本文所建立的孔隙水压力模型实用性及有效性。算例的边坡几何模型即网格如图4所示。模型全长为1 610 m，高505 m，其中边坡长680 m，边坡高度306.5 m。模型所选饱和粉质黏土的主要力学参数详见表1。基于本文所建立的孔隙水压力模型，算例中孔隙水压力设定为与粉质黏土分形维数相关的函数。算例中边坡内部孔隙水流速变化很小，因此可假设单位时间内土体中孔隙水流量变化为一常数ΔQ，则边坡中任一高度h处孔隙水压力可表示为

(19)

式中粉质黏土相关参数n和D可从土壤的水分特征曲线数据的拟合结果[10]中得出，代入后即可在粉质黏土边坡不同高度施加相应的孔隙水压力。

图4 算例几何模型Fig.4 Geometric model of the computation example

弹性模量/MPa泊松比内摩擦角/(°)黏聚力/kPa重度/(kN·m-3)1 5000.35343323

图5为算例中边坡饱水状态下的塑性区位移分布，其位移分布规律基本与文献[11]中相似条件下的位移分布规律一致。边坡失稳通常发生在降雨过程中，随土中含水率升高，土体逐渐趋于饱和，孔隙水压力上升致使土体有效应力降低，抗剪强度下降。边坡中饱水的土体区域持续产生塑性变形，区域内不平衡力范围逐渐扩大，最终会导致边坡失稳产生滑坡。

图5 塑性区位移分布Fig.5 Contours of plastic zone displacement

本文中孔隙水压力计算模型考虑饱和土中孔隙水渗流时流量交换产生的孔隙水压力增量，与以往基于传统孔隙水压力计算公式所做的边坡稳定分析相比，更能反映边坡产生滑动破坏的真实情况。基于分形理论得到的孔隙水压力计算公式，可通过不同类型土体对应的分形维数研究不同工况中各种土体内部孔隙水压力的数值及分布情况，对实际工程具有重要的参考价值。

4 结 语

基于孔隙数目-尺寸分形模型建立了土体材料孔隙度与分形维数之间的关系，并且进一步推导了饱和土体孔隙水压力状态方程，建立了饱和土体中任一位置土颗粒所受孔隙水压力与分形维数及孔隙间流量变化之间的关系。并将所建立的饱和土孔隙水压力函数代入ANSYS命令流中，对饱和状态下的边坡进行有限元稳定性分析，以验证建立的饱和土孔隙水压力计算模型的准确性及实用性。本文所得公式可用于饱和土体的有效应力及抗剪强度计算修正，并可应用于饱和土体宏观-微观的多尺度液相-固相耦合渗流分析。