0.   引言

围堰是指为建造永久性水利水电工程设施,修建的临时性围护结构。其作用是防止水进入建筑物的修建位置,以便在围堰围成的基坑内排水、开挖和修筑大坝等工程设施。深厚覆盖层上的土石围堰的防渗体系一般由堰体下部和覆盖层中的混凝土防渗墙、堰身土质心墙或防渗膜、覆盖层底部与岸坡风化基岩中的水泥灌浆防渗帷幕构成^[1-3]。研究围堰及其基础的渗流、应力变形和稳定性,防止围堰发生危及防渗体系安全的过大变形、堰基渗透破坏和滑坡失稳^[1-3],是围堰设计的重要内容。

拉哇水电站上游围堰高60 m,坝址覆盖层总厚度近70 m,其中低渗透性的堰塞湖沉积层厚度约50 m。堰体在一个枯水季节填筑形成。低渗透软土地基上填筑,面临孔压过高,变形过大,抗力过小的难题^[4]。拉哇堰基中低渗透土层的固结系数较低,围堰填筑过程中堰基内产生的超孔隙水压力的消散速度缓慢,需要在覆盖层中设置碎石桩加速堰基排水固结,提高抗滑稳定性和减小变形。围堰填筑和运用过程中堰基中的孔隙水压力、应力变形及稳定安全性的评估是围堰合理设计的关键。高52 m的务坪水库心墙堆石坝的心墙上游侧坝基存在最大深度32 m的软弱土层,其地基处理的研究进行了离心机模拟试验、渗流变形耦合计算等研究工作,为软基上修建大坝和围堰积累了经验^[4-5]。拉哇上游围堰的条件更加复杂,需要考虑包括水中填土、碎石桩和防渗墙施工及堰基开挖全过程的仿真,对渗流与变形耦合计算提出了新挑战。

基于拉哇围堰工程的设计需要,本文发展了深厚覆盖层上的土石围堰填筑与基坑开挖全过程中的渗流与应力变形强耦合动态仿真方法,并在本单位自主研发的岩土工程分析软件系统LinkFEA中实施。本文所发展的方法、对拉哇围堰分析的成果,对低渗透饱和地基上的围堰、码头堆场、排土场等的分析和设计有一定的参考意义。

1.   拉哇上游围堰工程概况

拉哇水电站是金沙江上游河段13级开发方案中的第8级。枢纽主要由混凝土面板堆石坝和位于右岸的溢洪洞、泄洪放空洞、地下厂房等建筑物组成。施工平面布置如图1所示,采用土石围堰拦断河床的隧洞导流方式。上游围堰最大高度约60.0 m,大坝基坑开挖坡高约70 m,下游围堰最大高度约24.0 m。导流建筑物属于临时性建筑物,其级别为3级。

图  1  施工平面布置图

Figure  1.  Construction plan of Lawa Hydropower Station

下载: 全尺寸图片 幻灯片

围堰两岸弱风化基岩裸露,岩性为绿泥角闪片岩（P_txn^a-1）。左岸地形坡度约60°,右岸地形坡度35°～45°。两岸强卸荷带埋深10～15 m,弱卸荷带埋深40～45 m。河床部位基岩岩性为绿泥角闪片岩（P_txn^a-1）,弱风化下限铅直埋深55～75 m,其岩体厚度5～15 m。

上游围堰结构如图2所示。堰基及基坑边坡区域河床覆盖层厚度65～68 m,从上至下分别为表层砂卵石层（Q^al-5）、中部堰塞湖沉积层（Q^l-3层、Q^l-2层）、底部砂卵石层（Q^al-1）,各层详细情况如下：

图  2  上游围堰结构剖面图（基坑开挖后）

Figure  2.  Structural profile of upstream cofferdam

下载: 全尺寸图片 幻灯片

（1）表层Q^al-5砂卵石层,为河床冲积砂卵石层夹少量漂石,在堰基区域其厚度为1.4～4.6 m。

（2）中部Q^l-3堰塞湖相沉积层,含淤泥质粉砂、黏土质砂,在堰基区域其厚度为14.7～18.1 m、分布高程约2514.00～2521.00 m。

（3）中部Q^l-2堰塞湖相沉积层,厚度约为31.4 m,自上而下可分为Q^l-2-③、Q^l-2-②、Q^l-2-①3个亚层,其中,Q^l-2-③层以低液限黏土为主,多呈流塑状,厚度4～8.5 m；Q^l-2-②层以低液限粉土和砂质低液限粉土为主,多呈可塑—软塑状,厚度10～15 m；Q^l-2-①层以低液限黏土为主,局部为低液限粉土,多呈可塑—软塑状,厚度15.2 m。

（4）上游围堰下游堰脚处河床底部靠右岸发育有Ⅰ号透镜体,物质组成为崩石、块石,顺河向长度为325 m、宽度约70 m,最大厚度为32.1 m,分布高程为2489.50～2487.00 m。

（5）河床底部Q^al-1砂卵石层,为卵石、块石夹砂,在堰基区域钻孔揭露的厚度4.6～5 m。

采用LinkFEA软件,对围堰设计初步方案的三维渗流分析和二维渗流与应力变形全过程耦合研究后,确定的优化设计方案如图2所示。堰基混凝土防渗墙位于堰轴线上游147.5 m,墙厚1 m、与防渗墙下部及两岸趾板下的灌浆帷幕和堰表面的防渗膜构成上游围堰的防渗体系。防渗墙下游侧堰基梅花型布置直径1.2 m碎石桩（有效直径1.0 m）。防渗墙下游177 m范围内碎石桩的桩中心线间距和排距均为3 m；防渗墙下游177 m至350 m（下游堰脚）范围的桩间排距均为2.5 m。在围堰下游坡脚处布置两排共5口降水管井以降低基坑开挖过程中底部Q^al-1层中的水压力,防止基坑底部堰塞湖沉积层因Q^al-1层中的水压力过大而隆起失稳。基坑开挖时可以将井中水位降至最下层覆盖层Q^al-1层顶部（高程约2486 m）。优化后的设计方案是否满足防渗体系安全性和堰基抗滑稳定安全性,也需要进一步的渗流与应力变形耦合分析来论证。

2.   耦合计算理论与模拟技术

渗流与变形耦合模拟是在围堰防渗方案论证后进行的,本文重点介绍渗流与应力变形耦合分析的研究成果。

2.1   渗流与变形耦合的基本方程

围堰的渗流存在非饱和区,土的饱和度变化影响渗透系数和含水量。固结试验结果表明,堰塞湖沉积层的渗透系数在0～3.2 MPa的固结压力作用下,变化可达1.6～2个数量级。因此,饱和渗透系数的变化要考虑孔隙压缩的影响。参考非饱和渗流的Richards方程^[6],假定孔隙气压力为0时,容易推导出以指标符号系统表示的变饱和度耦合变形的渗流微分方程：

式中 s为土体的孔隙水饱和度,$k_{\text{r}}(s)$为相对渗透系数（非饱和渗透系数与饱和渗透系数的比值）,是饱和度s的函数；$K_{ij}\text{(}\varphi \text{)}$为渗透张量,其主值随着土体的孔隙率而变化；$p_{\text{w}}$为水压力；${\gamma }_{\text{w}}$为水的重度；z为正方向与重力方向相反的坐标轴,可以理解为位置水头；$\varphi$为孔隙率；$s^{\prime }(p_{\text{w}})$是饱和度对孔隙水压力的偏导数,$s^{\prime }(p_{\text{w}})=\partial (s)/\partial p_{\text{w}}$；t为时间；u为位移向量；i,j为下标,表示坐标轴,i,j=1,2,3；下标中逗号‘,’表示求偏导数,重复下标表示求和。

方程（1）左边的第1项是渗流速度向量的散度,中括号内的表达式是广义达西定律表示的渗流速度；第2项和第3项分别表示因饱和度和孔隙率变化引起的孔隙水体积含量变化率（体积应变与孔隙率变化量相等）。方程中除了渗流基本变量$p_{\text{w}}$和变形基本变量u以外,还包含饱和度s、相对渗透系数$k_{\text{r}}(s)$和渗透张量函数$K_{ij}\text{(}\varphi \text{)}$这3个未知量,需要补充3个关系式。饱和度s由饱和度-吸力（即非饱和土中的气相与水相的压力差）关系函数$s(p_{\text{w}})$计算；相对渗透系数$k_{\text{r}}$由饱和度与相对渗透系数关系函数$k_{\text{r}}(s)$计算。相对渗透系数常采用Mualem's^[7]公式描述。饱和渗透系数与孔隙率的关系,则由固结试验确定。

变形耦合渗流的微分方程由弹性力学的力平衡方程（应力采用弹性力学的符号系统,以拉应力为正）（2）、有效应力公式（3）、本构关系（4）和几何方程（5）这4个方程组成

式中 ${\sigma }_{ij}$为总应力张量；$f_i$为体积力；${{\sigma }^{\prime }}_{ij}$为有效应力张量；$\mathrm{sgn}(p_{\text{w}})$是孔隙水压力$p_{\text{w}}$的符号函数,表示忽略负的孔隙水压力对变形的影响,即当孔隙水压力为负值时,有效应力取值同总应力；$D_{ijkl}$为弹性矩阵张量,由本构关系计算；${\epsilon }_{kl}$为应变张量。下标i,j,k,l=1,2,3。

以上4个方程合并可得耦合渗流的变形微分方程

标量方程（1）和向量方程（6）,基本变量为孔隙水压力$p_{\text{w}}$和位移向量u_i。采用有限元方法耦合求解。方程（6）中的本构关系采用邓肯E-B非线性弹性模型按中点应力法计算。两个方程联立求解,迭代计算收敛后即可获得一个时间步的有限元节点位移增量和当前时刻的孔隙水压力结果,并依据位移增量计算单元高斯点的应力增量。

2.2   围堰全过程耦合仿真若干问题探讨

除了上节的基本理论以外,还需发展若干技术方法,以实现围堰工程全过程渗流与变形的数值仿真。下面撷取模型概化、渗透性随压密变化等问题进行初步探讨。

（1）碎石桩在二维模型中的概化模拟方法

二维模型能以较低的成本快速了解堰基全生命周期的渗流与变形情况。稳定性分析也主要基于二维有效应力和孔隙水压力模拟结果。因此,二维计算十分重要。

二维模型中,碎石桩区域可以分成桩和基土两种材料相间分布,按照与三维排水速度等效来确定基土的宽度。碎石桩的等效宽度取有效置换率与基土宽度的乘积。

以地基中等间距布置的碎石桩为例,在平面上将碎石桩与基土复合地基概化成如图3所示的一维模型。碎石桩的半径为r₀,桩间距为a,复合地基可以分割成形状完全相同的三角形区域,如三角形ADE。扇形ABC是碎石桩,圆心角为$\text{π}$/6,BDEC区为基土。一维模型中AB段是碎石桩,BC段是基土。

图  3  碎石桩二维模型概化图

Figure  3.  Schematic diagram of 2D model for gravel piles

下载: 全尺寸图片 幻灯片

碎石桩对基土起水平径向排水作用,将三角形ADE区域等效为半径R的扇区。半径为R的园即1个碎石桩的影响范围。R-r₀即基土的最大排水距离。对于间距为a,排距为b的碎石桩,R可按如下公式计算：

圆心角为$\theta$的二维扇形区域CBGF土体自FG向CB的径向排水,区域中任意一点到碎石桩中心A点的距离为r,过该点的等势线长度为$\theta r$。等效为矩形区域时,该点的等效水平渗透系数为

式中,k为原基土渗透系数,整个区域的等效水平排水渗透系数为

区域的排水固结因数与渗透系数成正比,与排水距离的平方成反比,若基土的渗透系数仍取土本身值进行计算,扇形区域CBGF按照固结时间因数等效的矩形区域排水距离$H^{*}$为

二维模型中桩间土的长度等于c=2$H^{*}$,碎石桩的宽度d根据有效置换率计算。

（2）二维模型中包含三维绕渗效应的方法

本项目研究中,先进行三维渗流分析,然后进行河谷中心剖面的二维渗流与应力变形耦合分析。三维渗流研究中发现,围堰填筑完成且上游水位2594.6 m、下游水位2541.0 m的稳定渗流工况,防渗墙上游覆盖层中的水头降落要少于下游覆盖层中的水头降落。二维模型中,如果仅包含模型底部的绕渗,水头分布结果与三维差异很大。因而,二维模型中应增加考虑三维绕渗效应的方法。

首先依据三维模型典型工况的稳定渗流计算,获得防渗剖面处覆盖层底面以上基岩中的绕渗流量与高程的关系,计算出各高程的绕渗流量与上游水位和防渗墙后水位或计算点高程之差的商,得出绕渗系数。在二维模型中防渗剖面后沿着高程取1列节点作为绕渗入渗面,设置1列绕渗单元,绕渗单元的宽度可取1m,渗透系数取绕渗系数与覆盖层平均宽度之商,绕渗单元的上游面节点的边界条件为上游水头。从而在二维模型中,自动实现不同水位条件下防渗平面下游侧的绕渗流量叠加。详细的方法将另文投稿。

二维模型围堰正常高水位2594.6 m、下游水位2541.0 m时稳定渗流工况中是否包含绕渗处理的等水头线对比结果如图4所示。图中可见上游侧和下游侧的等水头线分布均有较大差异。设置绕渗单元,对等水头线结果影响还是比较大的。

图  4  正常高水位时的二维渗流等水头线对比

Figure  4.  Contrast of results of iso-head line at normal high water level in 2D models

下载: 全尺寸图片 幻灯片

设置了绕渗单元后,二维全过程耦合仿真中自动将渗流的三维绕渗效果纳入进来,使模拟结果与实际更接近。

（3）渗透系数的时空变化模拟

一般低渗透土层的渗透系数,随着压密变化很大。从固结试验中,不但可以获知土体的压缩特性、固结系数,还能推知土体的渗透系数。按如下公式可以从固结试验中获得渗透系数与固结压力的关系。

式中,$c_{\text{v}}$为固结系数（m²/s）,k为渗透系数（m/s）,e为孔隙比（无量纲）,${\alpha }_{\text{v}}$为压缩系数（Pa^-1）,${\gamma }_{\text{w}}$为水的重度（N·m^-3）。

按照式（11）和固结试验中的固结压力、压缩系数、固结系数、孔隙比关系,推算出孔隙比与渗透系数之间的对应关系。有限元计算中根据高斯点的孔隙率结果,即可获得各计算时步的渗透系数。

拉哇堰基堰塞湖沉积土层的渗透系数与固结压力的关系如表1所示。表中堰塞湖沉积土层中的4个土类的固结垂直压力范围由0～50 kPa增加到100～200 kPa时,渗透系数降低到（1/4.0～1/12.6）；由0～50 kPa增加到400～800 kPa时,渗透系数降低到（1/7.6～1/31.9）。可知同一土层不同深度处和同一土体在荷载变化过程中渗透系数的变化很大。耦合计算中模拟低渗透土层的渗透系数时空变化是必要的。

表  1  固结试验中土层渗透系数与固结压力关系

Table  1.  Soil permeability coefficient vs. consolidation pressure in consolidation test

土类型	各固结压力下的渗透系数/(m·s-1)
	0~50 kPa	50~100 kPa	100~200 kPa	200~400 kPa	400~800 kPa	800~1600 kPa	1600~3200 kPa
粉土质砂	2.38×10-9	8.83×10-10	4.39×10-10	2.94×10-10	2.03×10-10	1.31×10-10	8.60×10-11
砂质粉土	2.31×10-9	1.26×10-9	5.75×10-10	4.83×10-10	3.06×10-10	1.77×10-10	6.23×10-11
粉土	5.92×10-9	1.28×10-9	4.70×10-10	3.02×10-10	1.86×10-10	9.04×10-10	5.94×10-11
黏土	2.41×10-9	7.20×10-10	5.16×10-10	3.86×10-10	1.77×10-10	1.27×10-10	2.92×10-11

下载: 导出CSV 

| 显示表格

（4）施工过程的仿真

在天然地基中的防渗墙和碎石桩,在全过程耦合仿真中其所在部位原为覆盖土层,是在仿真过程中因施工而变换为新材料,可依据其实际形成的计算级通过置换单元对应的材料编码实现。置换时单元的应力需要给初值,且防渗墙这种刚性结构物与周围土体之间还需要设置接触面单元。建模与网格剖分时需要规划材料替换和接触面设置要求。在材料替换的计算级,实现单元材料和类型的转换。单元应力的初值,可以采用原位替换前的应力,原普通实体单元转换为接触面单元时,需要进行应力张量的坐标变换,从而求出接触面的法向应力和切向应力作为接触面单元的初始应力。两种材料的重度之差,作为自重荷载在替换计算级施加于模型之中。

（5）算法的收敛性问题

渗流和变形的耦合求解,既包含有限元渗流计算和变形计算中存在的算法问题,还包含二者耦合产生的问题。非饱和区的渗透系数和外边界溢出面范围需要在迭代求解过程中确定,围堰渗流场中还存在内部溢出面,其范围和传递流量也需要在迭代过程中确定,是非线性问题,迭代的收敛性问题突出^[8-9]。土的变形求解过程中,变形模量矩阵需要由应力状态确定,也是非线性迭代问题。即使只是饱和多孔弹性介质的Biot固结计算,也还存在收敛性问题^[10-11]。本文的渗流计算方法参见文献[8,9],包含了非饱和渗流参数、内部溢出面和边界溢出面的非线性迭代。全过程按时间轴分成若干计算级,计算级又分成若干计算时步。迭代分成2层：外层是本构模型模量矩阵迭代计算,采用中点应力法确定；内层是渗流计算迭代,该层中单元高斯点的变形模量矩阵不变,渗透系数随着高斯点的孔隙率和饱和度变化。内层中两次计算的节点正孔隙水压力差的最大值作为收敛变量,小于给定值或迭代次数达到设定值完成一个内层迭代。研究过程中发现,荷载增加速率大的计算时步往往在给定的迭代次数中达不到收敛标准,且相邻两次迭代计算的误差并不一定随迭代次数增加而减小。

3.   计算模型与条件

3.1   计算模型与条件

围堰全过程模拟共分为61个计算级进行,每个计算级依据其时间间隔长短又分成1～5个计算步。各计算级信息如表2。第1计算级计算天然地基的初始应力场和渗流场；第2级进行碎石桩施工；第3～9级戗堤施工；第10～17级为Ⅱ、Ⅲ区域的填筑；第19～41级进行Ⅳ区填筑,其间防渗墙在第26级形成,第32级降水井形成,其后降水井水位与下游水位或基坑开挖底面高度保持基本一致,直到基坑开挖至2486 m高程时,降水井中水位保持抽水至2486 m高程不变。第42级,进行土工膜铺设,上游水位上升至2566 m,下游水位不变,历时1 d。第43级上游水位上升至最高水位；第44～60级上游水位维持不变,基坑逐级开挖直至覆盖层底部。

表  2  计算级信息表

Table  2.  Computation-level information

计算级	计算时刻/d	上游水位/m	下游水位/m	填筑高程/m	备注
—	0	2545.74	2545.74		初始条件
1	0.5	2545.74	2545.74		天然地基形成,初始应力场
2	1	2545.74	2545.74		碎石桩形成
3~8	1.5~4	2545.74	2545.74	2540~2550	I区戗堤自堰基面填筑至2550 m
9	47	2547.11	2541.00		无填筑,上下游水位改变
10~17	47.5~77	2547.11	2541.00	2553	II、III区自堰基填筑至防渗墙施工平台
18	78	2551.18	2541.00		无填筑,上游水位上升
19~25	84.18~121.27	2551.18	2541.00	2555~2567	IV区碾压石渣自2553 m填筑到2567 m
26	127.45	2551.18	2541.00	2569	防渗墙材料替换,绕渗单元添加
27~41	133.64	2551.18	2541.00	2571~2597	IV区填筑到堰顶,第32级降水井形成
42~3	215.00	2566.00	2541.00		复合土工膜铺设,上游水位上升
43	225.00	2594.60	2541.00		上游水位上升至2594.6 m
44~60	336.60~733.86	2594.60	—		覆盖层逐层开挖,至基岩面
61	1000.00	2594.60	—		开挖后挡水期

下载: 导出CSV 

| 显示表格

3.2   有限元网格

桩间距3 m区域的桩间土等效宽度3.14 m,桩宽0.30 m；水平方向桩间土剖分成15个单元,桩剖分1个单元。桩间距2.5 m区域的桩间土宽度2.31 m,宽度0.33 m；水平方向桩间土剖分11个单元,桩剖分1个单元。防渗墙（厚1 m）水平方向剖分3个单元,墙两侧设置0.1 m厚接触面单元（置换前为普通单元）。防渗膜厚取为0.1 m,其后侧设置0.1 m宽的接触面单元。防渗墙和防渗膜下游面设为可能的渗流内部溢出面,其后部接触面单元同时作为渗流接触面单元。防渗墙和防渗膜后设置宽度1 m的绕渗单元,其下游侧边与接触面单元的下游侧边重合,其上游侧边列入渗流的上游水位边界（绕渗单元不做变形计算）。降水管井的模拟是在第一排管井处（堰下游坡脚外8.17 m）覆盖层底部沿深度方向设置水头边界。填筑和开挖通过有限元中通用的“生”、“死”单元来实现。开挖边界面单元的应力边界条件按边界面上正应力和剪应力为0的条件在边界面施加反向的分布力。二维整体有限元网格如图5所示,模型共有58320个节点,59102个单元。高程作为垂直坐标,水平坐标原点位于围堰顶部的中轴线上,方向指向围堰下游。

图  5  模型有限元网格图

Figure  5.  Finite element mesh

下载: 全尺寸图片 幻灯片

3.3   水压力与渗流边界条件及模拟方法

以往水中土体按浮重度计算,不计算边界上的水压力。全过程耦合过程中,水中填筑体外边界也因填筑而变化,水中填筑体内部也有超孔隙水压力产生,需要新方法来模拟。如第3～17计算级中,堰身戗堤和Ⅱ区、Ⅲ区填筑,其中部分单元是水下填筑形成,有些水下外部边界消失,也产生新的水下外部边界。本文的方法是对于水中填土消失的边界,在当前计算级的第一个时间步,在消失的水中原边界上施加反向的水压力,而水中填筑产生的新边界,水压力荷载增量是边界上的水位作用下的全部水压力（老边界上只计算水位变动引起的增量水压力）。填筑体的重度,则取饱和重度,初始孔隙水压力取0。实现了水中填筑过程土体的渗流边界和外力作用的直接模拟。

3.4   材料参数

计算用到的土的参数包括干密度、孔隙率、填筑饱和度、莫尔库仑强度指标、本构模型参数、渗透系数、非饱和相对渗透系数与饱和度关系、饱和度与吸力关系曲线。土体的渗透破坏判断用到允许渗透坡降值。第1级天然覆盖层应力计算时覆盖层初始应力的修正用到天然覆盖土层的侧压力系数（依据计算所得的垂直正应力和侧压力系数修正水平向正应力）。低渗透土层考虑渗透系数与固结压力的关系,用到固结试验中的水平固结系数和垂直固结系数（考虑水平与垂直向渗透系数的差异）,压力–压缩模量关系。因篇幅限制,本文仅列出主要参数。

土体的本构模型采用邓肯E-B模型,其密度、孔隙率、强度指标与邓肯E-B模型参数列于表3。防渗墙、防渗膜和基岩采用线弹性模型。防渗墙的弹性模量取1500 MPa,泊松比取0.2。防渗膜的弹性模量取100 kPa,泊松比取0.49。各风化程度基岩的渗透系数和弹性参数因篇幅限制从略。接触面单元应力–应变关系采用邓肯-克拉夫模型,法向模量在受压时取10 GPa,受拉时取1 kPa。防渗墙与两侧土体之间、防渗膜与堰体之间在防渗墙和防渗膜生成后设接触面单元,其强度参数取接触土体的参数值,模型中所有的接触面的$k_1$,n和$R_{\text{f}}$均取值100,0.57,0.68。

表  3  填筑料和覆盖层土体邓肯E–B模型参数

Table  3.  Parameters of soils for Duncan E–B model

材料	干密度/(g·cm-3)	孔隙率	c/kPa	$\phi$/(°)	K	n	Rf	Kb	m	Kur
抛填石渣	1.90	0.30	0	38	900	0.25	0.85	393	0.22	1500
抛填砂砾料	1.60	0.38	0	29	1000	0.28	0.75	400	0.22	1200
碾压石渣	2.05	0.25	0	21	900	0.25	0.85	393	0.22	1500
碎石桩	2.05	0.25	0	38	900	0.25	0.85	393	0.22	1500
Qal-5	2.05	0.25	0	35	1000	0.35	0.8	340	0.20	1200
Ql-3	1.40	0.48	28.7	22	125	0.57	0.68	90	0.56	150
Ql-2-③	1.36	0.50	45.0	20	87	0.58	0.62	60	0.58	105
Ql-2-②	1.38	0.49	31.0	21	100	0.56	0.65	73	0.56	120
Ql-2-①	1.36	0.50	42.0	20	85	0.57	0.63	60	0.57	102
Qal-1	1.95	0.25	10.0	36	1000	0.35	0.80	340	0.20	1200

下载: 导出CSV 

| 显示表格

河床覆盖层的渗透系数和固结系数列于表4,其中固结系数取定值。低渗透土层的压缩模量与固结压力的关系列于图6,程序中依据这一关系和土层的侧压力系数,折算出固结试验中体应力与体积模量的关系。再依据渗透系数与体积压力、固结系数的关系和高斯点的有效应力状态,计算低渗透覆盖土层单元中各高斯点的渗透系数,从而模拟渗透系数在空间和时间上的变化。复合土工膜实际膜厚0.75 mm,材料渗透系数取5×10^-11 cm/s,建模时土工膜厚度0.1 m,计算时等效渗透系数6.67×10^-9 cm/s。抛填石渣渗透系数5.00×10^-1 cm/s,砂砾料和碾压石渣渗透系数5.00×10^-2 cm/s,其他材料渗透系数见表5。

表  4  河床覆盖层渗透与固结参数

Table  4.  Permeability and consolidation parameters of soils

土层	侧压力系数K0	渗透系数/(cm·s-1)	垂直固结系数/(cm2·s-1)	水平固结系数/(cm2·s-1)
Qal-5	0.35	5.5×10-1	—	—
Qal-1	0.36	3.0×10-2	—	—
Ql-3	0.50	6.8×10-5	4.4×10-3	4.5×10-3
Ql-2-③	0.60	2.9×10-6	3.1×10-3	3.9×10-3
Ql-2-②	0.55	3.5×10-6	3.8×10-3	4.2×10-3
Ql-2-①	0.66	2.0×10-6	3.1×10-3	3.9×10-3

下载: 导出CSV 

| 显示表格

图  6  侧限固结试验压缩模量与固结压力关系

Figure  6.  Compression modulus vs. consolidation pressure in consolidation tests

下载: 全尺寸图片 幻灯片

表  5  覆盖层以外的材料渗透系数

Table  5.  Permeability coefficient of materials beyond overburden

材料	碎石桩	防渗墙	帷幕	抛填石渣
渗透系数/(cm·s-1)	5.00×10-2	1.00×10-7	1.00×10-5	5.00×10-1

下载: 导出CSV 

| 显示表格

4.   计算结果与讨论

获得了61个计算级的节点位移、孔隙水压力和高斯点应力结果。从这些结果中,可以了解到围堰填筑、挡水和基坑开挖全过程中的孔隙水压力、变形和应力过程,并可以根据各级的应力结果,计算抗滑稳定安全系数。因篇幅限制,仅介绍戗堤填筑形成和填筑完成2个典型计算级的部分渗流、位移和应力结果,并分析渗流应力的耦合作用和碎石桩的排水效果。

4.1   戗堤填筑形成时的计算结果

戗堤填筑前,地基中已经施工了碎石桩,构成了复合地基。

（1）地基中的孔隙水压力

第8计算级戗堤自建基面3 d均匀上升到2550 m高程。不同高程的水平位置–水头关系如图7所示。Q^l-2-1土层的底部（结合图9看土层的位置）,即覆盖层底部强透水的Q^al-1层顶部,水头随水平位置略有变化,最大值2546.42比静水位2545.74 m大0.68 m,且与图中范围的最小值2546.02 m（位于x=50 m处）水头相差仅0.40 m。而Q^l-3土层顶部,即强透水的Q^al-5层底部,戗堤中部的水头显著大于两侧,最大值2548.22 m,超过静水位2.48 m。此处水头大是由于下部低渗透土层中的孔隙水因孔隙压缩向上排出（主要通过碎石桩）后,在强透水层中向两侧流动,而非Q^al-5层孔隙压缩产生超孔隙水压力引起。Q^l-3土层的中部最大水头2554.51 m,超过静水位8.73 m,超孔隙水压与填筑土的垂直有效自重荷载（堰基表面上砂砾石柱的垂直有效土重170.51 kPa,垂直荷载速率57 kPa/d）之比为50.2%。戗堤下覆盖层中虽然布置了3 m间距1.2 m直径碎石桩,低渗透的堰塞湖沉积粉土层中仍然超孔隙水压力消散不足一半。

图  7  戗堤填筑完成时不同部位水平位置–水头关系

Figure  7.  Horizontal position and water head at different parts at completion of berm filling

下载: 全尺寸图片 幻灯片

图  8  戗堤中部碎石桩与桩间土的水头高程关系

Figure  8.  Water head vs. elevation in gravel piles and in soils at middle part of two piles

下载: 全尺寸图片 幻灯片

图  9  Q^l-3土层中部水平位置与桩和桩间土的正应力关系

Figure  9.  Horizontal position vs. normal stresses in piles and soils in Q^l-3

下载: 全尺寸图片 幻灯片

戗堤中部位于水平坐标-76.59 m的碎石桩边线和-74.95 m处的桩间土中线高程–水头关系如图8所示。碎石桩在Q^l-3土层顶部水头比Q^l-2-1土层底部高1.79 m,说明碎石桩向底部强透水层Q^al-1排水。桩间土中线的高程–水头关系的变化很有意思。堰塞湖沉积层z轴方向网格第1个内部点和最底一个内部点的值,远远大于其他内部点的值,不符合土层接近上下强透水边界超孔隙水压力消散更多的规律。这不是物理现象本身,而是数值模拟中产生的问题。在对饱和多孔线弹性介质的一维固结计算的测试算例中,四周水平位移约束底部水平垂直位移约束,底部和四周均不透水,顶部排水,顶部施加垂直荷载。当模型仅为1个1次单元,且加载持续时间与渗透系数的乘积足够小时,底部节点的超孔隙水压力接近垂直荷载的2倍。出现这个结果的原因是,饱和土单元在快速荷载作用下因为来不及排水其孔隙压缩量接近于0,而顶部节点的孔隙水压力因为是排水边界而为固定值,在整个单元中的超孔隙水压力平均值要达到接近于垂直荷载这个条件,不可避免地就计算得出了底部节点的超孔隙水压力接近荷载的2倍这个结果,比实际的物理结果大100%。由此可见,耦合计算中近强透水边界处的超孔隙水压力可能存在较大误差。因为水平向桩间土剖分单元较多,网格尺度仅为0.2 m,仅为垂直向的网格尺度2.3 m的10%,因此,水平方向的水头连线比较平顺,数值计算本身造成的误差较小。近强透水层的土层边界单元越密,误差越小。因此,对于快速荷载来说,靠近排水边界的第一个内部点的结果,尤其是其值为最大值的时候,可能大大高于真实物理情况,不能将此节点的结果,作为超孔隙水压力的特征值。若要减少数值模拟本身的这种误差,低渗透性土层的网格在与高渗透性土层接触处应取较小的尺度。Q^l-3土层内部的第3个点水头值2554.51 m,作为该水平位置高程方向上的最大值。

（2）地基中的有效应力

Q^l-3土层中部土单元或桩单元高程约2528.3 m的1排高斯点的正应力与水平坐标关系连线如图9所示。土中的垂直向正应力${\sigma }_z$和水平正应力${\sigma }_x$水平向波形变化,在桩边处于波峰,桩间土中部处于波谷。土中的正应力${\sigma }_z$和${\sigma }_x$在两桩边线范围的变化幅度基本上等于孔隙水压力的变化幅度。桩与土的应力比较,在戗堤的堤脚内桩的垂直正应力${\sigma }_z$比土中大,在堤脚外比土中小；而水平正应力${\sigma }_x$桩内与其临近的土中的值是比较接近的。垂直向应力${\sigma }_z$在戗堤中线最大,向堤脚两侧减小,而水平正应力${\sigma }_x$则在戗堤中线部位出现极小值,堤脚处出现极大值。究其原因,是垂直向的总应力（指弹性力学的应力,相对有效应力而言）基本与垂直荷载相等,水平与垂直方向正应力比值${\sigma }_x$/ ${\sigma }_z$看成是侧压力系数,这个系数是小于1.0的。由于超孔隙水压力的影响,低渗透土层填筑区域的垂直向正应力增加,水平向正应力减小。Q^l-3土层中部水平坐标与土的剪应力水平关系如图10。戗堤中心的剪应力水平在桩间土中部区域已经达到了1.0,即剪应力已经达到甚至超过了抗剪强度（未进行应力迁移计算）,而桩间土边缘的应力水平较低,可见超孔隙水压力对应力水平的影响是很大的。低渗透土层上快速填筑荷载作用下戗堤中部区域桩间土中部土体的有效应力路径是水平正应力减小,垂直正应力增加,剪应力水平急剧增加；戗堤中部区域桩侧土则是垂直正应力与水平正应力同步增加,但垂直正应力增加幅度较大,剪应力水平也增加。堰基桩间土的应力路径沿水平向变化急剧。随着后续荷载的施加和超孔隙水压力的消散,应力路径变化极其复杂。

图  10  Q^l-3土层中部水平位置与土的剪应力水平关系

Figure  10.  Horizontal position vs. shear stress level of soils in Q^l-3

下载: 全尺寸图片 幻灯片

桩对软土地基排水固结以外,还有桩体材料力学性能改善所带来的加强作用。从图9可以看出,在x∈（-104.1,-48.5）的荷载作用区域,桩中的垂直正应力比土大,能发挥桩的垂直向承载加强作用,也能提高水平向的摩擦抗剪强度,而荷载作用区域以外,桩的垂直正应力比桩间土小。由此可知,其摩擦抗剪强度,在荷载区域以外是很难发挥出来的。一般在边坡加固中,常采用在坡脚处设置阻滑桩,以起到“固脚”的作用。软土地基上在填筑体坡脚外打设碎石桩,排水固结作用显著,而碎石桩抗剪加强作用有限。基于此新认识,本工程围堰坡脚以外和基坑开挖边坡上原来拟设置的碎石桩,在优化方案中全部取消。

（3）堰体和地基中的位移

戗堤填筑完成时位移如图11所示,最大向上游水平位移0.26 m,向下游水平位移0.24 m；最大沉降量0.71 m,位于覆盖层表面。最大隆起0.07 m,位于x坐标-154.4 m处的覆盖层表面。堰基位移的态势是向下部和两侧挤压。没有碎石桩的天然地基的计算结果,则填筑体外侧向上隆起很严重,水平位移也更大。对于应力路径在时空上如此复杂变化,且实际上大量局部达到抗剪强度的情况,其位移模拟结果要达到定量的程度是很困难的。本文采用邓肯E-B这种比较简单的非线性弹性模型,且不考虑土体的应力超过其抗剪强度以后的应力迁移计算,其位移结果本文认为也是定性的程度大于定量。若采用如修正剑桥模型等复杂的弹塑性本构模型,在剪应力水平达到1.0时的塑性迭代对于本文的这个问题,计算收敛是很难的。因此,在现有的模拟水平下,拉哇上游围堰这种复杂的情况,位移模拟结果是定性多于定量。

图  11  戗堤填筑完成时的位移

Figure  11.  Displacement at completion of berm filling

下载: 全尺寸图片 幻灯片

4.2   围堰填筑到顶时的结果

堰体填筑213 d到顶（第41计算级）,此时上游水位2551.18 m,下游水位和底部覆盖层减压井水位2527.66 m（基坑中已抽水）,上下游水位差23.52 m。

第41级围堰中的等孔隙水压力和等水头线如图12所示,其中图12（a）为孔隙水压力等值线,图12（b）为水头等值线。图13为第41级堰塞湖沉积层中水平坐标–水头关系。图12,13结合来看,堰基碎石桩区域的孔隙水压力和水头在桩间土和桩中的水平方向是波动变化的。堰基下覆盖层的水头高于上游水位,表明堰基中仍然有比较大的超孔隙水压力。图12（a）还可见堰基x=120 m处的孔隙水压力,显著受到了抽水井的影响。图12（a）中有两条0孔隙水压力线存在,查看局部等值线图,证实是绘图软件插值误差造成。

图  12  围堰填筑到顶时的等水头线和孔隙水压力

Figure  12.  Pore water head and pressure when the cofferdam is filled to the top

下载: 全尺寸图片 幻灯片

图  13  围堰填筑完成时堰塞湖沉积层中水平坐标–水头关系

Figure  13.  Horizontal coordinate vs. water head in sedimentary layer at cofferdam filling to top

下载: 全尺寸图片 幻灯片

图14为碎石桩边线位于x=-76.9 m和x=0.9 m两处桩边线与桩右侧桩间土中线的水头对比。前者位于堰上游坡的中部,后者靠近堰轴线。碎石桩的顶部水头高于底部,可知底部Qal-1层起到了良好排水作用。图12,13与图14结合,可以获得堰塞湖沉积层的水头的整体分布情况。Q^l-3层的顶部,即强透水的Q^l-5层的底面,在防渗墙的上游侧,水头与上游水位相等,防渗墙以右,在x=-50 m左右,水头达到峰值2564.5 m,高于上游水位,比下游坡脚强透水Q^l-5层的底面2534.1 m高30.4 m,说明Q^l-5层中沿着水平方向还是存在较大的水头梯度,也就是说Q^l-5层及其上部堰体填筑料的渗透系数,对堰基中的孔隙水压力还是存在较大的影响,如果实际渗透系数比计算渗透系数大,则堰基上部强透水层中的水头梯度会降低,堰基中的水头总体上也会降低,反之,则还会提高。Q^l-3层的中部最大水头2595.6 m,位于该层节点x=2.53 m和x=2.83 m的节点上,比右侧与碎石桩接触的节点的水头2561.3 m（左侧节点水头2561.8 m）高出34.3 m；近防渗墙处最左1排碎石桩的右侧桩间土中的水头,比碎石桩中仅高1.8 m左右,一方面是其所处位置在填筑高程超过防渗墙平台以后,堰体填筑的后续附加垂直正应力较小,另一方面其较早时间的填筑附加应力作用产生的超孔隙水压经历了较长时间的排水固结而消散程度较大。间距2.5 m区域桩间土与桩中的水头差,即波动曲线的波峰与波谷差值,要大大低于碎石桩3 m桩距区域。在2个区域的分界线处,3 m间距区域侧的水头差为32.0 m,而2.5 m间距区域的水头差是20.4 m,可见超孔隙水压力的消散程度对碎石桩的间距是十分敏感的。

图  14  填筑完成时沉积层中碎石桩与桩间土中线水头对比

Figure  14.  Water head in gravel piles and in soils at cofferdam filling to top

下载: 全尺寸图片 幻灯片

5.   结论

针对拉哇水电站上游围堰的设计方案,实现了包含碎石桩和防渗墙施工、土体水中填筑、水位动态变化和基坑开挖全过程的渗流与变形耦合仿真,包含了低渗透土层渗透系数随孔隙压密变化的模拟。二维模型中还包含了三维绕渗影响的模拟。直接服务于拉哇水电站的围堰方案设计论证与优化。

对拉哇的工程分析表明：

（1）深厚低渗透覆盖层在围堰填筑过程中存在较大的超孔隙水压力,其累计与消散对堰基碎石桩的间距十分敏感。

（2）堰基中的水头、垂直和水平有效正应力、剪应力水平在碎石桩处理区沿水平方向是波动变化的,孔隙水头在碎石桩中处于波谷、在桩间土中部处于波峰。

（3）填筑体自重作用下,堰基中超孔隙水压力的产生,使垂直向正应力增加,水平向正应力减小,堰基桩间土局部的剪应力水平快速提高,堰基中容易产生较大的沉降和水平位移,对围堰的变形和稳定性均不利。

（4）碎石桩在填筑体下部同时起到加速排水固结和桩体加强作用,在填筑体坡外加强作用有限。