摘要：在开裂单元中引入不连续形函数，然后基于虚功原理推导出开裂单元的广义刚度矩阵和广义荷载的具体表达式，以及有限元平衡方程，从而建立了一种三维裂缝扩展的不变网格有限元分析方法。运用该方法对已建的二滩拱坝和待建的小湾工程拱坝在相同的工况下进行了对比计算。计算结果表明：在目前所考虑的影响因素下，小湾拱坝和二滩拱坝的开裂范围、开裂程度及对网格疏密、混凝土抗拉强度和水力劈裂等影响因素的敏感程度都非常相似。考虑到二滩工程已建成，且运行正常，可以初步判定小湾工程的拱坝坝踵开裂危险性不大。
关键词：拱坝 开裂 有限单元法
坝踵开裂是高拱坝设计中需要考虑的一个重要问题。奥地利的科恩布赖拱坝建于1977年，坝高200.2m，在蓄水前坝下游底部水平施工缝张开，1978年蓄水接近满库时上游坝踵开裂穿通廊道，引起廊道内漏水，后来进行了十几年的漫长修补加固，并在坝的下游坝趾处建造了高约70m的混凝土支撑体。俄罗斯的萨扬舒申斯克拱坝建于1987年，坝高242.0m，施工期产生大量裂缝，当1990年水位首次达到正常水位时，河床坝段上游坝面靠近坝基处产生水平裂缝，引起廊道内漏水。这些高拱坝的损坏事故引起了世界各国坝工界的极大关注。
研究拱坝开裂分析的主要方法可以分为模型试验和数值计算两大类。可用于拱坝开裂分析的数值计算方法很多，如拱梁分载法、边界元法、流形元法、无单元法、有限单元法等[1～7]，其中有限元法是较为成熟且广为应用的数值方法。但有限元法在本质上是一种连续介质力学方法，必须加以改进才能被运用来分析裂缝扩展这样的不连续问题。改进的方法可以概括为两类：变网格法和固定网格法。变网格法随着裂缝的扩展而重新划分网格[8，9]，这种方法的优点是不需要开发新的有限元软件，但是存在计算量大和前处理难等问题。固定网格法则是保持有限元网格不变，通过修改开裂单元的插值关系和本构关系来反映裂缝的影响[10～12]。相比较而言，固定网格法的应用更为方便。本文在开裂单元中引入不连续形函数，然后基于虚功原理推导出开裂单元的广义刚度矩阵和广义荷载的具体表达式，以及有限元平衡方程，从而建立了一种三维裂缝扩展的不变网格有限元分析方法。本文方法在小湾工程拱坝的分析中获得了应用。由于拱坝开裂的影响因素复杂，计算方法也多种多样，很难做出完全定量的分析评判，且目前也没有相应的开裂深度安全准则，故研究中采用“校准法”的思路：取一个与小湾工程相类似的已建且运转正常的工程，在相同的工况下，采用同一计算方法进行对比分析，通过开裂深度与坝底宽度比值的比较，评价小湾高拱坝的开裂危险性和安全度。根据设计部门的建议，作为校准的已建工程选为二滩工程。
1 三维裂缝扩展的不变网格有限单元分析方法
1.1 开裂单元的位移模式[12] 常规有限单元法中，单元位移场{u}可以表示为：
|
{u}=[N]{u}e |
(1) |
|
|
|
|
|
|
式中：[N]和{u}e分别为单元形函数矩阵和单元结点位移向量。
式(1)表示的位移场在单元内是连续的，不能描述单元开裂以后沿开裂面的不连续位移模式。为求改进，需构造不连续位移场插值函数。
设单元e被开裂面分为e+和e-两部分(图1)，它们的位移场分别记为{u+}和{u-}，可以利用原单元e的形函数矩阵表示为：
|
{u+}=[N]{a}e {u-}=[N]{b}e |
(2) |
式中：{a}e和{b}e是单元结点的广义位移向量。
令：
|
|
(3) |
则：
|
{a}e={c}e+{d}e；{b}e={c}e-{d}e |
(4) |
将式(4)代入式(2)可得：
|
{u+}=[N]{c}e+[N]{d}e；{u-}=[N]{c}e-[N]{d}e |
(5) |
或：
|
{u}=[N]{c}e+[N]H(x){d}e |
(6) |
式中：[N]H(x)称为不连续形函数矩阵，其中H(x)是一个不连续函数，定义如下：
|
|
(7) |
式(6)右端的第一项与式(1)的右端项在形式上是一致的，表示的是单元中的连续位移场，而第二项表示的则是由于裂缝的出现而附加的不连续位移场。根据式(6)，开裂单元的位移模式是在连续单元位移模式的基础上，增加了反映不连续位移的插值项。这些附加的插值项与单元结点对应，但其形函数不连续，故称之为开裂单元的广义结点。
开裂单元的应变可以表示为：
|
{ε}=[B]{c}e+[B]H(x){d}e |
(8) |
式中：[B]是应变矩阵。
应力可以表示为：
|
{σ}=[D]{ε}=[D][B]{c}e+[D][B]H(x){d}e |
(9) |
式中：[D]是弹性矩阵，它是单元材料的弹性模量和泊松比的函数。
引入不连续形函数和广义结点的概念后，开裂单元的位移模式和应变模式都可以看作是在连续单元的模式上增加了广义结点的影响项，从而使得固定网格下的裂缝扩展分析变得简便易行。
1.2 有限元平衡方程 不连续单元的虚功原理可以表示为：
|
|
(10) |
式中：v和s分别表示单元的体积和边界(包括裂缝面)；{qv}和{qs}分别是单元的体积力和边界面力(包括裂缝面上的分布力)。
将式(6)、式(8)和式(9)代入式(10)，整理可得：
|
[K]{u}e={f}e |
(11) |
|
|
(12) |
|
|
(13) |
|
|
(14) |
式中：[K]称为单元广义刚度矩阵，{u}e称为单元结点广义位移向 量，{f}e称为单元结点广义荷载向量。[K]由四种类型的子矩阵组成，它们的表达式分 别为：
|
|
(15) |
|
|
(16) |
{f}e由两部分组成，分别为：
|
|
(17) |
|
|
(18) |
其中式(15)和式(17)在形式上与常规有限元法是一致的。
考虑Ⅰ型裂缝的情况，且直接采用裂缝前沿附近单元高斯点的应力，根据最大拉应力准则判断裂缝是否扩展，并且假定开裂单元被裂缝面贯穿。开裂分析中裂缝的扩展将产生应力释放，可以按下式计算其等效荷载：
|
|
(19) |
式中的积分沿裂缝面进行。
1.3 水力劈裂作用的模拟 在拱坝开裂计算中究竟如何考虑水力劈裂的作用是一个复杂的问题。随着裂缝的扩展，水流渗入裂缝，在裂缝面上产生劈裂作用，同时也引起渗流场的改变。如果考虑渗流场与裂缝扩展的耦合作用，将会使计算过程变得非常复杂，甚至可能无法实现。因此，宜对水力劈裂的作用进行适当的简化。对于高拱坝坝踵的张裂缝，可以不考虑水流渗入裂缝引起的渗流场改变，而只考虑水流在裂缝面的劈裂作用，作用力大小近似取上游全水头压强。这样简化是偏于安全的。以上算法已在有限单元法程序CORE3实现，并且通过算例考核证实了其合理性[12]。
2.1 工程概况 小湾水电站位于云南省大理州的南涧县和临沧地区凤庆县的交界，是澜沧江中下游河段规划8个梯级电站中的第二级，电站装机4200MW，水库库容151.32×108m3。大坝为混凝土抛物线双曲拱坝，最大坝高292m，坝顶高程1245m，是即将开工的世界第一高拱坝。小湾拱坝拱冠梁顶部宽度12m，拱冠梁底部宽度72.91m，厚高比为0.250。
二滩水电站位于四川省攀枝花市境内，雅砻江下游，电站装机3300MW，年发电170亿度。二滩大坝为混凝土抛物线双曲拱坝，最大坝高240m，坝顶高程1025m，是已建的中国第一高拱 坝。二滩拱坝拱冠梁顶部宽度11m，拱冠梁底部宽度55.74m，厚高比为0.232。
2.2 计算条件 高拱坝在坝踵、坝趾和基础高应力区会出现应力集中，并在坝踵和部分基础内产生裂缝。随着裂缝的开展，库水会进入裂缝，产生水力劈裂作用，从而促使裂缝进一步扩展。如果坝体与地基连接很好，接触面具有很高的强度，开裂可能会深入地基，但是这种裂缝往往属于稳定扩展，对大坝不会构成太大的危害。一般情况下，坝体与地基的交接面附近为一弱面，此面在靠近上游坝踵处往往受一法向拉应力，从而导致裂缝沿坝体与地基的交接面或稍高处扩展，此种裂缝对拱坝安全的危胁最大。本研究将以沿坝体与地基交接面或稍高部位的裂缝扩展作为重点。
2.2.1 计算模型 计算模型的建立对三维裂缝问题的求解有显著的影响，网格划分密实，单元形状相对正态，能够充分提高求解精度。为了分析不同计算模型对求解精度的影响，并对小湾拱坝和二滩拱坝进行分析比较，共建立了四套有限元模型，其中小湾拱坝粗网格模型和精细网格模型各一套，二滩拱坝粗网格模型和精细网格模型各一套。
|
|
|
|
|
|
2.2.2 计算工况 分别对考虑水力劈裂作用和不考虑水力劈裂作用两种情况进行了计算。考虑水力劈裂作用时，在已开裂缝面上施加上游全水头。此外，对坝体混凝土的单轴抗拉强度分别取0MPa、1MPa和2MPa进行敏感性分析。以上情况共组合成24种工况。
2.3 计算结论 图4和图5分别为小湾拱坝和二滩拱坝建基面裂缝扩展范围示意图(混凝土抗拉强度1MPa，考虑水力劈裂)，图6和图7分别为在相同工况下小湾拱坝和二滩拱坝拱冠梁断面开裂后应力矢量图。
|
|
|
|
|
|
(1) 计算裂缝的扩展范围大小与网格形式有较大的关系。网格加密后，开裂范围普遍有所减小，小湾模型的开裂深度占底厚的百分比平均减小约8.68个百分点，二滩模型平均减小约7.14个百分点。
(2) 随着混凝土抗拉强度的增加，裂缝的扩展范围逐渐缩小。混凝土抗拉强度由0MPa提高到1MPa时，小湾模型的开裂深度占底厚的百分比平均减小约10.42个百分点，二滩模型平均减小约9.52个百分点。混凝土抗拉强度由1MPa提高到2MPa时，小湾模型的开裂深度占底厚的百分比平均减小约13.02个百分点，二滩模型平均减小约9.52个百分点。
|
|
|
|
|
|
(3) 在其它条件相同情况下，考虑水力劈裂后的裂缝扩展范围比不考虑水力劈裂时大，小湾模型的开裂深度占底厚的百分比平均增加约7.99个百分点，二滩模型平均增加约10.32个百分点。
(4) 混凝土抗拉强度取1MPa时，小湾拱坝的最大绝对开裂深度为31.90m，最大相对开裂深度(开裂深度:坝基宽度)为0.4375；二滩拱坝的最大绝对开裂深度为25.73m，最大相对开裂深度(开裂深度:坝基宽度)为0.4615。
(5) 在目前所考虑的影响因素下，小湾拱坝和二滩拱坝的开裂范围、开裂程度及对网格疏密、混凝土抗拉强度和水力劈裂等影响因素的敏感程度都非常相似。考虑到二滩工程已建成，且运行正常，可以初步判定小湾工程的拱坝坝踵开裂危险性不大。
3 结语
本文介绍了三维裂缝扩展的不变网格有限元分析方法。拱坝开裂分析结果的影响因素很多，其中荷载、边界条件、开裂力学模型与参数、计算方法等都很重要却又很难准确把握。而且，作为一种新方法，配套的安全准则取值也需进一步研究。基于拱坝开裂问题的复杂性，在对小湾工程拱坝进行应用研究时，取已建的二滩工程为校准工程，在相同的工况下进行对比 计算，比较其裂缝扩展范围与坝基宽度的相对比值，从而可对小湾高拱坝的开裂危险性和安全度进行评价。
能对拱坝裂缝扩展进行三维分析的方法很多。笔者体会，包括本文介绍的方法在内，各种模型和方法都尚不够成熟以可靠地模拟裂缝发生位置、发展过程和的真实安全度。以下几个方向的研究可能会有助于问题的解决：(1) 考虑真实荷载与边界条件的施工和运行过程仿真模拟；(2) 在拉剪和压剪状态下混凝土及岩体的开裂机理和准则(包括起裂条件、开裂长度与方向)以及对应的参数；(3) 任意形状与分布的三维裂缝的定位与追踪技术；(4) 开裂深度安全准则。
参 考 文 献：
[1] 朱伯芳，栾丰.拱与梁产生裂缝后的失效角[J].水力发电学报 ，1997，(3)：55-60.
[2] 杨海霞，杜成斌，王德信.拱坝非线性开裂分析的分载位移法[J].水力发电，1997，(7)：19-22.
[3] Linsbaurer H N， Ingraffea A R. Simulation of caracking in large archdam. Part 1 [J]. J. of Struct.Eng.， 1989， 115 (7): 1599-1615.
[4] Linsbaurer H N， Ingraffea A R. Simulation of caracking in large archdam. Part 2. [J]. J. of Struct. Eng.， 1989， 115 (7): 1616-1630.
[5] 石根华.数值流形方法与非连续变形分析[M].北京：清华大学出版社，1 997.
[6] Belytschco T， Lu Y Y， Gu L. Element-free Galerkin Methods [J]. Int. J. for Num.Meth. in Eng.， 1994， 37: 229-256.
[7] 冠晓东，周维垣.应用无单元法近似计算拱坝开裂[J].水利学报，2000，(10)：28-35.
[8] Dhondt G. Automatic 3D mode I crack propagation calculations with finite elements [J]. Int. J. for Num. Meth. in Eing.， 1998， 41: 739-757.
[9] 王光纶，张楚汉，王少敏.混凝土重力坝的非线性断裂分析[J].水利学报，1997，(增刊)：182-187.
[10] Marek Klisinski. Finite element with inner softening band [J]. J. of Eng. Mech.， 1991， 117 (3): 575-587.
[11] Moes N， Dolbow J， Belytschco T. A finite element method for crack growth without remeshing [J]. Int. J. for Num. Meth. in Eng.， 1999， 46: 131-150.
[12] 汪卫明，徐明毅，陈胜宏.三维裂缝扩展的不变网格有限元分析方法[J].岩石力学与工程学报，2001，20(增2)：1608-1612.