面板坝垂直缝失效后垫层区渗流细观研究

乔 蓓，柳玉兰，江德军，黎 惠

(1.成都大学，四川 成都 610000； 2.国能大渡河流域水电开发有限公司，四川 成都 610000)

垂直缝、周边缝与水平缝等止水结构失效后的集中渗流是面板堆石坝(尤其是在深覆盖层、高寒及强震等一些复杂地质条件下的面板坝)的一个重要现象，当面板接缝止水出现裂缝，即使很小的裂缝也会导致面板坝局部或者整体渗透特性改变，严重会致使大坝产生溃坝，所以研究者们在这方面做了很多研究。王瑞俊等[1]用层流等宽运动特性，研究了面板坝极端工况的渗流问题，从而实现有限元法算面板坝的集中渗流；潘少华等[2]结合金川面板坝，采用无厚度平面单元研究止水接缝失效时渗流场规律特性，并结合解析法求解以得出准确的浸润线。但目前研究基本围绕基于连续介质的有限元法，其在求解离散颗粒物质的系统现象存在一定局限性，为此，周健等[3]运用颗粒流法定义流体和压力方程，基于流固耦合成功的运用颗粒流法模拟土体渗流特性；罗勇等[4]运用颗粒流PFC软件模拟流体中颗粒自由下落规律、渗流量及梯度在不同压力差下变化特性；曾锃等[5]结合某大坝实例，用有限元方法确定反滤层的薄弱部位及坡降，为细观研究提供边界条件，采用颗粒流PFC2D软件从细观角度研究反滤层渗透特性；游碧波等[6]采用颗粒流法模拟管涌口周围地层，基于流固耦合方法模拟管涌口管涌变化的全过程，细观研究结果与宏观研究基本吻合，研究结果验证了PFC2D在渗透破坏细观研究方面具有可行性，也为管涌现象的细观研究提供了一个新方法。与此同时，其他学者也采用可求解颗粒介质流动及固体力学大变形的颗粒流理论，证实了颗粒流理论可以用来研究渗流问题，但仍处于初级研究阶段[7-10]。

在上述研究基础上，按照某面板坝垂直缝止水失效后垫层运行实际工况，运用颗粒流理论，综合垫层孔隙率、土颗粒之间的接触关系以及密度建立细观模型，采用粗网格流法对面板垂直缝失效后垫层区渗流特性及发展过程进行模拟，分析不同水力梯度下垫层区颗粒流失量、渗透系数、流速及孔隙率等要素随时间的演化规律，研究止水失效后垫层区的渗流发展及水土相互作用过程，相关成果可为垫层区渗透稳定、管涌临界值等确定提供参考，并对实际工程中集中渗流、管涌、流土等机理研究、数值仿真奠定基础，为类似工程提供参考。

1 颗粒流基本理论

1.1 基本假定及计算原理

PFC2D软件可用来模拟土颗粒间作用及运动(1979年Cundall与Strack)[11]。其基本假定：(1)假设不考虑颗粒间孔隙率，假定用颗粒间接触处的平行板路径模拟流体的渗流路径，定义其为“管道”；(2)颗粒流模型中土颗粒用颗粒代表，流体单元并不存在，假定模型中存在可作用并储存水压的单元体，引入Domain域，定义该单元体为域，压力在计算中不断迭代更新作用于颗粒。

PFC2D存在接触刚度、滑动与粘结三种典型模型[12]，其将接触本构模型用于材料本构特性的模拟。其采用显示求解法，对所有接触的平行板路径应用液相流动方程，对于所有域交替应用压力方程，求解临界计算时步。目前，PFC2D已是求解颗粒介质的流动及固体力学的一种有效方法[13]，其计算基本原理图如图1所示。

图1 计算原理图

1.2 液相流动方程

根据颗粒流的基本假定，假设管道孔径为R，路径长为b，根据达西定理[14]计算沿垂直于平面向的单宽流量为：

(1)

式中：k为渗透系数；Pn-Pn-1为相邻域间压力差；q为单宽流量。

1.3 流固耦合方程

流固两相流模型可用连续性方程来表述，用Navier-Stokes方程表述该模型的液相流动，视流体是密度不会改变、不可压缩的。两个公式里，压力梯度均受液相与固相影响。根据这两个公式得出作用于颗粒的力fdij(i=1-np,j=x,y,z)表示为：

(2)

式中：n为局部孔隙率;βintj为颗粒和流体间摩擦系数。

1.4 流体计算单元

利用中心单元Scalar cells中心点计算压力及孔隙率，错列单元相对Scalar cells分别沿X、Y方向平移半格Scalar cells，用这种方式离散动力方程，通过Scalar cells边界计算各方向流体速度，流固耦合作用力求解过程中的颗粒半径及速度通过错列单元的平均值法计算。在模型周围加上额外虚拟单元以此来施加边界条件。

图2 流体计算单元

2 垫层料模型

2.1 工程概况

某混凝土面板堆石坝坝顶高程为3 997.00 m，最大坝高为122 m，坝顶宽为10 m，上游坝坡坡比为1∶1.4，下游坝坡综合坡比为1∶1.79，下游设置9 m宽的“之”字马道共五层，上游正常蓄水位为3 992.0 m，下游河床水位为3 897.92 m，大坝典型剖面图及分区如图3所示。

2.2 模型生成

颗粒流计算中采用理想颗粒模拟实际颗粒，将单个颗粒当做球体考虑，但在实际计算中很难建立起与实际级配曲线完全一致的数值仿真模型。考虑到一般垫层料均采用砂砾石或人工砂石等，均属于无黏性土类型，故模型建立时可不考虑颗粒之间的粘结力。利用某混凝土面板坝垫层区设计级配曲线(见图4)建立颗粒流模型，运用Fish函数随机生成土体颗粒(见图5)。垫层区密度考虑设计孔隙率、相对密度试验及室内比重实验结果等因素综合确定。

图3 面板坝典型剖面图

图4 垫层级配图

考虑到将整个坝体离散为颗粒集合体，在计算机处理能力高度发达的今天，也是不可能的，颗粒流只适用于计算小尺寸数值模型，一是由于受计算机现实条件的限制，二是粒径的上限与下限差距太大也会影响流固耦合法计算的收敛性，因此在本项目研究中只选取自坝基起25 m坝高位置开裂面板下游侧垫层区，模型尺寸0.30 m×0.15 m，此外，考虑对垫层的细观模拟主要关心的是垫层发生渗透破坏时的破坏机理，而不是具体量值的大小，因而模型尺寸的选取是合理的。为了提高模型的计算速度，粒径离散化为1～2 mm、1～3 mm、3～6 mm、6～8 mm、8～10 mm、10～12 mm共计六部分。采用粒径区段半径下限、上限及区段颗粒占比确定不同粒径段颗粒个数，并在区域内按照半径扩展法随机生成孔隙率模型，颗粒集合体在自重应力下平衡，颗粒流模型中包含共450个流体单元，共25 270个颗粒数。实际计算中，为消除颗粒间及与墙体间由于过多的重叠而引起较大的不平衡力，故将不平衡力率的均值及极值控制在1×10-4之内。

图5 计算模型

(3)

(4)

2.3 参数及边界条件选取

模拟垫层料这种人工砂石料或者加工后的砂砾石等无黏性土时，一般选择接触刚度模型，根据相关学者研究工作[3,7-8,15]及本工程实际情况确定模型参数如表1。

表1 垫层料参数表

本文模拟大坝在正常蓄水位工况下的受力条件，在模型上游侧施加不同的水头压力，下游侧为透水边界(压力为零)，上下两侧加墙体不透水边界，则模型两侧产生水头压力值，并根据宏观尺寸与细观模型尺寸比例计算得出模型两侧的水头压力值。根据图6有限元计算结果显示某面板坝垫层料在单缝失效时上下游水头差约为60 m，因垫层实际水平尺寸为4 m，沿渗径向细观模型尺寸为0.3 m，根据相似性原理计算，应施加细观模型压力差约45 kPa。考虑实际运行过程中的压力变化，在模型上游侧分别施加40 kPa、45 kPa、50 kPa三组水头压力，从而使上下游形成压力差，编写FISH函数，经过多次迭代计算，并利用流体单元记录全过程颗粒流失量、流速及孔隙率等要素随时间及压力变化的演化过程及规律。

图6 单缝失效时垫层区总水头变化曲线

3 数值模拟结果分析

3.1 渗流要素与计算时步间演化规律

对不同压力梯度下颗粒流失量随计算时步的演化规律进行分析。图7为压力差分别为40 kPa、45 kPa、50 kPa时模型下游侧颗粒流失量随时步的变化趋势，从模拟数据可以看出，在一定压差范围下，颗粒流失量整体较小，最高损失量在5%以内。这三种压力差下，开始计算阶段，颗粒流失量几乎为0，处于静止状态，随着时间的推移，颗粒流失量不断增大，直至流失量达到稳定状态，所以当作用一定时间的渗透力后，垫层才会慢慢产生侵蚀。同时，从不同压差颗粒流失趋势可以看出，颗粒流失量随着压力差(即水力坡降)的增大而增大，符合一般规律。

图7 颗粒流失量与时步曲线

垫层区集中渗流是水土相互耦合作用的过程，在渗透压力的作用下，土体骨架中的颗粒随水流产生运动，土体颗粒发生重组，流速也会随着渗透压力及时间变化。图8、图9给出了平均渗流流速及流量与时步的关系，可以得出：初始施加压力水头不久，流速及渗流量迅速增大，直至达到峰值，随着计算时步的增加，流速先是小幅减小，最终达到稳定状态。计算时初始颗粒流失量为0，这是因为初始施加压力时水流流速较小，不足以致使颗粒移动，仍然保持静止状态，当渗流流速达到颗粒起动速度时，土体骨架中的小颗粒随水流出现缓慢运行，运动启动后，土颗粒运动速度将低于启动速度，并持续保持颗粒移动。

图8 流速与时步曲线

图9 渗流量与时步曲线

图10、图11给出了不同压力差垫层区孔隙率及渗透系数随计算时步的演化趋势。可以得出：孔隙率随着计算时步的推移而逐渐增大，最终达到平衡，在不同的压力作用下平衡状态的孔隙率相差不大，与初始状态相比，在三种不同压力下孔隙率增大均约11%；同时，渗透系数也随着计算时间的推移呈现缓慢逐渐增大趋势，增大比例约20%左右，其中相比孔隙率，渗透系数随着压力差的变化较大。产生上述现象是因为：垫层在渗透力作用下，细颗粒流失导致垫层孔隙率逐渐增大，在未产生明显渗透破坏前，通过颗粒的重新排序、重组，达到新的平衡状态，孔隙率的增大进而导致渗透性增大。

图10 孔隙率与时步曲线

图11 渗透系数k与时步曲线

3.2 流速随压力差变化规律

图12给出了当渗流达到稳定时垫层区渗流流速与压力差的关系，可以得出：在水力坡降较小时，渗流缓慢，流速小，处于层流状态，符合达西定理，两者线性正相关；但随着水力坡降的不断增大，渗流流场呈现一定的紊流状态，流速与压力差逐步表现出非线性关系。

图12 流速随压力差变化曲线

3.3 不同时步流速矢量图

图13跟踪记录压力差45 kPa下平均流速矢量，初始时刻，垫层料流速较小且基本成均匀分布，随着时间推移，流速不断增大，当达到一定时刻，下游侧流场发生偏移，细颗粒开始流失，但是在后续时步，流场基本不变，所以在该压力差下，垫层区并未出现局部流速突增或者模型上下游形成贯穿通道，其并没有产生渗透破坏。同时，在50 kPa及55 kPa压力差下，也未出现明显的渗透破坏。

图13 压力差45 kPa下平均流速矢量图

4 结 论

本文按照某面板坝垂直缝止水失效后垫层运行实际工况，运用颗粒流理论，综合垫层孔隙率、土颗粒之间的接触关系以及密度建立细观模型，采用粗网格流法对面板垂直缝失效后垫层区渗流特性及发展过程进行模拟，通过写入FISH函数并利用流体cell记录不同压力差下垫层区颗粒的流失量、流速、渗透系数及孔隙率等要素规律，分析不同水力梯度下垫层区颗粒流失量、渗透系数、流速及孔隙率等要素随时间的演化规律，研究止水失效后垫层区的渗流发展及土水相互作用过程，得出如下结论：

(1) 该面板坝止水失效后，从各水力要素随时间演化规律看，垫层区未发生明显渗透破坏，颗粒流失量总体较少，均在5%以内。

(2) 渗流关键要素颗粒流失量、渗透系数、流速及孔隙率等均随压力梯度呈正相关，且随着计算时步的增加，各渗流要素均有增大，但增幅大小各异，最终达到平衡状态。

(3) 采用基于颗粒流研究面板坝垂直缝失效后垫层区的渗透特性具有可行性，且从细观出发，采用颗粒流方法对渗透特性机理等进行研究能够更加直观解释渗流发展及水土相互作用过程，相关成果可为垫层区渗透稳定、管涌临界值等确定提供参考，并对实际工程中有限元法模拟渗流问题奠定基础。