基于GeoStudio的石渣坝渗流场与应力场耦合分析

基于GeoStudio的石渣坝渗流场与应力场耦合分析
作者：李晔艾丛芳金生丁伟业
来源：《南水北调与水利科技》2015年第08期
摘要：选取香水水利枢纽工程的黏土心墙石渣坝的横剖面作为研究对象，利用GeoStudio 的耦合分析功能对其进行非饱和渗流场-应力场的耦合分析及边坡抗滑安全系数的计算，得到相应的非饱和渗流场分布、水平和竖向位移分布、抗滑安全系数值。

通过与非耦合情况下的应力场及渗流场对比，表明边坡位移在深层区变化缓慢，在浅层坡面位移变化比较大。

与非耦合情况下的边坡抗滑安全系数值对比，表明上游坡脚稳定性不够，应采用放缓边坡以及铺石压坡等方法提高边坡稳定性。

进而说明耦合分析更具有实际意义。

关键词：石渣坝；渗流场；应力场；耦合分析；抗滑安全系数
中图分类号：TV314文献标志码：A文章编号：1672-1683（2015）-002-0001-04
石渣坝受到荷重作用，会改变石渣坝体内的渗流场分布[1]，这是由于荷重作用降低了坝体内相应的孔隙度，从而引起渗透系数的改变。

渗流场分布的改变会直接导致渗流压力及渗流体积力也随之变化[15]，这两个参数的变动进而会影响到石渣坝体的外荷载，最终导致石渣坝体的应力场分布发生变化。

可见，互相作用以及互相影响的渗流场与应力场间存在一种耦合关系[15]，一种使得应力场分布和渗流场分布都会发生变化的关系。

正是基于此原因，当前工程建设中愈发重视应力场与渗流场的耦合分析。

GeoStudio作为全球最知名的岩土工程分析软件之一，最常用的三个基本模块分别是SEEP/W、SIGMA/W、SLOPE/W，所有模块可以在同一环境下运行。

SIGMA/W与SEEP/W耦合分析时，首先要通过SEEP/W得到瞬时孔隙水压力，再将结果导入SIGMA/W中，从而得出在渗流场作用下不同时刻的应力应变。

本文基于GeoStudio软件的耦合分析功能[15]，即SIGMA/W模块与SEEP/W模块以及SLOPE/W模块的耦合功能，进行石渣坝渗流场与应力场耦合分析，同时分析出石渣坝边坡的抗滑移稳定性。

[BT2][STHZ]1工程概况
[BT3][STHZ]1.1地质地貌
香水水利枢纽工程位于牡丹江右岸支流大石河的中上游[6]，大石河流域属山区性河流，河源位置在敦化市翰章乡富尔岭东麓，流域内最高山顶高程为1 227 m。

坝址地处于富尔岭和新开岭余脉的低山丘陵台地区[6]，地形连绵起伏，多山地，绝对高程700～1 000 m。

[BT3][STHZ]1.2地层岩性
香水水利枢纽工程勘察区地层划属松花江地层区，吉林延边分区，延边小区区内出露的地层有[6]：中生界侏罗系中统帽儿山组（J2m）；白垩系上统龙井组（K2l）；新生界第三系中新统土门子组（N1t）和上新统船底山玄武岩（βN2C）；第四系下更新统白金玄武岩（βQ1b）和上更新统二级阶地及冰川、冰水堆积物（Q3）、马连河玄武岩（βQ3m）；第四系全新统地层（Q4）。

[BT3][STHZ]1.3地质构造
本区内华夏式构造体系为主要的构造体系，自新生代时期开始活动，于第三纪时期形成了以西清真至高松树一带为沉降中心的沉降盆地[6]，本测区位于其西南边缘地带。

另外，本区北东向压扭性断裂较发育，主要有敦化-镜泊湖断裂、色洛河断裂。

新生代以来[6]，本区受喜山运动的影响，玄武岩浆沿着断裂带频繁喷溢，形成了以玄武岩为特征的地层结构和受玄武岩控制的地貌轮廓。

本区地震基本烈度在《中国地震动参数区划图》（GB 18306-2001）上标定为小于Ⅵ度区，对应的地震动峰值加速度小于005 g。

[BT3][STHZ]1.4工程概述
香水水利枢纽工程规模为中型水库，工程等别为Ⅲ等[6]，水库正常蓄水位为69410 m，相应库容为3 993万m3；校核洪水位为69543 m，总库容为4 679万m3；香水水利枢纽由黏土心墙石渣坝、岸边开敞式溢洪道、坝下输水涵洞、供水管道等组成（图1）。

大坝坝型为黏土心墙石渣坝，最大坝高421 m，最大坝长1 3376 m。

图1香水大坝示意图
[BT2][STHZ]2计算原理
[BT3][STHZ]2.1非饱和土Darcy定律
非饱和区Darcy′slaw的应用近似于饱和区 [7]，只不过此时的渗透系数k不再是常量，而是一种关于水的容量以及孔隙水压力的函数，也就是所谓的渗透系数函数，它可以准确表示非饱和区的渗透系数随基质吸力的变化。

q=-k（θ）[FK（W*9。

1*3] [XZ（180][JX1*3]Δ[XZ）][JX-*9/9][FK）] [KG-1*7]H[JY]（1）
式中：q为流量；k（θ）为渗透性函数；[FK（W*9。

1*3] [XZ（180][JX1*3]Δ[XZ）][JX-*9/9][FK）] [KG-1*7]H为总水头梯度。

[BT3][STHZ]2.2渗流控制方程
体积水容量的变化与孔隙水压力的变化可以用下式表达[11]：
θ=mwuw[JY]（2）
式中：uw为孔隙水压力；mw为土水关系函数斜率。

总水头定义如下[10]：
H=[SX（]uw[]γw[SX）]+y[JY]（3）
式中：y为高度；γw为水的容重。

[JP+2]根据非饱和土Darcy定律可以建立瞬态渗流方程[710]：[JP]
[SX（][]x[SX）][JB（（]k（θ）[SX（]H[]x[SX）][JB））]+
[SX（][]y[SX）][JB（（]k（θ）[SX（]H[]y[SX）][JB））]+Q=mwγw[SX（]H[]t[SX）][JY]（4）
式中：H为总水头；k（θ）为渗透性函数；Q为边界流量；θ为体积水容量；t为时间。

令λ=mwγw，根据虚位移原理对式（4）进行有限元离散得[712]：
t∫[WTB1X]v[WTBX][B]T[C][B]d[WTB1X]v[WTBX]{H}+
∫[WTB1X]v[WTBX]λ〈N〉T〈N〉d[WTB1X]v[WTBX]{H}，t=
q∫A〈N〉TdA[JY]（5）
[JP+1]式中：[B]为梯度矩阵；[C]为单元渗透系数矩阵；{H}为节点水头矢量；{H}，t为水头随时间的变化；〈N〉为形函数矢量。

[JP]
二维渗流分析模型中单元节点的厚度在整个网格中一般被认为是常数，因此式（5）可改写为[712]：
t∫A[B]T[C][B]dA{H}+t
∫Aλ〈N〉T〈N〉dA{H}，t=
qt∫L〈N〉TdL[JY]（6）
式中：t为单元厚度。

将式（6）简化为：
[K]{H}+[M]{H}，t={Q}[JY]（7）
式中：单元特征矩阵[712]
[K]=t∫A[B]T[C][B]dA；单元质量矩阵[712]
[M]=t∫Aλ〈N〉T〈N〉dA所用流量向量[712]
{Q}=qt∫L〈N〉TdL。

[BT3][STHZ]2.3渗透性函数
地球表面的陆地有相当大的一部分是处在干旱或者半干旱区的[13]，所以在实际工程中经常会遇到非饱和土的情况，非饱和土与饱和土有很多不同之处，因为在非饱和土中存在着负的孔隙水压力，使得非饱和土产生了独特的土力学性质。

非饱和土的两个最独特的土力学参数就是渗透性函数和体积含水量函数[1314]，如图2与图3所示，前者反映出渗透系数随着基质吸力的变化情况，后者是一种土水特征曲线，模块正是通过定义这两个参数，从而分析出非饱和土中的负孔隙水压力，以便进行非饱和问题的渗流分析。

稳态情况下不考虑孔隙水压力随时间的变化，因此无需定义体积含水量函数，相反瞬态情况则需要考虑。

SEEP/W2007提供了一种用体积含水量方程来估计渗透系数方程，而且提供了不同类型土的几种典型的含水量函数，在使用这些样本函数时，需要输入饱和含水量和残余含水量，再选择类似的土样，就可以很快建立非饱和土计算模型。

本文采用的是软件自带的Fredlund&Xing 函数模型。

饱和渗透系数和饱和含水量等相关实验参数是由延边[JP+1]水利设计院在《延边香水水利枢纽工程地质报告》中给出。

残余[JP]
图2不同材料的渗透性函数
图3不同材料的体积含水量函数
[JP2]含水量一般取饱和含水量的10%，基质吸力范围选为001～[JP]1 000。

SEEP/W2007首先通过筛分法或者沉降分析法[1314]，可以定义出粒径分布曲线，再由粒径分布曲线、饱和状态体积函数、体积压缩系数和给定孔隙水压力范围可以得到相应的体积含水量函数。

最后利用体积含水量函数，结合饱和状态渗透系数、孔隙水压力范围以及适当估算方法得到相应的渗透性函数。

由此可以看出渗透系数函数是在体积含水量函数的基础上派生出来的，非饱和区内渗透系数在不断变化，并且当越来越多的空气进入土体内部孔隙时，这种变化会更加明显，实际表现就是渗流会变得更加困难，进而土体就会出现负孔隙水压力。

[BT3][STHZ]2.4耦合方程
考虑时间增量的应力应变分析的有限元方程表示如下[16]：
∫[WTB1X]v[WTBX][B]T[C][B]d[WTB1X]v[WTBX]{a}=
b∫[WTB1X]v[WTBX]〈N〉Td[WTB1X]v[WTBX]+
p∫A〈N〉TdA+{Fn}[JY]（8）
式中：[B]为应变-位移矩阵；[C]为本构关系矩阵；{a}为节点增量的x，y位移列向量；〈N〉为插值函数的行向量；A为沿单元边界的面积；[WTB1X]v[WTBX]为单元体积；b为单位体力强度；p为表面压强增量；{Fn}为集中节点增加荷载。

[JP]对于这种时间增量分析，在每一个时间段上的增加荷重[1116]，SIGMA/W模块都能计算出与之相对应的位移增量，把这个时间段上计算出的所有位移增量全部叠加到上一时间段上，[JP]不过这种分析情况下只能使用一次单元单位质量力。

[JP2]式（6）与式（8）联立即得渗流场与应力场耦合分析控制方程[1116]。

在耦合问题的平衡有限元方程中[1116]，一边是外力荷载（既有体力又有外荷载作用），一边是内力荷载（既有有效应力又有孔隙水压力），方程两边的共同作用才能实现平衡。

[JP]
[BT3][STHZ]2.5极限平衡法
考虑弯矩平衡的安全系数方程[1116]：
Fm=[SX（]∑（c′βR+（N-μβ）Rtanφ′[]∑Wx-∑Nf±∑Dd[SX）][JY]（9）
考虑力平衡的安全系数方程[1116]：
Ff=[SX（]∑（c′βcos（α）+（N-μβ）tanφ′cosα）[]∑Nsinα-∑Dcosω[SX）][JY]（10）
式中：c′为有效凝聚力；φ′为有效摩擦角；μ为孔隙水压力；N为条块底部法向力；W为条块重量；D为线荷载；β，R，x，f，d，ω为几何参数；α为土体底部倾斜角
普通条分法Ordinary不考虑土条两侧面上的作用力[17]，在接近水平的滑动面上的法向和侧向底面剪切力很明显无法平衡，导致力多边形无法很好的闭合，可见这种方法得不出准确最佳的安全系数值。

Bishop 法不考虑土侧间剪切力，而且不满足所有的水平方向的力的平衡条件，只考虑力矩平衡，即土侧间法向力，以及除了水平方向的所有的力矩平衡条件；Janbu 法不考虑条间剪力，只考虑力矩平衡，而且仅满足所有水平力的平衡，其他方向并不满足；Morgensternprice（MP）法取决于土条间相互作用力的函数，可采用任何近似函数，并且同时满足力和力矩平衡[17]，考虑条间剪力和正应力；Finiteelement method考虑了稳定分析下的应力-应变关系[17]，满足了位移的兼容性，可以得到更加真实的应力应变分布情况。

[BT2][STHZ]3模型计算与分析
[BT3][STHZ]3.1建立模型
计算模型取高53 m，剖面长度218 m，该剖面模型共划分7 716个单元，7 892个节点，假定模型为线弹性，渗透性函数估算方法选取Fredlund&Xing模型，体积含水量函数估算方法选取Sample Functions，样本材料是黏土（Clay）。

模型材料的物性参数见表1。

计算上游水位5088 m，溢出边界压力水头为0 m，时间步长设为12步，步长增加为指数型，持续时间为400 d。

计算所得最大渗流速度为4468 8×106 m/s，不同时刻的最大渗流速度及渗透流量见表2。

计算12个时间步长，即达到400 d，从表2的数据可以看出，随着坝体受到荷载作用力，土体不断压缩，压缩减慢后，孔隙水压力消散也减慢，进而坝体渗流量逐渐趋于稳定。

说明随着时间推移[1820]，孔隙水压力逐渐消散减缓，坝体渗流量逐渐减小，并最终达到稳定状态。

[BT4][STHZ]3.2.2渗流场与应力场耦合分析
从表3的数据可知，不考虑耦合作用下坝体最大位移是1064 m，考虑耦合作用下坝体最大位移是1226 m，坝体位移最大处发生在坝体上游面与坝基接触间，即坝踵处，气孔玄武岩层之间。

观察水平方向位移图（图4、图5），可以发现不考虑耦合作用下最大位移为0179 4 m，发生在坝体上游坡面混凝土两侧石渣料处，考虑耦合作用后的最大位移为0206 2 m。

在坝体上游坡面，并且由坡面向心墙过渡，坝体的水平位移逐渐减小，坝体内部水平位移以坝轴线为界，上游坝体向下游移动，下游坝体向上游移动。

边坡位移在深层区变化缓慢，在浅层坡面位移变化比较大。

观察竖向位移图（图6、图7），发现不考虑耦合作用下最大位移是1064 m，发生在坝体顶部，考虑耦合作用下的最大位移为1226 m，且心墙下部位移分布明显不同，这是因为心墙下部位移变化较快，主要集中在心墙与气孔玄武岩层接触区域，因为在静水压力作用下，
[HJ1.8mm]上游坡面向下游移动，使得坝体向下沉降，覆盖层受到压缩，这是符合实际情况的。

3.2.3边坡稳定性分析
采用不同的极限平衡方法分别计算在考虑应力场与不考虑应力场下的上下游边坡抗滑安全系数。

通过对比表4和表5的数据可以看出，耦合后的上游边坡抗滑安全系数下降很多，表明上游边坡发生滑坡的可能性很大，这是因为上游边坡的表层在孔隙水压力及应力场作用下，很快达到饱和，此时边坡的基质吸力减小，从而使得土体抗剪强度降低，因此容易形成滑坡。

4结论
渗流场与应力场耦合分析下，边坡位移在深层区变化缓慢，在浅层坡面位移变化比较大，耦合工况下的位移最大值要大于不耦合工况下的最大值，而且在边坡抗滑移稳定计算中，耦合后的抗滑安全系数明显比不考虑耦合情况下的要小，更不利于坝体安全。

说明渗流场与应力场耦合分析在土石坝渗流稳定计算中更具有实际意义。

针对香水大坝的稳定分析表明，上游坡脚稳定性不够，应采用放缓边坡以及铺石压坡等方法提高边坡稳定性。

参考文献：
[1]李宗坤，王鹏飞，赵凤遥.基于流固耦合理论的土石坝稳定性分析[J].郑州大学学报：工学版，2009（3）：4447.
[2]王环玲，徐卫亚，晏鄂川.坝区裂隙岩体渗流场与应力场耦合分析[J].兰州理工大学学报，2004（5）：112115.
[3]孙兆涛，张剑，杜研岩等.基于GeoStudio软件的某尾矿坝边坡稳定性分析[J].露天采矿技术，2013（5）：3435。