地震作用下土工格栅加筋高土石坝的稳定性分析

马 驰1，朱亚林，彭雪峰

(1.宿州学院 管理学院，安徽 宿州 234000； 2.合肥工业大学 土木与水利工程学院，合肥 230009)

1 研究背景

随着水资源开发利用进程的推进，我国还将兴建一批高土石坝，其中双江口、古水以及如美水电站坝高都将近或超过300 m。从高土石坝的分布位置可以看出，高土石坝主要分布在西部等高地震烈度区。地震发生导致边坡失稳破坏危害性极大且分布较广，通过5·12汶川地震资料发现，次生灾害中由地震滑坡造成的损失占地震所有损失的1/3[1]。

目前，分析地震作用下坝体稳定性问题的方法主要有拟静力法和完全非线性动力反应分析法等。其中，拟静力法应用比较广泛。 Belghali等[2]在加筋岩石斜坡的地震作用分析中采用拟静力法，栾茂田等[3]在土石坝边坡动力稳定性研究中分析滑移量时采用拟静力法，但都没有分析可能发生滑裂面的位置。朱亚林等[4]对面板坝的加固效果进行初步探讨，研究铺设土工膜和筋材的加固效果。周扬等[5]用大连理工大学开发的有限元软件GEODYNA计算出动力最小安全系数时程，研究土工格栅加固后的坝坡稳定性，计算时用钢筋代替格栅建立二维模型。这些研究通常用拟静力法代替动力时程，或是在二维模型的基础上分析。

地震作用力是和时间有关的作用力，所以采用非线性动力时程分析法是很有必要的。用非线性动力分析输入加速度波模拟地震作用，利用FISH语言的二次开发得到安全系数时程图，分析对比土工格栅的加固效果。在软件FLAC3D中二次开发得到计算安全系数时程的分析方法，通过该时程法可以更全面地分析土石坝动力稳定性。

2 三维模型的建立

2.1 计算模型

本文采用三维心墙坝，坝高200 m，上下游坝坡坡比均为1∶1.5，心墙上下游坡比均为1∶0.15，共分为主堆石区、反滤层和心墙3个区域。基本模型单元共47 000个，如图1所示。采用位移边界约束条件，模型底面约束水平向和竖向位移。

图1 坝体几何模型Fig.1 Geometric model of core wall dam

2.2 本构模型和材料参数

为准确合理地模拟高土石坝的动力反应分析特性，本文计算中堆石料本构模型采用修正的Mohr-Coulomb模型，在数值模拟计算中利用FISH语言对材料剪切模量做了修正[6]：

(1)初始剪切模量随围压的变化由公式(1)确定，计算在FLAC3D中通过FISH语言来实现，修正后坝体的剪切模量云图如图2。

(1)

图2 心墙坝体最大横断面初始剪切模量分布Fig.2 Distribution of initial shear modulus of the maximum cross section of core wall dam

(2)每一计算时间步，对Mohr-Coulomb模型的切线模量按照等效线性化模型的思路进行调整，即由当前时间步的剪应变参照模量衰减曲线得到模量衰减因子，用模量衰减因子乘当前时间步的切线模量。本文数值模拟计算时所用的剪切模量衰减曲线如图3。

图3 剪切模量衰减曲线Fig.3 Degradation curves of shear modulus

FLAC3D软件中的混合离散模型可以进行岩、土体的三维结构特性模拟和塑性流动分析。数值模拟时坝体材料静力参数见表1，动力参数见表2[7]。

表1 坝体的静力计算参数Table 1 Static parameters of dam

表2 坝体的动力计算参数Table 2 Dynamic parameters of dam

2.3 土工格栅

FLAC3D中土工格栅单元受力时主要由土工格栅材料本身的抗拉性能和土工格栅与网格之间接触面的摩擦力来承担。在FLAC3D中土工格栅单元体只抵抗薄膜类荷载(切应力和正应力)，不承受弯曲荷载。本文土工格栅将看成各向同性并且不会到达屈服极限的线弹性材料。在FLAC3D软件中铺设geogrid单元时，用sel geogrid确定格栅单元，geogrid单元是由一个交叉横条组成的粗网格，如图4(a)。孔宪京院士等[6]提出的高土石坝加固范围，在坝高1/4～1/5范围内坝顶部永久位移较大，所以在基本模型算例中土工格栅沿坝轴线对称布置，从坝高160 m开始，间距为4 m，共10层，铺设土工格栅效果如图4。

图4 铺设格栅效果Fig.4 Layout of geogrid

本文中土工格栅的弹性模量由15%应变下的抗拉强度通过线性关系计算得到，厚度和泊松比均按实际格栅尺寸进行数值计算。耦合切向刚度根据抗剪强度和面积的比得到，界面耦合摩擦角与土体摩擦角之间有个界面摩擦系数，一般取土体内摩擦角的0.85～0.92[8]，格栅参数见表3。

表3 格栅相关参数Table 3 Parameters of geogrid

2.4 水压力

FLAC3D中使用water table设置水面，会自动计算出静水压力分布图，如图5所示，上游水位设置在坝高180 m处，下游水位在50 m处[9]。

图5 孔隙水压力等值线Fig.5 Contours of pore pressure

3 坝坡的动力时程稳定性分析

对于坝坡的稳定性分析，计算方法多采用拟静力法。时程法计算坝体动力安全系数主要有2种方法：①事先假设滑裂面存在的位置对其进行分析；②从单元安全度的角度分析每个单元的稳定性[10]。

在已有文献的基础上利用单元屈服条件，在地震作用的每个时刻搜索潜在滑裂面求解安全系数。坝体的动力计算使用修正的Mohr-Coulomb模型。从强度破坏的角度来分析土体的稳定性[11]，根据单元受力状态，计算出单元安全系数。

根据研究采用Mohr-Coulomb准则，它的屈服函数为

f=σ1-σ3-(σ1+σ3)sinφ-2ccosφ=0 。(2)

可以写成

式中：σ1为第一主应力；σ3为第三主应力；c为黏聚力；φ为内摩擦角。

定义某单元i剪切破坏时安全系数Fsi计算表达式为[12-13]

(4)

式中：τfi可定义为单元i的抗剪强度，见式(5)；τi为i单元的最大剪应力，见式(6)。

(5)

(6)

式中ci,φi,σ1i,σ3i分别为单元i的黏聚力、内摩擦角、第一主应力、第三主应力。

图6 安全系数计算 示意图Fig.6 Schematic diagram of safety factor calculation

对于受地震作用的坝体来说，地震作用的每一个时刻坝体都会有相应的变化。也就是说每个时刻坝体发生变形，相应单元的应力都在变化。所以地震作用的每一时刻潜在滑裂面的位置都在随地震发生时间发生变化。本文计算安全系数的特点是不假定最危险滑裂面的位置，而利用FISH语言搜索每一时刻下游坝坡最大剪应变处的单元，这些单元连接形成的面，即为可能发生滑动的滑裂面。通过这个面上单元每一瞬时的动应力计算出坝坡的安全系数，安全系数计算如图6，地震作用时每个时刻都会搜索出一条滑弧，利用该滑弧上的单元并采用式(7)以加权平均的方式计算安全系数。

(7)

FLAC3D动力时程分析时，基于显示有限差分法，通过周围真实网格疏密度所得到的集中节点质量来求解运动方程，这个方程可以与结构单元耦合，所以能够模拟土界面与土工格栅的动力相互作用，这个方程还可以与流体计算耦合模拟坝体孔隙水压力。

模型计算选用的地震波以人工波为主，采用经拟合的糯扎渡坝址区100 a超越概率为2%的人造波[14]，顺河向峰值加速度为2.83 m/s2，竖向地震动输入为顺河向的2/3(如图7)。

图7 100 a超越概率2%地震加速度时程Fig.7 Seismic acceleration input with an exceedance probability of 2% in 100 years

图8 坝坡安全系数时程Fig.8 History of safety factor of dam slope

用上述计算安全时程的方法，选用前述的基本模型，分析格栅是否起到加固作用，结果如图8。由安全系数时程可以分析得到，无格栅时安全系数<1.0的时间累计量较多。地震波作用时，随着地震峰值加速度的到来，坝坡安全系数出现<1.0的情况，说明此时坝坡可能开始出现滑动状况。铺设格栅后，安全系数有所增加，但趋势较小，可见铺设土工格栅可以增加高土石坝的坝坡安全系数。可以看出从安全系数的角度分析格栅的加固效果是可行的。为得到合理铺设格栅的方法，达到最佳的加固效果，还需进一步研究。

4 格栅加筋效果分析

4.1 格栅弹性模量的影响

弹性模量是衡量土工格栅强度大小的关键指标，土工格栅弹性模量的强弱可以直接影响到土石坝的加固效果。以前述大坝为例，坝体材料参数不变，改变土工格栅弹性模量，寻找最合适的加筋土工格栅材料参数。弹性模量分别为0.25，2.50，25.00 GPa,得到的安全系数时程如图9。

图9 不同格栅弹性模量对坝坡安全系数的影响Fig.9 Safety factor histories of dam slope in the presence of varying elastic modulus of geogrid

土工格栅弹性模量为0.25 GPa时相比无格栅时安全系数增加，安全系数<1.0的时间开始减小；弹性模量为25.00 GPa时安全系数<1.0的时间开始消失，但安全系数增加幅度较小；弹性模量为2.50 GPa时安全系数明显增加，可以得到只有选择最合适的弹性模量才能起到提高安全系数的作用。建议弹性模量选择在2.50 GPa左右。

4.2 格栅长度的影响

以前述大坝为例，坝体材料参数不变，改变格栅的铺设长度，分别取1.3a和2.6a(a表示无加固措施情况下坝坡的滑裂面至边坡的距离)和满铺(由反滤层外缘向两岸延伸至坝坡处)。图10为不同格栅长度的安全系数时程，格栅弹性模量为2.5 GPa。

图10 不同格栅长度坝坡安全系数时程Fig.10 Safety factor histories of dam slope with different geogrid lengths

选择合适的弹性模量后，为确定格栅长度对坝体安全系数的影响，改变格栅长度后得到安全系数时程图。格栅长度越长安全系数就越大，格栅长度为1.3a和2.6a时提高幅度较小，但已经基本解决了安全系数较小的问题，此时安全系数<1.0的时刻还存在；只有满铺时安全系数相对较大，安全系数全部>1.0，满铺时坝坡安全性较好。

4.3 格栅间距的影响

以前述大坝为例，坝体材料参数不变，在坝高160 m处，分别间隔6，4，2 m铺设格栅，相应的格栅分别为7层、10层、20层，得到安全系数如图11。

图11 不同格栅间距坝坡安全系数时程Fig.11 Safety factor histories of dam slope with different geogrid spacings

从安全系数时程可以看出，不同间距下安全系数全部都>1.0，理论上坝坡是安全的。格栅铺设的间距越小，安全系数提高幅度就越大。格栅间距为4 m时安全系数增加幅度比格栅间距为2 m和6 m时都大，由于土工格栅单元并不抵抗弯矩，土工格栅过密并不能提高土石坝的相对刚度，所以过密铺设土工格栅效果并不是很明显，铺设格栅要选择合适间距。

5 实例验证

双江口水电站的心墙堆石坝的坝顶高程为2 510.00 m，最大坝高为312 m，坝顶宽度为16.00 m，坝顶长度为648.66 m。上游坝坡比为1∶2.0，下游坝坡比为1∶1.9。双江口心墙堆石坝设计成防渗体分区坝，坝体共分为堆石区、过渡层、2层反滤层和防渗体4个主要区域。其中，上游2层反滤层各4 m，下游2层反滤层各6 m；过渡层在反滤层和堆石体之间。防渗心墙顶高程为2 058.00 m，心墙顶宽为4 m，心墙上、下游坡比为1∶0.2。坝体参数见表4。

表4 双江口坝体材料的计算参数Table 4 Material parameters of Shuangjiangkou dam

自2008年汶川地震后，高烈度震区的土石坝采用土工格栅提高土石坝的抗震性能。双江口坝体的三维有限元网格共有单元82 360个，坝体采用六面体单元，三维模型如图12。双江口坝体网格及格栅铺设如图12(b)，对称铺设，从坝高约4/5处开始，每隔4 m铺设一层，格栅参数与基本模型的参数一致，详见表3。

图12 双江口坝体网格和加格栅效果Fig.12 Meshing and location of geogrid for Shuangjiangkou dam

大连理工大学双江口计算报告指出坝址区地震基本烈度为7度，70 a超越概率2%的基岩地震动峰值加速度为2.83 m/s2，数值模拟计算时采用的地震波加速度时程如图13[15]。

图13 70 a超越概率2%地震加速度时程Fig.13 Seismic acceleration inputs with exceedance probability of 2% in 70 years

计算结束后，对比有无格栅结果后处理得到坝体安全系数，如图14。

图14 坝体的安全系数对比Fig.14 Safety factor histories of dam slope in the presence and in the absence of geogrids

由图14可以分析得到，未加固时最小安全系数为0.447，安全系数<1.0的时间累计量为0.4 s；铺设格栅加固后最小安全系数均>1.0 s，安全系数也普遍增高。铺设土工格栅可以增加双江口土石坝的安全系数，与坝体报告结果一致。

6 结 论

本文采用水平和竖直2个方向的地震输入，用安全系数时程的分析方法，分析土石坝动力稳定性，并得到最佳的格栅铺设方法。具体结论如下：

(1)铺设土工格栅加固高土石坝后，通过动力安全系数时程分析，可得到地震作用的每一个时刻坝体安全系数变化，分析每个时刻坝体稳定性变化情况。

(2)通过安全系数时程图得到，改变土工格栅弹性模量,安全系数不同程度地增加，建议弹性模量选择在2.5 GPa左右；改变格栅长度后，长度越长，安全系数就越大，满铺时安全系数<1.0的情况已经消失，此时效果最佳；格栅铺设的间距越小，安全系数提高幅度就越大，但格栅铺设过密后加固效果不明显，格栅间距在4 m左右时最理想。

(3)通过坝高312 m的双江口心墙土石坝的有限差分数值分析，比较了土工格栅加固措施对坝坡动力稳定的影响。铺设土工格栅后，双江口土石坝坝坡安全系数提高，稳定性提高，土工格栅对坝体起到加固作用，增加了坝体抗震性能。并由此验证了土工格栅铺设方法的合理性和安全系数时程法分析的正确性。