堤(坝)基层状粗粒土发生渗透变形的颗粒运移规律

黄海均,严新军*,毛海涛,王正成

(1.新疆农业大学水利与土木工程学院，乌鲁木齐 830052；2.河北工程大学水利水电学院，邯郸 056002；3.重庆三峡学院土木工程学院，万州 404100)

堤(坝)基中常出现互层的粗粒土透水层，是渗透薄弱区域，极易发生渗透破坏，具有代表性的粗粒土互层主要由砂砾石层、细砂层和砂砾石层组成，渗透时各土层微观的颗粒运移规律对于揭示堤(坝)基渗透变形和破坏机理至关重要[1]。因此，需要深入研究层状粗粒土渗透变形颗粒运移规律。

近年来，中外学者针对管涌发生机理开展了大量的研究工作。刘昌军等[2]、梁越等[3]利用室内砂槽模型试验对单、双层堤(坝)基管涌的发生及发展过程进行了模拟，将管涌破坏的发生及发展过程划分为4个阶段；贾恺等[4]利用自行设计试验装置对双层堤(坝)基管涌通道扩展机制进行试验，得出了双层堤(坝)基管涌通道是否扩展取决于通道内冲刷水流、边壁渗流力和砂层颗粒条件等因素。Shamy等[5]运用颗粒流方法，对上游水位上升时层状堤(坝)基发生管涌导致水工建筑物出现较大位移和沉降的过程进行了模拟；文献[6-7]给出了单一砂层堤(坝)基上临界水头的理论解答。不难发现，目前中外对此类问题的研究主要集中于管涌临界渗透破坏条件的判别和单双层堤(坝)基管涌通道发展机理方向，而对此类堤(坝)基管涌发展过程中颗粒运移特点研究极少，因此很有必要对层状覆盖层开展相关研究。

针对典型堤(坝)基的层状粗粒土进行研究，利用颗粒流软件PFC3D建立对应数值模型，模拟堤(坝)基颗粒运移特点及渗透变形发展过程并分析颗粒流失量、下沉量等参数的变化，为进一步认识层状覆盖层堤(坝)基渗透破坏提供参考。

1 颗粒流程序计算原理

1.1 流体相方程

在固液两相体系中，相对于固体颗粒间孔隙体积变形量，液体的体积变形量几乎可以忽略不计，因此可将颗粒孔隙中的流体视为密度不变的不可压缩流体，在求解渗流场的每一时步中，满足流体运动的连续性方程和平均Navier-Stokes[8]方程：

(1)

(2)

式中，n为颗粒孔隙率；t为时间；为拉普拉斯算子，为流体的速度矢量；ρf为流体密度；p为压力梯度；τ为流体黏性应力张量；g为重力加速度矢量；fint为单位体积内颗粒与流体的相互作用力矢量。

1.2 固体相方程

对于固体颗粒，考虑流体作用时，其主要受到流体对颗粒、颗粒间相互作用及外应力等作用力，并通过牛顿第二定律来实现模拟颗粒的运动，此时颗粒运动方程为

(3)

(4)

式中，mp为颗粒的质量；vp为颗粒速度矢量；fg为重力加速度矢量；fc为接触c(c=1，2，…，n)处颗粒间的接触力；fd为颗粒在流体作用下受到拖曳力，包含浮力及固液之间的相互作用力；Ip为颗粒的转动惯量；wp为颗粒的转动速度矢量；rc为颗粒间接触处指向颗粒中心的方向矢量。

1.3 流固耦合方程

固体颗粒在单个流体网格内的流动情况如图1所示。

图1 流体单元渗流示意图Fig.1 Seepage diagram of fluid unit

流体单元网格体积为ΔV=ΔxΔyΔz，假定流体单元网格内包含np个颗粒，则流体单元网格内的孔隙率

(5)

式(5)中，dpi(i=1，2，…，n)为颗粒粒径。

考虑流体单元网格内土体颗粒的受力平衡条件，可得单个固体颗粒受到的作用力

(6)

式(6)中，fdij为单个颗粒所受作用力分量；pj为流体压力梯度分量。

对于稳定无分叉流，假设颗粒流体间的相互作用力只来源于压力梯度，即：

fintj=npj

(7)

将式(7)代入式(6)，可得单个固体颗粒受到作用力的一般表达式：

(8)

当雷诺数为1～10的层流时，German将水力半径理论公式代入Darcy公式得压力梯度方程为

(9)

对于雷诺数大于10的紊流，Ergun[9]基于水力半径理论提出了Ergun公式：

(10)

1.4 流固耦合过程

在颗粒流程序中，整个流体区域被划分为若干个固定大小的流体计算单元，采用基于半隐式算法(SIMPLE法)结构的计算流体动力学(CFD)计算技术对液相流动方程进行求解。在计算过程中，首先根据力-位移定理计算得到颗粒间的接触力，而后将该力与流体对颗粒的拖曳力一起代入下一时步颗粒运动方程中，以此确定颗粒新的位置、速度及接触力等。同时颗粒的运动又将引起流体单元内压力、流速、孔隙率的变化，进而影响到颗粒间接触力及流体对颗粒拖曳力的变化，可见流固耦合模拟过程为一个动态循环平衡过程，其计算过程如图2所示。

图2 PFC3D流固耦合循环计算过程Fig.2 PFC3D fluid-solid coupling cycle computation process

2 数值模型建立

2.1 建模依据

中国西南和西北的河谷覆盖层在纵向结构上多具有强、弱透水层交替层叠的特点，各层岩性、物理力学特性等差异较大，较有代表性是砂卵砾石和细砂互层结构如图3(a)，各层取样如图3(b)、图3(c)所示。

图3 典型粗粒土及粗细相间的层状覆盖层Fig.3 Typical coarse grained soil and coarse and fine layered overburden

典型砂卵砾石层和细砂层的基本特征及透水性等级如表1所示。

表1 粗细粒土基本特征及透水性等级Table 1 Basic characteristics and permeability grade of coarse and fine-grained soil

根据上述典型的层状粗粒土层，运用离散元软件PFC3D来建立渗流水-土耦合模型，厘清颗粒运移规律。

2.2 数值模型生成及细观参数标定

2.2.1 数值模型生成

为了更准确地对实际堤(坝)基中层状粗粒土结构的渗透变形过程进行模拟，模型生成的颗粒必须能反映实际土体颗粒级配关系，实际土体中颗粒粒径范围(0.075～50 mm)较宽，如按照实际土体颗粒级配建立模型，将生成数量十分巨大的颗粒，导致计算效率降低甚至无法完成计算。因此，需对实际土体级配进行一定的近似处理，本文数值模型中颗粒最大最小粒径分别取为10 mm和0.25 mm，由工程地质图3及表1可知，上部土层为间断级配的砂砾石土层，以2 mm作为粗细颗粒的区分粒径，同时保持生成颗粒质量百分比与实际土体粗细颗粒百分比基本一致，近似处理后的土体级配如表2所示。应注意的是，该颗粒粒径处理方法可能对试验精度造成一定的误差，但对渗透变形中颗粒运移过程的研究是可行的。此外，为降低颗粒模型中总的自由度，这里运用土工离心机试验中的相似性原理，建立长×宽×高为120 mm×60 mm×60 mm的颗粒流模型如图4所示。同时，为使建立的模型与实际土体的工程地质分层情况基本一致，利用PFC3D内置的fish语言编写程序以生成不同土层，模型从上至下依次为：上部砂砾石层厚26 mm，中间细砂层厚6 mm，下部砂砾石层厚26 mm。此外，为防止上覆土层与模型上部墙体发生接触冲刷，除管涌口外，在墙体上部边界生成颗粒粒径为2 mm规则排列的固定颗粒边界。

表2 PFC数值模型土样参数Table 2 Soil sample parameters of PFC numerical model

图4 堤(坝)基中多层粗粒土结构颗粒流模型Fig.4 Particle flow model of multi-layer coarse-grained soil structure in embankment (dam) foundation

2.2.2 细观参数标定

为使所建模型能更有效的反映实际渗透变形过程，需要准确找出宏观材料与细观颗粒之间物理力学属性的联系。在利用PFC3D来模拟渗流破坏时，颗粒运动规律主要受到细观参数摩擦系数的影响[10]，因此在数值模拟之前，必须对颗粒的摩擦系数进行标定，但颗粒宏-细观参数之间存在较大的差异。研究表明，在散粒体材料中，材料的内摩擦角是其摩擦特性的综合体现[11]，因此采用“反演模拟法”[12]来对颗粒的细观摩擦系数进行标定，室内试验测得砂砾石和细砂的宏观内摩擦角分别为35°10′、26°30′，然后利用PFC3D进行三轴数值模拟试验(图5)，通过不断调整颗粒细观摩擦系数来获取颗粒的细观内摩擦角，将得到的细观内摩擦角与宏观内摩擦角进行对比，直到两者接近或相等为止，最终得到粗砾和细砂的摩擦系数为1.25与0.79。模型具体细观参数如表3所示。

图5 三轴数值模拟试验Fig.5 Triaxial numerical simulation test

表3 PFC模型具体细观参数Table 3 Specific meso-parameters of the model

参照表2及表3中模型细观参数，利用PFC3D中fish语言编写函数生成相应土层颗粒，各土层之间赋予线性接触模型，采用半径扩大法生成土颗粒，先将颗粒粒径缩小，再将颗粒粒径放大到指定的孔隙率，而后对土颗粒施加重力并进行循环计算，直至土颗粒间不平衡力较小或消除为止，以确保土层达到稳定状态。

2.3 流体网格划分及监测单元设置

PFC3D中运用固定粗糙网格法对流体进行处理，需将流体区域划分为若干流体单元网格，且应确保所划分的流体网格中均包含一定数量的颗粒。故考虑将整个模型沿长、宽、高方向分别划分为8、4、4份，单个流体单元网格尺寸为15 mm×15 mm×15 mm。再进行流体计算之前，需对边界条件进行设定，模型左侧设定为施加水头边界，水头压力大小为100 kPa，上部、底部及侧壁均设定为刚性不透水非滑移边界，将模型上部X=12 mm、Y=39 mm、Z=24 mm处设置为零压力流体出流口，出流口长宽为3 mm×3 mm[见图6(a)]，该流体单元网格模型图详见图6(b)。应注意的是，在左侧边界施加水头之前，为避免流域内形成紊乱的流场，应首先将颗粒固定，然后施加水头，直至流场趋于稳定，此后释放颗粒，让固体颗粒与流体进行作用。

图6 计算边界及流体单元Fig.6 Computational boundary and fluid unit

在模型开始计算之前，为清晰记录各土层不同部位颗粒流失的情况，分别在各土层设置以下监测区域，上部砂砾石层及细砂层从压力上游端至下游端分别编号为S1、S2、S3、S4及Z1、Z2、Z3、Z4，如图7所示，由于细砂层对下部砂砾石层起阻挡作用，导致下部砂砾石层颗粒很难发生流失，因此不对下部砂砾石层中颗粒流失进行监测。

图7 模型不同监测区域划分Fig.7 Division of different monitoring areas in the model

3 数值试样结果分析

3.1 上部砂砾石层渗透破坏过程

3.1.1 上部砂砾石层细颗粒迁移过程

该粗粒土层状结构的上部砂砾石层中细颗粒迁移过程如图8所示。

从图8可以看出，因上部砂砾石为管涌型土，粗颗粒之间存在一定的孔隙，在初始水头压力作用下，上部砂砾石层中细颗粒从上游端逐渐向管涌口运移，上部砂砾石层中细颗粒在管涌口陆续流失，但并未出现颗粒整体流失及细砂层颗粒侵入上部砂砾石层中的现象 [见图8(a)]。在计算时间步达到30万步时，上部砂砾石层粗颗粒中细颗粒持续从上游端向管涌口处移动，导致在上游水头S1区出现较大孔隙，在细颗粒移动过程中，有部分细颗粒因遇到较小孔隙而被阻挡形成阻塞区，此时，细砂层Z3和Z4区已有部分颗粒进入上部砂砾石层中，且由于S1区有较大孔隙出现，细砂层Z1区部分颗粒也开始进入上部砂砾石层中[见图8(b)]。随着计算时间步达到50万步时，上部砂砾石层中细颗粒沿着粗颗粒间孔隙不断从上游端向管涌口处移动，细砂层中颗粒进入上部砂砾石层中含量进一步增加，有部分颗粒已进入上部砂砾石层中部，以至于细砂层出现较大变形[见图8(c)]。

图8 上部砂砾石层中细颗粒迁移过程Fig.8 Migration process of fine particles in upper gravel layer

3.1.2 上部砂砾石层不同区域细颗粒流失量

上部砂砾石层中不同区域细颗粒流失量随计算时间步的变化曲线如图9所示。颗粒流失量定义为该区域细颗粒流失的总体积与水头压力施加前该区域细颗粒总体积之比。

图9 上部砂砾石层不同区域细颗粒流失量情况Fig.9 Fine particle loss in different areas of upper gravel layer

由图9可知，在数值试验过程中，上部砂砾石层上游端S1区细颗粒流失量最大，约20%，导致S1区粗颗粒间有较大孔隙形成。S2、S3区细颗粒流失量相对较少，当时间步增加至20万步左右时，S2、S3区颗粒流失量分别增至3.5%、3%，而后直至计算结束，流失量分别降低至2.5%、1%，这是由于细颗粒从上游端流向管涌口处时逐渐在S2和S3区形成阻塞区所致。当时间步计算至25万步时，S4区细颗粒流失量增至5%左右，之后随时间步的增加，S4区颗粒流失量几乎没有变化，这是因为在试验初期管涌口颗粒不断流失，随着时间步的增加，上游端细颗粒持续流向管涌口，同时，细砂层颗粒陆续进入上部砂砾石层中，导致S4区细颗粒流失量逐渐减少直至几乎不变。

通过对上部砂砾石层中细颗粒迁移现象及流失曲线分析可知，在渗透力作用下，上部砂砾石层中细颗粒在粗颗粒孔隙间移动而后被带出土体外，逐渐形成稳定的管涌通道，属于典型的管涌破坏。

3.2 细砂层渗透破坏过程

3.2.1 细砂层颗粒迁移过程

该层状结构中的细砂层中颗粒迁移过程如图10所示。

图10 细砂层中颗粒迁移过程Fig.10 Particle migration in fine sand

从图10(a)可清晰观察到，在施加水头压力初期，细砂层颗粒几乎没有发生。如图10(b)所示，随着计算时间步增加至30万步，细砂层中上游端Z1区及管涌口附近Z3及Z4区颗粒已部分侵入上部砂砾石层中，甚至在管涌口正下方细砂层Z4区已有少许颗粒进入上部砂砾石层中部，这是由于上部砂砾石层中细颗粒的逐渐流失，给予了细砂层进入上部砂砾石层的空间。当时间步进一步增加至50万步过程中，细砂层Z1、Z3及Z4区处颗粒进入上部砂砾石层中越来越多，且在细砂层Z4区已有一些颗粒流失至管涌口附近，随着颗粒持续流失，导致细砂层逐渐出现小范围的变形[见图10(c)]。

3.2.2 细砂层不同区域颗粒流失体积

细砂层不同区域进入上部砂砾石层颗粒流失体积随计算时间步的变化曲线如图11所示。

图11 细砂层不同区域颗粒流失体积变化Fig.11 Variation of particle loss volume in different areas of fine sand bed

由图11可知，在数值试验计算过程中，细砂层在管涌口下方Z4区处颗粒持续流失并且流失颗粒含量最多，约3.5×10-8m3，其余区域颗粒流失量相对较少。在数值试验初期，仅有曲线S4随时间步的增加而增加，而其余曲线无明显变化，表明细砂层最先在管涌口正下方Z4区处出现颗粒流失。随着计算时间步的增加至30万步，曲线Z1、Z2及Z3均随着计算时间步的增加而增加，说明Z1、Z2及Z3区颗粒都已发生流失，且在此过程中，可见Z2区颗粒是最后发生流失的，这是由于上部砂砾石层S2区有阻塞区形成，阻碍了细砂层颗粒的进入。当计算时间步增加至50万步时，曲线Z1随时间步缓缓增加，对应细砂层Z1区颗粒流失量逐渐增多，这是因为上部砂砾石层中上游端S1区细颗粒逐渐流失，以至形成较大孔隙，给予了Z1区颗粒进入上部砂砾石层空间。

通过对中间细砂层中颗粒迁移现象及流失曲线分析可知，细砂层中颗粒最先在管涌口正下方Z4区发生流失，其余区域颗粒流失相对较晚，且颗粒流失量均随着计算时间步的增加而增加，导致细砂层出现小范围的变形。

3.3 上部砂砾石层下沉量

如图12所示为上部砂砾石层下沉量随计算时间步的变化曲线，下沉量定义为上部砂砾石层颗粒垂直向下方向运动距离总和。从图12中可以看出，随着计算时间的增加，上部砂砾石层的下沉量逐渐增加至0.012 m，这是由于细砂层中颗粒不断流失所致，这种情况下当上部砂砾石层细颗粒含量流失达到一定程度，随着管涌通道的贯通而发生破坏，将对上部建筑物产生重大危害。通过颗粒流模拟，上部砂砾石层及细砂层所表现出的渗透破坏现象，这与刘杰[13]对间断级配土渗透破坏机理试验研究所得结论相吻合。

图12 上部砂砾石层下沉量Fig.12 Subsidence of upper gravel layer

4 结论

在PFC3D细观参数标定时，采用“反演模拟法”，较准确地建立了材料宏观参数与颗粒细观参数间联系，有效的模拟了多层粗粒土层渗透变形的发展过程，获得了堤(坝)基渗透变形中颗粒流失及下沉量等情况，为此类地基的渗透破坏的研究提供了一定参考，并得到以下结论。

(1)渗透破坏主要发生在上部砂砾石层中，随着渗透破坏的持续发生，逐渐影响下部土层；作为渗流出口的上部砂砾石层，该层中细颗粒在粗颗粒孔隙间移动而后逐渐被带出土体外，逐渐形成稳定的管涌通道，具有明显的溯源性。

(2)随着计算时间的增加，上部砂砾石层的下沉量逐渐增加，当上部砂砾石层细颗粒含量流失达到一定程度，堤(坝)基发生破坏，将对上部建筑物产生重大危害。

(3)中间细砂层在上部砂砾石层管涌破坏后，其颗粒最先在管涌口正下方Z4区发生流失，其余区域颗粒流失相对较晚，且颗粒流失量均随着计算时间步的增加而增加，导致细砂层出现小范围的变形出现。

(4)由于几何条件不满足，在上部土层不发生明显破坏的情况下，不会对深层的砂砾石层颗粒运移造成影响。

综上，上部强透水层是渗透破坏的关键区域，在实际工程中应重点处理；深层次的强透水层对堤坝基渗透破坏的贡献较小，在上部土层不发生破坏的情况下，深部土层一般不会发生颗粒流失现象。