刘昌军

(中国水利水电科学研究院,北京 100038)

尾矿坝的坝体是由透水性强弱不等的粗细尾矿砂堆筑而成,而坝坡浸润线是尾矿坝的生命线,它是直接影响坝体安全的一个非常重要的因素。我国多数尾矿坝采用上游法筑坝工艺,其显著特点是随着尾矿堆筑高度的不断加大,尾矿坝坝体内浸润线也不断提高,导致坝面溢出水位升高,易引起坝体渗流破坏,危及坝体稳定性。因此需要对尾矿坝设计一定的排水系统以降低坝体浸润线的高度。

在尾矿坝设计与施工中,虹吸井、导渗盲沟和辐射井等排渗措施都有广泛的应用,但虹吸井和导渗盲沟往往径向尺寸很小而纵向尺寸很大,数量众多,且空间分布复杂。因此在用有限元模拟数量众多的虹吸井和导渗盲沟作用下的渗流场分布就非常困难。目前堤坝渗控分析中对密集排水孔的模拟较多[1-5]。对尾矿坝虹吸井、导渗盲沟以及辐射井等排水措施的排水效果和渗流场进行数值模拟研究的文献尚不多见。

为了解决难以精细模拟密集虹吸井和导渗盲沟的排水效果的难题,笔者提出一种适合求解虹吸井、导渗盲沟和辐射井等布置有复杂排水系统的尾矿坝三维渗流场的虹吸井子结构方法,该方法较好地解决了尾矿坝排水系统的精细模拟问题,可用于多种排水措施共同作用下多种排水方案的计算比较。以某尾矿坝为例,采用虹吸井子结构法研究了密集虹吸井和导渗盲沟作用下的渗流场分布,并对不同排渗方案的排渗效果进行对比分析,给出了该尾矿坝排水系统的合理布置方案。

1 渗流场求解的关键技术

1.1 基本方程

非均质各向异性多孔介质的稳定饱和渗流问题的控制方程为

式中,kij为二阶对称的达西渗透系数矩阵;h为水头;h1为已知水头函数;ni为渗流边界面外法线方向余弦,i=1,2,3;Γ1为第一类渗流边界条件;Γ2为第二类渗流边界条件;Γ3为渗流自由面;Γ4为渗流逸出面;qn为法向流量,以流出为正;x3为z方向的位置高程。

1.2 渗流场求解的有限单元法

对式(1)～(5)可采用固定网格的有限单元法进行求解,根据变分原理,求解泛函数和支配方程为

式中:∏(h)为泛函数;K,P,F分别为饱和区的传导矩阵、节点水头列阵和已知节点压力列阵;Ω为积分域。

1.3 改进的截止负压法

文献[6]提出了有自由面渗流场的求解方法——截止负压法,文献[7]在文献[6]的基础上对截止负压法进行了改进。本文采用文献[7]中改进的截止负压法对渗流场进行求解。改进的截止负压法的详细公式推导见文献[7]。

2 虹吸井工作原理及其渗流行为

2.1 虹吸井的组成及工作原理

虹吸井排渗系统主要由虹吸井、观测井、水封槽和虹吸管路组成[8],见图1。其中,虹吸井是关键的降水设施,其成井质量要求较高。一般虹吸管内水头差约8m左右。水封槽的作用在于保证虹吸系统的绝对真空度,即在虹吸系统启动后应保证排水管出口始终位于一定的水位高度下。

2.2 虹吸井的渗流行为

虹吸井是靠虹吸作用实现自流排水。虹吸井渗流行为在算法上分为2种情况:一是虹吸井穿过自由面;二是虹吸井不与自由面相交。虹吸井穿过自由面又分为2种情况:一是虹吸井达到稳定状态,井内水位和水封槽水位相同;二是虹吸井内水位大于水封槽水位,虹吸井内出水量较大,虹吸井以恒定流量出水。有自由面穿过时,虹吸井内边界面上的渗流行为如图2所示。

图1 虹吸井排渗系统布置图

图2 虹吸井渗流行为

a.如果自由面水位低于cc′,则整个虹吸井起不到排水作用,虹吸井失效。因此可以在每个虹吸井水封槽水位高程处虚构一个数学开关器,详见文献[9]和文献[10]。虹吸井计算时,先假设开关器打开,孔内边界全部作为可能逸出边界,每一步迭代后对孔内节点压力进行甄别。

若虹吸井底面孔口节点压力大于零,认为虹吸井全部位于非饱和区,此时虹吸井完全失效。其边界条件数学表达式为

若虹吸井cc′处节点压力大于零,认为虹吸井水封槽高程以下部分位于饱和区,此时边界为定水头边界。其数学表达式为

式中:x3c为c点处z方向的位置高程。

若虹吸井cc′处节点压力小于零,且底面节点压力大于零,认为虹吸井自由面位于水封槽底部,此时虹吸井失效。

b.如果自由面水位位于bc或b′c′内,且逸出流量小于虹吸井设计恒定流量 Q,则 ac和a′c′面水头为定水头,水头值为水封槽的高程,此时 bc或b′c′面为隔水边界,因为虹吸井起动直至水位稳定,虹吸井内水位和水封槽高程保持一致。边界条件数学表达式为

c.如果自由面水位位于bc或b′c′内,且逸出流量大于虹吸井恒定出水量,则虹吸井以恒定抽水量抽水,此时虹吸井边界条件为流量边界条件(但实际工程中,虹吸井设计恒定流量一般大于虹吸井内逸出流量)。其虹吸井内边界条件数学表达式如下:

虹吸井在改进的截止负压法计算中,作为二类边界(可能逸出边界)来处理,即初始迭代是把孔壁内节点高程作为边界条件,然后逐步判别各节点是在饱和区还是在非饱和区。根据逸出点位置进行虹吸井孔壁逸出流量计算,计算逸出流量并与虹吸井设计恒定流量进行比较。然后根据上述情况转化虹吸井边界条件,进而实现虹吸井渗流行为的真实数值模拟。

对于导渗盲沟和辐射井的处理,同样采用上述子结构方法,子结构内边界按可能逸出边界处理。

2.3 虹吸井子结构法的计算原理

按照上述虹吸井渗流行为,结合改进截止负压法原理,将子结构作为主网格整体结构的一部分一起来考虑[2,3,9,11],则虹吸井子结构总传导矩阵及相应流量列阵写成分块形式,如式(12)所示。

式中:Kii为子结构内部未知节点水头相互作用传导矩阵;K ib,K bi分别为子结构内部未知节点水头与出口未知节点水头之间相互作用的传导矩阵;K bb为子结构出口未知节点水头相互作用传导矩阵;p i,p b分别为子结构内部和子结构出口未知节点水头列阵;Fi,Fb分别为子结构已知节点水头对内部未知节点水头和出口未知节点水头贡献流量(包括内部源汇项及非零流量边界的贡献流量)的列阵。

由式(13)～(14)可求得相应未知节点出口传导矩阵及相应流量贡献列阵K*bb,F*b如下:

用K*bb,F*b参与总平衡方程系数及常量矩阵的组装,求解 pb后,回代求解 pi。

2.4 子结构边界条件的处理

对于子结构内部已知水头节点,其节点压力为已知,因而其每步迭代求得的压力增量为零,故可将渗透矩阵中该节点对应的对角元素乘以一个大数即可。

对于子结构内部可能逸出面上的节点,如果发现其节点压力大于零,可以采用处理截止负压的方法进行类似处理,即将渗透矩阵中该节点对应的对角线元素置为一个大数 λ,同时将其对应的节点不平衡力赋为-λp(p为该节点的压力),因此可以迫使p=0。

2.5 虹吸井子结构的剖分模式及剖分方法

在尾矿坝工程中,为精确反映虹吸井的三维尺寸效应,虹吸井通常是在特定区域内的母单元中进行剖分。虹吸井的布置方式和模式因工程特点和排水要求各不相同,在虹吸井子结构法中,虹吸井的子结构形式和子结构网格剖分模式参照文献[2,9,11]中介绍的方法。笔者采用IDL语言编写了有限元数值模拟软件GWSS,该软件提供了有限元网格剖分、虹吸井子结构网格剖分、渗流计算和渗流结果等值面(线)的显示等功能。在虹吸井子结构网格剖分方面,该软件实现了不同方向的虹吸井子结构单元网格的二次剖分和弯管子结构单元的剖分。其剖分原理和剖分步骤如下:①根据整体网格单元和事先控制的虹吸井和导渗盲沟的超单元,确定需要剖分子结构单元的母单元的节点信息和单元信息,同时搜索其相邻子结构单元的母单元串。②找出最小单元截面周长l,虹吸井等效母单元边长 a由公式a=πd/l(d为虹吸井的直径)算出,进而计算母单元中子结构单元的局部节点坐标,形成局部节点和单元信息。③根据母单元新增单元和站点信息组装子结构信息,进而得到整体网格和子结构网格的单元信息。

3 工程应用

3.1 工程概况

河北某尾矿坝,尾矿砂为中等颗粒,干密度为1.4 t/m3。该尾矿坝初期设计为碾压土坝,坝高20m,坝顶宽度为3.5m,坝顶高程为450m,坝底高程为430m。外坡坡比为1∶2.5,内坡坡比为1∶2.0。现初期坝顶已被碎石覆盖,碎石覆盖后的坝顶高程为459.6m,外坡坡比为1∶1.36。初期坝坝坡多处有清水渗出。

3.2 工程地质条件

可将尾矿坝地层分为3个部分,即天然地层、初期坝和尾矿堆积层。天然地层表面为残坡积含砾粉质黏土,其下部为粉土及粉质黏土充填的碎石土层(强风化层)。初期坝为碾压土坝。由于尾矿沉积的分选性不强,故砂性尾矿中普遍存在黏性尾矿夹层(3～10cm厚),且在同平面上不连续,呈透镜体状,但从整体上看仍存在一定规律:垂直层序上存在上粗下细的规律;水平层序上由堆积子坝向库区尾矿砂由粗变细。根据室内筛分试验判定该区尾矿砂大部分为尾矿中砂。

3.3 模型参数及计算工况

表1 尾矿坝各地层渗透系数 m/s

3.4 模型范围及网格划分

计算区域采用八节点六面体单元进行剖分,共剖分节点13120个,单元10500个,其三维网格见图3,尾矿坝虹吸井和导渗盲沟等排水系统采用子结构进行剖分,其三维网格见图4、图5。

3.5 计算结果分析

为节省篇幅,本文只给出了工况2和工况4计算结果。图6和图7分别为工况2和工况4顺河向剖面水头等值线,表2给出了不同工况设置虹吸井后虹吸井出水量。

从计算结果可以看出:①各工况渗流场的水头分布规律合理,水头等值线形态、走向和疏密程度都

图4 虹吸井子结构网格(局部放大)

图5 导渗盲沟子结构网格(局部放大)

表2 各工况虹吸井出水量 m3/d

尾矿坝各地层渗透系数见表1。工况1和工况2分别为正常蓄水位和设计洪水位下的渗流场分布,工况3和工况4分别为下游侧虹吸井和上游侧虹吸井失效的渗流场分布。准确反映了相应区域防渗和排水渗控措施特点、渗流特性和边界条件,计算区域内虹吸井的排水效果得到了精细模拟,其渗控效果得到了准确反映。②采用虹吸井子结构方法可以精细模拟每个虹吸井的排水作用,可以计算单个虹吸井的排水量,如工况2,上排单井平均出水量为108.500m3/d,下排虹吸井平均出水量为7.590m3/d。③从工况3和工况4的计算结果看,下游侧虹吸井完全失效后,地下水位逸出点高程为510m,而上游侧虹吸井失效后,地下水位逸出点高程为523m。和工况2相比,下游侧和上游侧虹吸井失效后尾矿坝逸出点高程都相应增加,严重影响坝体稳定性,因此建议采用2排虹吸井设计方案,并确保虹吸井施工质量,防止虹吸井失效。④从工况1～4的计算结果可以看出,尾矿坝设置排水棱体的排水效果较差,只有最下层排水棱体起到排水作用,因此建议去掉上面2层排水棱体。

图6 工况2顺河向剖面水头等值线(单位:m)

图7 工况4顺河向剖面水头等值线(单位:m)

4 结 论

a.虹吸井排渗系统设计原理简单,可操作性强,排渗效果良好,国内尾矿坝内使用较为广泛,采用子结构技术模拟虹吸井的渗流行为效果良好,算法理论严密可靠。

b.在改进截止负压法中引入虹吸井子结构算法,较好地解决了虹吸井作用下渗流场的求解问题,进而可研究尾矿坝密集虹吸井排水降压效果。

c.工程算例的计算结果表明,采用虹吸井子结构法可以很好地模拟尾矿坝在虹吸井、导渗盲沟等排水措施作用下的渗流场分布,为尾矿坝渗流控制设计提供理论支持和技术指导。

[1]王镭,刘中,张有天.有排水孔幕的渗流场分析[J].水利学报,1992(4):15-20.

[2]朱岳明,陈振雷.改进的排水子结构法求解地下厂房洞室群区的复杂渗流场[J].水利学报,1992(9):79-85.

[3]赵坚,沈振中.尾矿坝复杂排水系统渗流计算方法的改进[J].河海大学学报:自然科学版,1997(3),25(2):110-113.

[4]关锦荷,刘嘉忻,朱玉侠.用排水沟代替排水井列的有限元法分析[J].水利学报,1984(3):10-18.

[5]ZHAN Mei-li,SU Bao-yu.New method of simulating concentrated drain holes in seepage control analysis[J].Journal of Hydrodynamics:SerB,1999(3):27-35.

[6]速宝玉,沈振中,赵坚.用变分不等式理论求解有渗流问题的截止负压法[J].水利学报,1996(3):22-29.

[7]张乾飞,吴中如.有自由面非稳定渗流分析的改进截止负压法[J].岩土工程学报,2005,27(1):49-54.

[8]王凤江.虹吸排渗系统在尾矿坝降水工程中应用[J].中国地质灾害与防治学报,2005,16(1):102-105.

[9]朱岳明,陈建余,龚道勇,等.拱坝坝基渗流场的有限单元法精细求解[J].岩土工程学报,2003,25(3):326-330.

[10]崔皓东,朱岳明,吴世勇.有自由面渗流分析中密集排水孔幕的数值模拟[J].岩土工程学报,2008,30(1):440-445.

[11]陈建余.有密集排水孔幕的饱和-非饱和渗流场全精细数值模拟分析[J].中国农村水利水电,2004(11):53-60.