摘要：本文使用刚体弹簧元模型来模拟拱坝-地基这类非连续介质和连续介质组成的系统，对一个悬臂梁的位移场和应力场的计算并与有限元方法计算结果的对比表明，这种方法具有良好的精度.应用这一模型对玛尔帕塞拱坝的溃坏进行了仿真，重现了拱坝的破坏过程，探讨了拱坝坝肩滑动失稳的原因，计算结果表明坝肩岩体裂隙和断层的强度不足是导致整个拱坝溃毁的关键因素.

关键词：三维刚体弹簧元 拱坝地基破坏 拱坝 仿真

　　近年来，用非连续介质力学模型研究高坝地基、高边坡与地下工程的变形稳定有长足进展，离散元[1～4]、dda[5]、刚体弹簧元[6]等的发展提供了求解非连续介质力学问题的有效工具.关于高坝-地基系统的破坏机理与仿真分析，由于问题的极端复杂性，迄今仍未见诸报导，本文利用改进的刚体弹簧元研究玛尔帕塞拱坝溃坏的机理与过程，是试图用连续-非连续介质统一的数值模型研究高坝-地基系统破坏过程的尝试.

法国玛尔帕塞拱坝1959年12月2日的溃坝事故是坝工史上的重大事件[7,8]，引起了坝工界的极大重视.世界各国专家和学者从未停止过对玛尔帕塞拱坝失事原因的探究，并提出了各自见解[9]，归纳起来有以下3种看法：(1)坝基变形，认为左岸坝基修建在弱风化顶层上，岩体弹性模量低，仅1000mpa左右，地基的过大变形使拱的推力向上向下转移，其中向上转移的推力使左岸重力墩超载而产生位移，导致拱坝失去支撑而破坏.(2)坝体上滑失稳，认为岩体浅层的裂隙在不利的条件下不能提供足够的抗滑力，致使坝体沿建基面向上和向下游滑动，拱圈拱弦被拉长，使坝沿一岸或两岸拱端转动，导致拱冠断裂，坝体崩毁.(3)坝肩岩体滑动，认为左岸坝基存在一个潜在的由下游断层、发育的节理面，以及下游临空面组成的滑动楔形体.由于地质条件在楔形岩体上游侧形成很高的静水扬压力使岩体滑移失稳.本文利用改进的三维刚体弹簧元对玛尔帕塞拱坝的失稳破坏过程进行仿真，对该坝的失稳机理作进一步的探讨.

1 拱坝破坏仿真模型

1.1 三维刚体弹簧元模型 假定三维刚体弹簧元的单元块体为不可变形的刚性体，在引入了刚度、强度等效准则后，将块体变形与强度指标统一集中于接触弹簧上.块体的运动符合刚体的平动和转动定律，块体与块体依靠相互之间连接的切向和法向弹簧和阻尼器来传递相互之间的作用，如图1.弹簧的作用力由块体间的相对位移确定.如图2，块体p对块体q的相对位移为δu，则他们之间的弹簧力变化量为

式中：ks、kn分别为切向和法向弹簧的刚度；δs、δn分别为相对位移增量的切向和法向分量；δs、δn分别为切向和法向弹簧力的变化量.

图1 接触模型

图2 相对位移增量示意

1.2 使用三维刚体弹簧元模拟非连续介质 三维刚体弹簧元模型基于离散元方法提出，是非连续介质模型的1种，由于块体之间可以错动、分离，因此可以方便地模拟非连续体如裂隙发育岩体的大变形、大位移问题.离散元方法中接触关系一般可分为角面、角边、角角和边边四类，关于各类接触关系的检索和其它离散元的详细内容已在cundall和鲁军等人的文章中有详细叙述，参见文献[1～4].

1.3 使用三维刚体弹簧元模拟连续介质

图3 刚体弹簧元模拟连续介质原理示意

图4 中心受拉杆件

1.3.1基本原理 如图3，通过将连续介质的变形和应力凝聚到单元之间连接的法向和切向弹簧上，可以用刚体弹簧元模型来模拟连续介质.即由弹簧的变形刚度(kn、ks)来等效连续介质的变形模量(e、g)，由弹簧力系统(fn、fs)来推求与之静力等效的连续介质的应力(σ、τ)分布，从而将刚体弹簧元系统和连续介质联系起来.(1)刚度等效.如图4所示，一连续的l×l1×l2方杆，弹性模量为e，泊松比为υ.将杆沿长度方向等分成n个l3×l1×l2的块体，其中l3=l／n，每个块体间的面面接触被视为四个角点处的点面接触，当n足够大，取弹簧的法向刚度kn=l1×l2×e/4l3时，使用三维刚体弹簧元计算它在中心作用力p下的变形与使用材料力学方法的结果是一致的.弹簧的切向刚度借用弹性模量和剪切模量的关系可取ks=kn/[2(1+υ)].单元应尽可能采用长方体，对于不规则的块体可以将其均化成长方体使用上述估算式，或另外推求相应的计算式.

(2)强度等效.如果刚体弹

簧元网格划分得比较密，则可以取单元面上每个弹簧力合力除以该面的面积求得其平均应力作为该面上的应力值.因为单元网格较密，可以认为当面上的平均应力达到材料的强度时该面上的弹簧全部断开，则单元在该面发生破坏.由于刚体弹簧元的计算工作量比较大，当为了节省计算工作量，网格划分得较稀疏时，可以按照材料力学方法假定单元面上的应力分布方式，由弹簧力根据静力等效原则求出该面上的应力分布.

(3)荷载等效.对于集中荷载，将荷载作用在相应的块体形心上，并将荷载对中心的力矩加入到该块体的合力矩中即与原荷载的作用等效.对于分布荷载也要等效成集中荷载加入到相应的块体中.

图5简单算例

1.3.2简单算例 如图5，一长40m，高4m，宽1m的悬臂梁受均匀荷载p=105n/m2作用，支座a固定，不考虑自重，坐标见图中所示.悬臂梁弹性模量e=1010n/m2，泊松比υ=03，密度ρ=2.4×103kg/m3，离散网格如图5(在厚度方向为一层，厚1m).使用上述的刚体弹簧元方法计算悬臂梁的变形和应力并与有限元的结果进行对比，图6为悬臂梁沿y轴的位移对比结果，图7和图8分别为x方向的正应力σx和剪应力τxy的对比结果.从图中可以看出两种方法计算的位移和应力结果除有限元法在固定端能反应局部应力集中外，整体应力具有很好的近似性.

图6 y方向位移等值线(单位：m)

图7 x方向正应力σx等值线(单位：mpa)

图8 剪应力τxy等值线(单位：mpa)

刚体弹簧元通过引入单元表面应力分布假设和变形近似等效来模拟连续介质，利用这些假设，刚体弹簧元通过上述的刚度等效、强度等效和荷载等效，能够把连续介质和非连续介质统一在一个模型之中，为研究坝肩失稳导致拱坝断裂破坏需要仿真连续与非连续耦合介质的大变形与破坏断裂的复杂过程提供了有力的工具.

2 玛尔帕塞拱坝失稳机理的研究

2.1 模型的建立

2.1.1 坝肩和坝体的离散和计算条件的设定 法国玛尔帕塞拱坝坝高66m，坝顶高程102.55m，坝顶长222.7m，长高比3.3，坝顶厚1.5m，底厚6.78m，厚高比为0.1.坝址岩体由带状片麻岩组成，含千节理倾向下游右岸，倾角30°～50°.该坝于1954年建成，蓄水历时5年至1959年12月蓄至最高水位约100m时，坝在短瞬间溃决，左岸与中部坝体全部冲走，右岸底部坝体残存.根据玛尔帕塞拱坝的地形和地质资料以及拱坝的体形和横缝设置将坝址处的岸坡简化成如图9所示，坐标系设置和拱坝离散网格如图10.在分析中边界条件设定为除“二面沟”内的岩体外，其他岩体均被固定.坝体、右岸坝肩和中间岩体依照连续介质进行模拟，坝体弹性模量取e=1010pa，坝体材料的抗拉强度取2.0mpa，抗压强度为20mpa，摩擦系数f=1.0；考虑到右岸坝基事后调查仍有坝体残存，表明该区岩体强度较高，故将坝体右岸及中间部分与基础的交界面强度较坝体提高1倍.左岸坝肩依照非连续体计算，重力墩和地基之间的摩擦系数设定为1.0，坝体和左岸坝肩之间以及“二面沟”内的岩体两组节理面与水平层面取相同的摩擦系数f1，根据不同的工况选取，上、下游断层摩擦系数为f2，亦根据不同的工况进行选择.计算中未考虑粘结强度c的作用.主要荷载为库水作用，稳定分析中库水水位均取100m，考虑左坝肩扬压力时，按静水压力作用施加在沿上游面坝基线的裂隙处，作用面见图9中虚线段所示.为了简化计算，所有荷载均一次施加.计算中阻尼参照cundall提供的阻尼模式选取刚度与质量成比例的rayleigh阻尼型式，阻尼系数α、β由试算确定，试算表明：如不出现过阻尼或阻尼过小情况计算结果均收敛于一稳态解的近似值，详见参考文献[3].

2.1.2位移分析 计算拱坝位移时左坝肩暂按连续介质考虑，并且不考虑扬压力作用，得到拱坝在1955年9月(水位79.75m)蓄水至1959年7月(水位94.19m)拱坝产生的位移，并与相应实测位移值进行比较，图11列出了顶拱和拱冠梁处的比较结果.可以看出二者最大值基本接近，但在拱冠梁坝基附近，实测的位移值较计算大1cm左右，顶拱计算结果较实测结果略偏大.这一比较结果表明选择的弹簧刚度参数kn、ks以及单元等效模型基本符合实际.

2.2 玛尔帕塞拱坝稳定性分析及失稳过程的仿真 通过改变坝体和左岸坝肩之间以及左岸坝肩“二面沟”内岩体的摩擦系数f1和上下游断层面的摩擦系数f2可以模拟和分析拱坝在各类参数条件下失稳破坏的情形.

图9岸坡离散

图10刚体弹簧元离散网格

图11拱坝位移计算结果与实测结果比较

2.2.1工况1——不考虑扬压力的浅层滑动 不考虑左岸坝肩处的扬压力，将坝体和左岸坝肩之间以及左岸坝肩“二面沟”岩体的摩擦系数f1和上、下游断层面的摩擦系数f2设为相同，从1.0开始折减，逐步计算，直到拱坝发生破坏.计算结果表明当摩擦系数f1和f2折减到0.6时坝肩失稳，拱坝破坏.图12列出了左坝肩附近的块体1(见图10)的位移过程.f=0.7时，块体位移量陡增，接近30cm，此时拱坝处于失稳的临界状态，f=0.6时，块体位移不再收敛，开始沿着坝肩向下游、向上滑动，整个拱坝随之破坏.

图12 块体1的位移过程(工况1)

图13 f=0.6时左坝肩重力墩位移过程(工况1)

分析图13～图15，拱坝的整个破坏过程可以分为两个阶段，第一阶段是大约3s以前拱坝整体位移，第二阶段是3s以后由于过大的变位使拱坝的工作方式发生改变，应力恶化，坝体溃毁.描述如下：在第一阶段，拱坝利用拱的作用将水压力转化为水平作用在坝肩上的推力，沿坝肩向上的分量使拱圈产生向上滑动的趋势，由于左坝肩岩体的摩擦系数较小，不能提供足够的抗力，巨大的拱推力转移到了左岸的重力墩上，使重力墩沿着墩轴线偏向下游移动了近2m(见图13，3s前重力墩主要以平动为主)，同时整个坝体沿建基面向下游、向上滑动，带动坝体以右岸为轴向下游转动(由图14中y方向的位移分布，坝体向下游位移由右到左依次增大，显示以右坝肩为轴向下游转动)和向上旋转(见图14中z方向的位移分布，坝体上抬量由右到左依次增大，显示以右坝肩为轴向上转动).此时，整个拱坝基本上保持整体发生变位，同时有沿右岸转动趋势.

(a)y方向位移等值线

(b)z方向位移等值线

图14 f=0.6，3s时拱坝坝体位移分布(工况1，单位：m)

拱坝在经过了第一阶段后，拱圈拱弦拉长，逐渐被压扁，承载能力极度降低，应力发生恶化，坝体破裂，拱坝的破坏进入到第二阶段——破坏阶段.由于中、下部块体受水压力较大，拱坝坝体从该部位偏左岸发生断裂，向下游弯折，首先被冲出，随后其周围坝块跟着折向下游，上部坝块受水压力较小，发生翻转后随在中下部坝块之后冲出，拱坝彻底破坏.与重力墩连接的坝段在被向下游冲出的过程中又带动重力墩发生了向下游的位移和旋转(见图13中3s以后的位移，向下游移动1m左右，在平面上旋转约3°).图15给出了上述过程

一些典型时刻拱坝的状态，拱坝破坏后最后残留状态见最后一幅.

图15 f=0.6时拱坝破坏过程(工况1)

上述拱坝破坏过程的计算结果和拱坝失事后的调查结果基本相符：(1)f值的大小，“实有摩擦角为30°～35°”，本文计算表明f值取0.7时拱坝处于临界状态，对应的摩擦角为35°，当拱坝左坝肩的实有平均摩擦角小于该值时，拱坝会发生上滑失稳，与拱坝破坏事实相符；(2)坝体的残留部分，拱坝失事后“右岸部分的坝体和中间部分的坝体底部基础留存了下来”，坝基由于嵌固作用，一般不易破坏，当拱坝上滑失稳时，应力发生恶化，拱坝沿横缝和工作缝等薄弱处断裂，同时底部的坝体得到卸载，从而保留了下来，计算中将坝体右岸及中间部分与基础的交界面强度较坝体提高了1倍，坝体的残留部分如图15，和拱坝破坏后的结果相符；(3)重力墩的位移，“左按重力墩在沿墩轴线方向上已移动约2m，沿垂直于墩轴线方向，向下游移动了约80cm，并略有倾侧”，本文计算的过程解释了上述重力墩先沿轴线方向移动(拱坝整体滑动时的推力所致)，又向下游移动并转动(拱坝破坏后，重力墩附近的坝段向下游破坏时的带动)的现象，位移量也基本接近；

(4)残留部分的位移规律，“ab坝段的位移很小，各坝段的位移值自右岸向左岸逐渐增大，至jk坝段处增达80cm；位移方向与坝轴线略呈斜交，并沿上游至下游方向略偏向于左岸.”，“左右岸，建筑物上游的岩石从混凝土与基岩接触点的几十cm处发生破裂，并有一条张开裂隙从底部沿岸坡脚向上延伸至接缝b处；这条裂隙在底部宽40～50cm，向上续渐减小.在右岸坝的下游可以看到若干在建筑物发生位移时所压碎的岩石区”，由于本文计算原来主要希望研究左岸坝肩滑动安全，未对中间坝基和右岸坝肩做离散和参数研究，因此未能计算出上游的裂缝和残留坝段80cm的位移量，但3s时的位移分布可以看出整个拱坝沿右岸的旋转和向下游的倾侧现象，这和上述调查结果相吻合，是产生上述裂隙和位移的原因，如果坝基岩体强度不够，整个拱坝沿右岸的旋转以及沿坝基向下游的倾侧必然产生很大的对上游面岩体的拉力和对下游面岩体的压力，使其开裂或破坏.

2.2.2 工况2——考虑扬压力时的浅层滑动在左岸坝肩处施加扬压力，重复工况1的计算.其结果与工况1基本一致，f=0.7时，拱坝处于失稳的临界状态，f=0.6时，拱坝破坏.表明扬压力对于浅层滑动的情况影响较小.

2.2.3 工况3——不考虑扬压力的深层滑动取坝体和左岸坝肩之间以及左岸坝肩“二面沟”岩体的摩擦系数f1为1.0并保持不变，在不考虑左岸坝肩处的扬压力情况下，将上、下游断层面的摩擦系数f2从1.0开始折减，逐步计算，直到拱坝发生破坏.计算结果表明：直到f2折减到0.1时坝肩仍然保持稳定，拱坝没有溃毁.因此，如果没有扬压力的作用，拱坝发生深层滑动的可能性极小，依靠岩体自身的重量基本上就可以保持稳定.

2.2.4 工况4——考虑扬压力时的深层滑动在左岸坝肩处施加扬压力，重复工况3的计算.计算结果表明当上、下游断层的摩擦系数f2折减到0.3时坝肩发生滑动失稳，拱坝破坏，见图16.分析图17、图18，与拱坝发生浅层滑动相对应，整个破坏过程也可以分为拱坝整体位移和坝体溃毁两个阶段.其过程如下：在第一阶段，作用在坝肩上的拱端推力向坝肩岩体深层传递，在最不利的滑动面(下游断层面)上，坝肩推力沿该面向上的分量和作用在坝基岩体上的扬压力有使坝肩岩体沿滑动面产生滑动的趋势，当岩体自重和滑动面上的抗滑力不足以阻止岩体的滑动时，整个坝肩岩体将沿着断层面产生滑动，巨大岩体的滑动带动整个左坝肩和重力墩运动(见图18中6.25s时的状态，坝肩岩体向上和下游滑动，重力墩折断破坏)，使重力墩发生平动和旋转(见图17，625s前重力墩的旋转)，拱坝坝体左边坝肩部位被抬起并推向下游，带动拱坝整体沿着右岸发生向下游转动.以后的拱坝运动、破坏的过程和工况1十分类似，不再赘述.

图16 块体1的位移过程(工况4)

图17 f2=0.3时坝肩重力墩位移过程(工况4)

深层滑动的计算结果和拱坝失事后的调查对照如下：(1)f2值的大小，“下游断层面中充填的粘土(摩擦角)的试验值只有30°”，本文计算的结果表明f2值取0.4时拱坝处于临界状态，对应的摩擦角为22°，当拱坝左坝肩上下游断层的平均摩擦角小于该值时，拱坝会发生深层滑动；(2)坝体的残留部分，深层滑动计算结果亦与其相似.(3)重力墩的位移，深层滑动计算的过程是重力墩一直在发生转动，稳定后在平面内的位移量和转动与实际调查结果接近；(4)残留部分的位移规律，这一现象和浅层滑动相同，不再重复.

图18 f2=0.3拱坝破坏过程(工况4)

3 结论

(1)计算表明玛尔帕塞拱坝在坝肩岩体抗滑强度较低时两种滑动失稳模式均有可能发生.当左坝肩岩体节理、裂隙的摩擦系数小于0.7时拱坝将发生浅层滑动失稳，而左坝肩上、下游断层的摩擦系数小于0.4，且坝肩作用有扬压力时拱坝会发生深层滑动失稳，因此，似乎拱坝发生浅层滑动的可能性更大一些.

(2)对拱坝沿坝肩岩体的浅层滑动和深层滑动失稳破坏过程的仿真计算均较好的解释了拱坝失事原因，再现了拱坝溃坏的过程和坝址处的残留情况.浅层滑动失稳过程似乎与调查结果符合得更好，然而由于计算资料并不十分充分，考虑的因素也不可能很全面，另外在计算中有所简化，深层滑动过程与浅层滑动过程亦有很多相似，因此很难排除或肯定拱坝失稳确切是那一种滑动模式.

(3)扬压力对浅层滑动影响很小，它主要是由于库水荷载的作用产生沿坝肩的滑动，但扬压力对深层滑动起着重要的影响，如果没有扬压力，对于玛尔帕塞拱坝发生深层滑动失稳的可能性很小.

(4)两种滑动有很好的近似性，尤其是坝体本身的位移和转动以及破坏过程十分相似.两种滑动方式下拱坝的整个破坏过程似乎都可以划分为两个阶段，第一阶段是拱坝整体发生位移，围绕右坝肩作向下游和向上转动，拱坝仍然以拱的方式工作，第二阶段仍是由于过大的变位使拱坝的工作方式发生改变，应力恶化，坝体溃决.

(5)本文计算从另一个方面也证实了拱坝坝体本身具有很高的变形适应性.在坝肩位移量高达30cm时(如工况1中，摩擦系数为0.7时，坝肩部位块体位移量达30cm)，拱坝由于其超静定的结构形式仍然是稳定的，充分说明了拱坝的超载安全储备.正像柯因在失事后所说：“我相信坝本身并未破坏，它是因地基破坏而被推走的”[7].