三维快速拉格朗日法进行小湾拱坝稳.doc

三维快速拉格朗日法进行小湾拱坝稳定分析 寇晓东，周维垣，杨若琼 清华大学 水利水电工程系 摘 要三维快速拉格朗日分析是基于三维显式有限差分法的数值分析方法，它可以模拟岩土或其他材料的三维力学行为，可以准确地模拟材料的屈服、塑性流动、软化直至大变形。由于采用显式差分法，它尤其适合于材料的弹塑性、大变形分析以及施工过程的模拟。本文首先介绍了三维快速拉格朗日分析的基本原理及其特点，然后用美国Itasca Consulting Group Inc.开发的三维快速拉格朗日分析程序FLAC-3D进行了小湾拱坝的应力变形分析和参数敏感分析。 关键词三维快速拉格朗日分析；显式有限差分法；大变形；小湾拱坝 收稿日期199-09-09 基金项目国家自然科学基金资助项目59909004 作者简介寇晓东1971-，男，陕西吴堡人，清华大学水利水电工程系讲师。 小湾水电站位于澜沧江中游与支流墨惠江汇口下游，是澜沧江中下游河段开发的龙头水库，是以发电为主的综合利用水利枢纽。其主要建筑物有高292m的抛物线型双曲拱坝，岸边泄洪隧洞，岸边全地下式引水-发电系统。总库容为153亿m3，总装机容量为4200MW.枢纽区河谷走向近南北，地层为中、深变质岩系及第四系，分布岩石以黑云花岗片麻岩和角闪斜长片麻岩为主，加少量的约5片岩，岩石坚硬完整性尚好，但经受多期构造活动，断裂构造比较发育，以近东西向横河向、陡倾上游的单斜构造为主，断层也是近东西向较发育，如上游的F7，下游的F11、F5、F19、F23等断层均横贯两岸，右岸还有F10、F27、F22等近东西向断层，顺河向断层较少，仅在左岸坝肩下游有F20断层，这些构造将影响到大坝的稳定。此外坝肩还分布有蚀变岩体以及倾向河床的缓倾角裂隙。由于存在上述各种不利的地质构造，为确保大坝的安全，需要对小湾拱坝进行坝体和坝基的三维整体稳定分析。 三维快速拉格朗日分析是一种基于三维显式有限差分法的数值分析方法，它可以模拟岩土或其他材料的三维力学行为。三维快速拉格朗日分析将计算区域划分为若干六面体单元，单元网格可以随着材料的变形而变形，即所谓的拉格朗日算法，这种算法可以准确地模拟材料的屈服、塑性流动、软化直至大变形，尤其在材料的弹塑性分析、大变形分析以及模拟施工过程等领域有其独到的优点。鉴于此，本文将对三维快速拉格朗日分析的基本原理作一介绍并用美国Itasca Consulting Group Inc.开发的三维快速拉格朗日分析程序FLAC-3D对小湾拱坝的应力、变形及参数敏感性进行分析。 1 三维快速拉格朗日分析的基本原理 三维快速拉格朗日分析使用了如下三种计算方法[1]1离散模型方法。连续介质被离散为若干六面体单元，作用力均被集中在节点上。2有限差分方法。变量关于空间和时间的一阶导数均用有限差分来近似。3动态松弛方法。由质点运动方程求解，通过阻尼使系统运动衰减至平衡状态。 1.1 空间导数的有限差分近似 在快速拉格朗日分析中采用混合离散方法[2]，即，将区域离散为常应变六面体单元的集合体，而在计算过程中，又将每个六面体看作以六面体角点为角点的常应变四面体的集合体，应力、应变、节点不平衡力等变量均在四面体上进行计算，六面体单元的应力、应变取值为其内四面体的体积加权平均，六面体内四面体的应力应变第一不变量由该六面体内所有四面体的体积加权平均得到，然后重新得到四面体的应力、应变。这种方法既避免了常应变六面体单元常会遇到的位移剪切锁死现象，又使得四面体单元的位移模式可以充分适应一些本构的要求，如不可压缩塑性流动等。 如图1所示一四面体，节点编号为1到4，第n面表示与节点n相对的面，设其内任一点的速率分量为vi，则可由高斯公式得 图1 四面体 ∫Vvi,jdV∫SvinjdS 1 式中V为四面体的体积；S为四面体的外表面；nj为外表面的单位法向向量分量。 对于常应变单元，vi为线性分布，nj在每个面上为常量，由式1可得 2 式中上标l表示节点l的变量；l表示面l的变量。 1.2 运动方程 快速拉格朗日分析以节点为计算对象，将力和质量均集中在节点上，然后通过运动方程在时域内进行求解。节点运动方程可表示为如下形式 3 式中Filt为在t时刻l节点的在i方向的不平衡力分量，可由虚功原理导出。ml为l节点的集中质量，在分析静态问题时，采用虚拟质量以保证数值稳定，而在分析动态问题时则采用实际的集中质量。 将上式左端用中心差分来近拟，则可得到 viltΔt/2vilt-Δt/2Filt/mlΔt 4 1.3 应变、应力及节点不平衡力 快速拉格朗日分析由速率来求某一时步的单元应变增量 Δeij1/2vi,jvj,iΔt 5 式中导数可由式2近似。 有了应变增量，即可由本构方程求出应力增量，各时步的应力增量叠加即可得到总应力，在大变形情况下，还需根据本时步单元的转角对本时步前的总应力进行旋转修正。然后即可由虚功原理求出下一时步的节点不平衡力，进入下一时步的计算，其具体公式这里不再赘述。 1.4 阻尼力 对于静态问题，在式3的不平衡力中加入了非粘性阻尼，以使系统的振动逐渐衰减直至达到平衡状态即不平衡力接近零.此时式3变为 6 阻尼力为 filt-α|Filtsignvil7式中α为阻尼系数，其默认值为0.8. 8 1.5 计算循环 由以上可以看出快速拉格朗日分析的计算循环如图2所示。 2 三维快速拉格朗日分析的特点 无论是静力还是动力问题，三维快速拉格朗日分析都利用动态的运动方程进行求解，这使得三维快速拉格朗日分析很容易模拟动态问题，如振动、失稳、大变形等。 图2 计算循环 三维快速拉格朗日分析采用显式方法进行求解，对显式法来说非线性本构关系与线性本构关系并无算法上的差别，对于已知的应变增量，可很方便地求出应力增量，并得到不平衡力，就同实际中的物理过程一样，可以跟踪系统的演化过程。此外，显式法不形成刚度矩阵，每一步计算所需计算机内存很小，使用较少的计算机内存就可以模拟大量的单元，特别适于在微机上操作。而且在求解大变形过程中，因每一时步变形很小，可采用小变形本构关系，这就避免了通常大变形问题中推导大变形本构关系及其应用中所遇到的麻烦，使它的求解过程与小变形问题一样。当然，由于算法自身的原因也造成了三维快速拉格朗日分析的一些固有缺陷，与其他方法如有限元相比，快速拉格朗日分析的计算效率较低，但随着计算机运算速度的加快，这与三维快速拉格朗日分析的优越性相比，已经是微不足道。 根据前述原理，美国Itasca Consulting Group Inc.开发了三维快速拉格朗日分析程序FLAC-3D，该程序具有强大的前后处理功能，能很好地模拟地质材料在达到强度极限或屈服极限时发生的破坏或塑性流动的力学行为，特别适用于分析渐进破坏和失稳以及模拟大变形。它包含了10种弹塑性材料本构模型，有静力、动力、蠕变、渗流、温度5种计算模式，各种模式间可以互相耦合，可以模拟多种结构形式，如岩体、土体或其他材料实体，梁、锚元、桩、壳以及人工结构如支护、衬砌、锚索、岩栓、土工织物、摩擦桩、板桩、界面单元等，可以模拟复杂的岩土工程或力学问题。 3 小湾拱坝的计算条件 3.1 体型 本次计算体型为抛物线型双曲拱坝，详见表1. 表1 小湾抛物线型双曲拱坝体型参数 拱冠梁 右拱端 左拱端 高程 Yc/m Tc/m 拱端厚度 半中心角 拱冠曲率半径 拱端厚度 半中心角 拱冠曲率半径 Tar/m φr/ Rcr/m Ta1/m φ1/ Rc1/m 1245 6.500 13.000 15.000 41.529 433.588 14.000 41.179 480.113 1210 5.446 27.084 35.725 43.998 379.036 27.520 42.881 402.737 1170 -15.016 38.000 46.600 45.777 341.115 39.000 43.856 333.045 1130 -20.391 44.904 58.253 46.584 294.264 49.641 44.253 280.728 1090 -21.735 49.378 65.311 46.264 246.864 60.816 43.994 227.868 1050 -19.213 53.000 71.000 44.757 211.590 69.000 42.933 195.094 1010 -12.990 57.352 72.147 40.673 179.797 71.465 38.859 163.829 975 -4.639 63.000 72.000 31.318 159.426 72.000 28.535 160.000 953 1.946 67.930 70.108 12.821 131.826 70.464 11.075 153.272 3.2 计算参数及荷载工况 计算域内共采用20种材料，其力学特性见表2. 表2 小湾拱坝整体计算材料特性 代号 岩石类别 变形模量E/103MPa f C/Mpa μ γ/t/m3 1 混凝土 20.0 1.5 3.75 0.167 2.40 2 Ⅰ类岩石 25.0 1.5 2.3 0.220 2.63 3 Ⅱ类岩石 20.0 1.4 1.8 0.260 2.63 4 Ⅲ类岩石 14.0 1.2 1.2 0.270 2.63 5 ⅣbE1 2.0 1.0 0.5 0.280 2.00 6 ⅣA 4.0 1.0 0.6 0.280 2.0 7 Ⅴa 2.0 0.9 0.35 0.280 2.0 8 Ⅴb 2.0 0.9 0.35 0.280 2.0 9 F7 0.07 0.6 0.2 0.350 1.90 10 F11 0.50 0.8 0.3 0.350 1.90 11 F5 0.35 0.8 0.3 0.350 1.90 12 F10 0.35 0.8 0.3 0.350 1.90 13 F27 1.40 0.8 0.3 0.350 1.90 14 F19 0.35 0.8 0.3 0.350 1.90 15 F22 0.70 0.8 0.3 0.350 1.90 16 F23 0.35 0.8 0.3 0.350 1.90 17 F20 0.70 0.8 0.3 0.350 1.90 18 E4-5 2.0 1.0 0.5 0.280 2.0 19 E1 2.0 1.0 0.5 0.280 2.0 20 4 10.0 1.0 0.9 0.280 2.63 3.3 计算方案 计算区域为上游近一倍坝高，下游近三倍坝高，左右岸均大于两倍坝高，即在上下游1000m，左右岸1400m范围内模拟各类岩体及所含断层，包括横穿两岸的F7、F11、F5、F19、F23断层，以及右岸的F10、F27、F22横河向断层和左岸的F20顺河向断层。岩体模拟深度为一倍坝高，模拟至基础下Δ600m高程。基础上下游边界采用顺河向即y向位移约束，两坝肩顺河向边界节点采用横河x向位移约束，底部边界节点采用垂向z向位移约束。在1245-1130m高程区，坝体沿厚度分为2层，1130m-990m高程区分为4层，990m-953m高程区分为8层。计算采用FLAC-3D的内嵌语言FISH生成网格，并输入单元参数，计算本构模型用rucker-Prager准则。断层采用FLAC-3D的Interface来模拟。 图3 顺河向位移等值线单位m 计算中模拟了实际的筑坝过程，首先计算在没有坝体情况下岩基的自重场，然后将其位移置零，再计算筑坝后考虑坝体自重的自重场，再将其位移置零，接着计算水荷载作用下的位移应力场，这中间充分利用了FLAC-3D便于模拟施工过程的优点。 4 计算结果 4.1 位移计算结果 图3为水载作用下的顺河向位移等值线图，位移基本对称，最大约24cm，出现在拱冠梁顶部。图4为横河向位移等值线图，坝体基本向河床变形，最大位移约4cm左右。 4.2 应力计算结果 图5为最大主应力等值线图，图6为最大主应力等值线图的坝基视图，可以看出，坝体和基础基本上处于压应力状态下，只在坝踵部位有一部分拉应力区域，最大拉应力约为1.5MPa左右。 图4 横河向位移等值线单位m 图5 最大主应力等值线单位Pa 图6 最大主应力等值线坝基视单位Pa 4.3 安全度等值线 图7为水载作用下的安全度等值线图。可见，大部分区域的安全度都大于2.0，但在坝基部分区域安全度小于2.0，这是由于在这个部位有拉应力及较大的压应力。 4.4 参数敏感分析 为了研究岩体参数对拱坝稳定性的影响，进行了参数敏感分析，即将f、c值按比例降低，然后计算位移及安全度。表3列出了f、c值分别降至原来的70、50、30、26、25、24时坝顶拱冠上游面点顺河向的位移值。其中参数为25、24时的值为中间结果，因为此时坝体已失稳，达不到平衡状态，位移将随着计算步数的增加而增大。 图7 点安全度等值线 表3 坝顶拱冠上游面点顺河向位移 参数 70 50 30 26 25 24 位移/cm 24.30 24.5 35.0 226 3060 15400 图8 f、c降低到25时坝顶拱冠梁上游面点的位移历程 单位m 图9 f、c降低到24时计算过程中某时刻的位移等值线 单位m 图8为参数为25时坝顶拱冠上游面点随计算步数的位移变化图，由其趋势可看出此时拱坝已失稳，位移达不到平衡状态。图9为参数为24时的位移等值线图，从网格的扭曲形态可以看出此时发生了大变形，这是有限元不能得到的，这也正是FLAC-3D的特长所在。 图10和图11为f、c值分别降至原来的26、25时的屈服区图。可以看出，当参数降至26时，坝体出现了一大片安全度屈服区，当降到25时，坝体基本上整体屈服。 图10 f、c降低到26时的屈服区 图11 f、c降低到25时的屈服区 5 结论 从本次分析的结果来看，在水载作用下，坝体在坝踵区域会出现安全度小于2.0的区域，整个结构基本处于弹性状态。另外，由参数敏感分析可得出，拱坝的参数安全储备大致为4倍左右。三维快速拉格朗日分析由于采用了基于显式有限差分法的拉格朗日算法，使得它在分析岩土工程结构的弹塑性力学行为、模拟施工过程等方面有其独到的优点，尤其在稳定分析中，当发生塑性流动或失稳的情况下，三维快速拉格朗日分析可以很方便地模拟结构从弹性到塑性屈报、失稳破坏直至大变形的全过程，这是其他一些连续介质方法无法比拟的。 参 考 文 献 [1] Itasca Consulting Group, Inc.USA.FLAC-3D[M]。Fast Lagrangian Analysis of Continua in 3 Dimensions, Version 2.0, Users Manual. [2] Marti J P Cundall. Mixed Discretization Procedure for Accurate Modelling of Plastic Collapse[J]。Int.J.Num. Analy.Ethods in Geomech.,1982,6.