水利工程论文-拱坝-地基破坏的数值模型与溃坝仿真.doc

水利工程论文-拱坝-地基破坏的数值模型与溃坝仿真摘要：本文使用刚体弹簧元模型来模拟拱坝-地基这类非连续介质和连续介质组成的系统，对一个悬臂梁的位移场和应力场的计算并与有限元方法计算结果的对比表明，这种方法具有良好的精度.应用这一模型对玛尔帕塞拱坝的溃坏进行了仿真，重现了拱坝的破坏过程，探讨了拱坝坝肩滑动失稳的原因，计算结果表明坝肩岩体裂隙和断层的强度不足是导致整个拱坝溃毁的关键因素.关键词：三维刚体弹簧元拱坝地基破坏拱坝仿真近年来，用非连续介质力学模型研究高坝地基、高边坡与地下工程的变形稳定有长足进展，离散元14、DDA5、刚体弹簧元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和鲁军等人的文章中有详细叙述，参见文献14.1.3使用三维刚体弹簧元模拟连续介质图3刚体弹簧元模拟连续介质原理示意图4中心受拉杆件1.3.1基本原理如图3，通过将连续介质的变形和应力凝聚到单元之间连接的法向和切向弹簧上，可以用刚体弹簧元模型来模拟连续介质.即由弹簧的变形刚度(Kn、Ks)来等效连续介质的变形模量(E、G)，由弹簧力系统(Fn、Fs)来推求与之静力等效的连续介质的应力(、)分布，从而将刚体弹簧元系统和连续介质联系起来.(1)刚度等效.如图4所示，一连续的ll1l2方杆，弹性模量为E，泊松比为.将杆沿长度方向等分成n个l3l1l2的块体，其中l3=ln，每个块体间的面面接触被视为四个角点处的点面接触，当n足够大，取弹簧的法向刚度Kn=l1l2E/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.4103kg/m3，离散网格如图5(在厚度方向为一层，厚1m).使用上述的刚体弹簧元方法计算悬臂梁的变形和应力并与有限元的结果进行对比，图6为悬臂梁沿y轴的位移对比结果，图7和图8分别为x方向的正应力x和剪应力xy的对比结果.从图中可以看出两种方法计算的位移和应力结果除有限元法在固定端能反应局部应力集中外，整体应力具有很好的近似性.图6y方向位移等值线(单位：m)图7x方向正应力x等值线(单位：MPa)图8剪应力xy等值线(单位：MPa)刚体弹簧元通过引入单元表面应力分布假设和变形近似等效来模拟连续介质，利用这些假设，刚体弹簧元通过上述的刚度等效、强度等效和荷载等效，能够把连续介质和非连续介质统一在一个模型之中，为研究坝肩失稳导致拱坝断裂破坏需要仿真连续与非连续耦合介质的大变形与破坏断裂的复杂过程提供了有力的工具.2玛尔帕塞拱坝失稳机理的研究2.1模型的建立2.1.1坝肩和坝体的离散和计算条件的设定法国玛尔帕塞拱坝坝高66m，坝顶高程102.55m，坝顶长222.7m，长高比3.3，坝顶厚1.5m，底厚6.78m，厚高比为0.1.坝址岩体由带状片麻岩组成，含千节理倾向下游右岸，倾角3050.该坝于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列出了顶拱和拱冠梁处的比较结果.可以看出二者最大值基本接近，但在拱冠