2018, Vol. 51文章信息
- 丁建新, 邱焕峰
- DING Jianxin, Qiu Huanfeng
- 基于不连续有限元的拱坝应力应变仿真分析
- Stress-strain simulation analysis of arch dam based on discontinuous finite element method
- 武汉大学学报(工学版), 2018, 51(3): 198-204
- Engineering Journal of Wuhan University, 2018, 51(3): 198-204
- http://dx.doi.org/10.14188/j.1671-8844.2018-03-002
-
文章历史
- 收稿日期: 2016-10-24
2. 武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉 430072;
3. 中国电建集团贵阳勘测设计研究院有限公司，贵州 贵阳 550081
2. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China;
3. Power China Guiyang Engineering Corporation Limited, Guiyang 550081, China
某水电站位于云南西部的澜沧江中游河段上，是澜沧江中下游河段规划的8个梯级电站之一，也是中下游河段的龙头水库和巨型电站，是以发电为主兼有其他综合利用效益的水利枢纽.挡水建筑物为混凝土双曲拱坝，坝体于2007年11月在1 060 m高程检查廊道发现拱向裂缝，分析成果表明坝体中下部存在数量较多的裂缝并呈基本贯通趋势.图 1显示了典型坝段的裂缝分布，这些裂缝明显属于温度裂缝，裂缝的存在影响大坝的正常工作.因此了解水库蓄水过程中及正常运行后的坝体应力状态、尤其是已产生裂缝部位的坝体应力分布状态，是各界所关心的重要问题之一[1-3].
|
| 图 1 典型坝段裂缝分布(单位：m) Figure 1 Distribution of cracks on typical dam section (unit:m) |
对于拱坝系统，建立一套网格单元形态好、能同时反映裂缝、横缝、诱导缝、坝体廊道、坝基断层等等诸多结构的有限元模型十分困难.同时该拱坝监测系统的九向应变计及坝基压应力计均埋置在建基面附近，单元的形态直接影响计算结果与监测结果之间的对比分析，因此，研究降低网格剖分的难度，同时保证计算精度的方法非常重要.
为减少网格剖分带来的困难，无网格法、数值流行元、扩展有限元等新型数值方法均避免了繁琐的单元网格生成，具有广阔的发展前景.然而现阶段，以上新型数值算法的研究和应用尚不成熟[4-5].
文献[6]借鉴块体元和有限元耦合分析方法的思想，初步提出了大坝建基面模拟的不连续有限元分析方法.本文继续进行更深入的研究，并将其应用于拱坝应力应变仿真分析.前处理分别建立拱坝坝体及地基模型，两套网格相互独立，最后通过坝体及地基建基面设置虚拟不连续结构面实现坝体与地基网格之间的耦合，仿真中考虑结构面的粘结、滑移效应.通过改变上游贴角与地基之间的虚拟不连续结构面的法刚和切刚来实现二者接触面在蓄水之后自动张开和闭合，并采用隐含裂缝单元模拟坝体裂缝，基于不连续有限元法进行了拱坝应力应变的仿真分析.
1 基本原理思想不连续有限元法由殷德胜、汪卫明等提出[6-7]，其最初的目的是为解决水工结构以及岩土工程数值计算中常遇的断层、节理、缝等不连续结构的模拟问题.其基本思想是将不连续面两侧的结构块单独划分网格，而不受结构面结点需对应的约束，从而减小网格划分的难度.图 2所示为含有一条不连续面k的结构，常规有限元分析中，用有厚度或者无厚度的节理单元模拟该结构面[8].不连续有限元则对上、下两部分区域进行单独的网格划分，两侧网格间的结点可以不对应，网格间由虚拟设置的不连续面k联接.通过在普通有限单元法的位移插值、本构方程以及平衡方程中考虑该不连续面的影响，求解方程，得到系统的位移和应力结果.不连续有限元法减小了网格划分的难度，其实质是把复杂的前处理网格剖分转化为并行处理.
|
| 图 2 不连续网格示意图 Figure 2 Sketch of discontinuous meshes |
不连续有限元方法中，单元内任意点的位移和应变计算同普通有限元，而不连续面k上任意点对Pi和Pj的相对变位为
(1)
其中：J(k)是单元i相对于不连续面k的上下盘系数；Lk是不连续面k的局部坐标转换矩阵，可以通过不连续面k的倾向φk、倾角θk计算得到；其他符号意义可参见文献[6].
2.2 本构方程不连续有限元法中，普通单元的本构方程与常规的有限单元相同.借鉴无厚度节理单元的本构方程，不连续面的弹性矩阵为
(2)
文献[6-7]已对式(2)中kn、ks的取值进行了比较详细的讨论，本文不再赘述.
考虑不连续面的粘塑性屈服，不连续面的粘塑性应变率计算如下：
(3)
粘塑性屈服流动取Mohr-Coulomb模型：
(4)
式中：φk、ck分别为不连续面k的内摩擦角和凝聚力；γk为流动参数.
考虑不连续面的张开，直接采用不连续面高斯点的平均法向应力，根据最大拉应力准则判断其是否张开，不连续面张开后将产生拉应力释放，可以按下式计算其等效荷载[9-10]：
(5)
在多步荷载的时空仿真分析中，不连续面的开裂状态根据其应力状态进行动态调整.
2.3 平衡方程根据虚功原理，推导系统的平衡方程：
(6)
其中, 整体刚度矩阵K含3种类型的子矩阵，其表达式分别为
(7)
(8)
(9)
其中:
(10)
式(6)中:ΔU为整体位移向量；ΔF为整体荷载向量.
2.4 数值积分利用普通的高斯积分对式(7)进行积分计算，并采用精确的分片高斯积分方法对式(8)和式(9)进行积分计算，计算精度满足要求.
实际施工中，各坝段各浇筑块为跳仓浇筑，因此，坝体多个坝段在整个时空仿真的不同时间出生.地基单元则在时空仿真的每一步均参与计算，当其上浇筑块在该仿真时间步出生时，激活该浇筑块与地基之间不连续面的高斯点，不连续面参与计算.
3 拱坝时空仿真对于大体积混凝土结构的温度以及应力应变仿真分析，目前最常用的是朱伯芳院士的温度应力理论[11].各科研院校依据该理论在ANSYS、ADINA等软件上进行二次开发，或者自编相应的仿真分析软件.例如，水科院张国新课题组开发的SAPTIS[12]，武汉大学陈胜宏课题组开发的CORE3D软件等，在国内一些工程的温度应力仿真分析中发挥了重要的作用，例如小湾、溪洛渡等工程[13, 14]，形成了一套较为成熟的理论体系，指导了工程建设实践.
拱坝的时空仿真分析主要包括温度场、应力应变场和渗流场的时空仿真分析.三物理场之间相互耦合，目前受分析难度、求解规模等限制，三场之间一般进行解耦分析，各物理场之间进行数据的插值传递.首先进行温度场的仿真，提取模型结点的温度结果，然后进行应力应变场的仿真分析，对于渗流场，一般是考虑上下游水压力，或者提供渗透荷载.
3.1 温度场仿真考虑空间域和时间域不耦合，通过有限元对空间离散，差分法对时间离散，温度场仿真的有限元求解方程如下：
(11)
式(11)中各表达式的含义可参考文献[11]，本文不再细述.式(11)结合初始条件和边界条件，即可求解获得温度场的结果.
在商业有限元软件中一般都有相应的温度场求解模块.温度场方程的特点是：求解比较简单，结点自由度为1，并且温度为标量.但是方程(11)为非稳定温度场，需要时间的迭代；其次在仿真分析中，随着混凝土浇筑块的出生，温度场的边界在不断变化；另外水库水温以及混凝土施工浇筑中的多期通水冷却等也需要进行模拟，从而才能得到比较精确的温度结果，以提供给应力应变场仿真.
3.2 应力应变场仿真应力应变场每一时间步的求解方程如式(6)，累加时间步后即得到各单元应力如下：
(12)
应力场仿真的特点是：荷载复杂，对于式(6)右端的整体荷载向量需考虑温度荷载、自重、水压力或者渗透体积力、徐变荷载、混凝土的自生体积变形荷载、粘塑性荷载等等.式(6)结点自由度为3个方向的位移，同时需要计算单元高斯点的应变以及应力结果.此外，应力应变场往往需要建立精细的复杂有限元整体模型.
3.3 裂缝单元模型拱坝在施工期以及运行期，受各种物理荷载的作用，可能在坝踵以及坝体产生裂缝，裂缝的存在改变了坝体的位移和应力分布规律，裂缝尖端应力集中.目前，常见的裂缝模拟理论和方法有：扩展有限元、数值流行元以及无网格法.以上新方法在工程实际的裂缝模拟应用中尚不成熟，商业软件ANSYS、ADINA则是通过非线性接触单元来模拟裂缝[15].借鉴岩土工程节理单元的思想，经过数学推导，开发新的隐含裂缝单元模型来模拟坝体裂缝.该模型具有以下两条基本原则：
1) 应变迭加原则：单元的应变增量等于混凝土块与裂缝的应变增量之和.
2) 应力一致原则：单元、混凝土块和裂缝的应力增量相等.
拱坝的应力应变场时空仿真中依据坝中缝的应力状态改变缝的刚度和开合状态实现缝的复杂力学行为模拟.
4 数值算例应用利用以上理论，以云南澜沧江某拱坝为例，建立相应的模型，制定仿真的方案和技术路线，进行拱坝应力应变场的仿真分析.
4.1 仿真模型不连续有限元网格模型如图 3~5所示.其中地基单元总数为276 765、结点总数为188 896，坝体单元总数为886 897、结点总数为753 712.
|
| 图 3 坝体网格 Figure 3 Mesh of arch dam |
|
| 图 4 地基网格 Figure 4 Mesh of foundation |
|
| 图 5 坝体-地基系统耦合模型 Figure 5 Coupled model of arch dam and foundation |
整体坐标系：横河向为X轴，指向左岸为正；顺河向为Y轴，指向上游为正；铅直向为Z轴，向上为正.
边界范围：左右岸方向分别取约738 m和715 m，上下游方向取935 m，铅直向下取到高程650 m，模型顶高程左岸为1 477 m，右岸为1 626 m.
边界条件：左右岸边界面X方向的法向位移约束、上下游边界面Y方向的法向位移约束、底面Z方向的法向位移约束.
模型中考虑了坝体横缝、诱导缝和裂缝以及F7、F11、F5、F10、F20等断层及其他地质因素.
典型河床22号坝段的浇筑过程如图 6所示，根据设计单位提出的“分期蓄水、监测反馈、逐步检测、动态控制”原则，水库于2012年10月31日首次蓄水至正常蓄水位1 240 m，水库的蓄水过程如图 7所示.
|
| 图 6 22号坝段浇筑过程线 Figure 6 Pouring process of dam monolith No.22 |
|
| 图 7 水库蓄水过程线 Figure 7 Reservoir impoundment process |
计算模型中主要的岩体力学参数及混凝土参数如表 1和表 2所示.
| 编号 | 岩体类别 | 变形模量 E/MPa | 泊松比μ | 粘聚力c /MPa | 内摩擦 系数f |
| 1 | Ⅰ | 25 000 | 0.22 | 2.2 | 1.50 |
| 2 | Ⅱ | 20 000 | 0.27 | 1.8 | 1.45 |
| 3 | Ⅲa | 14 000 | 0.28 | 1.2 | 1.20 |
| 4 | Ⅲb1 | 10 000 | 0.29 | 1.0 | 1.15 |
| 5 | Ⅲb2 | 6 000 | 0.30 | 0.7 | 1.10 |
| 6 | Ⅳ | 5 000 | 0.32 | 0.6 | 1.00 |
| 编号 | 弹性模量 E/GPa | 抗压强度 σc/MPa | 抗拉强度 σt/MPa |
| C40 | E(t)= 32.6t(9.79+t) | σc(t)= 48.49t(12.27+t) | σt(t)= 3.28t(9.09+t) |
| C35 | E(t)= 30.2t(10.3+t) | σc(t)= 42.65t(13.33+t) | σt(t)= 3.11t(13.76+t) |
| C30 | E(t)= 29.0t(9.07+t) | σc(t)= 36.7t(13.66+t) | σt(t)= 2.95t(13.65+t) |
| 注：t为混凝土龄期, d. | |||
按照设计单位提供的大坝浇筑(图 6)和水库蓄水过程(图 7)进行仿真模拟，其中，温度场仿真计算的时间步长为1 d，应力应变场仿真计算的时间步长为2~3 d.首先计算出坝体的温度场，然后插值到不连续有限元系统进行应力应变场的仿真.
4.4 仿真技术路线常规的三维弹粘塑性有限元程序可以模拟大坝的浇筑与运行过程中的应力/渗流/温度耦合仿真分析，并能实现裂缝的开合迭代.开发不连续有限元法程序，并利用其对拱坝系统进行仿真分析，仿真分析中考虑虚拟不连续结构面的粘结、滑移效应，并设置上游贴角与地基之间的虚拟不连续结构面在蓄水之后自动开裂和闭合.
5 计算成果分析为验证本文所提方法的正确性，同时进行有限元仿真计算，选取了典型的结果和有限元进行比较.
5.1 位移分析从顺河向位移分布图 8以及铅直向位移分布图 9上可以看出，不连续有限元与有限元的分布规律基本一致，量值接近.
|
| 图 8 22号坝段顺河向的位移对比图(水位：1 190 m) Figure 8 Displacement distribution of dam monolith No.22 along the river direction (water level: 1 190 m) |
|
| 图 9 22号坝段铅直向的位移对比图(水位：1 190 m) Figure 9 Displacement distribution of dam monolith No.22 along the vertical direction (water level: 1 190 m) |
选取典型坝段的九向应变计及坝基压应力计的不连续有限元结果同有限元及监测资料进行对比，可以反映出不连续有限元应力的精度及变化规律.
从图 10和图 11可以看出，建基面附近的九向应变计A22-S9-02以及坝基压应力计A22-C-02的结果都接近有限元的结果，并且二者反映的规律同监测值一致，这说明了不连续有限元应力结果满足仿真计算的要求.
|
| 图 10 22号坝段九向应变计A22-S9-02计算与监测对比 Figure 10 Historic curves of circumferential stress of point A22-S9-02 |
|
| 图 11 22号坝段坝基压应力计A22-C-02计算与监测对比 Figure 11 Historic curves of normal compressive stress of point A22-C-02 |
当前仿真条件下，九向应变计的应力结果规律如下：混凝土块出生后，横河向应力逐渐受压；封拱灌浆前通水冷却，压应力逐渐减小；此后形成整体，又逐渐受压；导流洞下闸之后，水库逐渐蓄水，横河向压应力有所减小，此后应力受水位变动而变化.顺河向应力规律类似横河向.坝体铅直向应力则主要由于混凝土的自重作用，出生之后一直受压，并且随着坝体混凝土的施工浇筑，压应力逐渐变大.坝基压应力计的规律类似九向应变计的铅直向应力规律，数值也接近，主要原因是九向应变计位于坝基压应力计的上部，两测点高差相差不大.
5.3 裂缝状态为了了解当前的坝体裂缝工作状态，需要对裂缝面的安全性进行评价，裂缝面的抗剪安全度采用“抗力/效应”进行度量，其定义公式如下：
(13)
式中：ns为裂缝单元总数；Si为裂缝单元面积；裂缝面的物理力学参数f′=0.92，c′=0.40 MPa(对应70%充填及其中20%粘结良好情况).
图 12给出了1 210 m水位仿真条件下的裂缝面工作状态，从图中可以看出，大部分裂缝处于闭合状态，右岸裂缝有少量压剪屈服.坝体裂缝基本处于受压状态，随着库水位的上升，缝面压应力增大，坝体裂缝不会扩展，裂缝处于压紧、闭合、稳定状态.
|
| 图 12 裂缝工作状态图(水位：1 210 m) Figure 12 Cracks station (water level: 1 210 m) |
拱坝坝体自发现裂缝后及时调整温控措施，并对坝体裂缝进行化学灌浆等措施处理，后续的检测工作中未发现新的裂缝，坝体原裂缝也未扩展.2012年10月31日首次蓄水至正常蓄水位1 240 m，各方面的科研以及监测成果显示拱坝整体工作情况正常.
6 结论由以上分析可以看出：
1) 不连续有限元将传统的大坝-地基系统前处理网格划分转变为并行的两套网格划分，基本将原来的时间减半，由于不考虑两套网格之间的相互影响，网格划分的难度也大大降低，适合于拱坝-地基系统的仿真分析.
2) 不连续有限元及有限元模拟实际浇筑过程的结果在位移及应力的分布规律及数值上接近，并且二者反映的规律同监测值一致.
3) 当前仿真条件下，典型剖面坝体顺河向向下游变形，最大位移约为80 mm；铅直向则竖直向下变形，最大位移约为30 mm.坝体混凝土出生后横河向应力逐渐受压，封拱灌浆前通水冷却，压应力逐渐减小，此后形成整体，应力受水位变动影响；顺河向应力规律类似横河向；铅直向应力随着混凝土的施工浇筑，压应力逐渐变大.
4) 计算条件下，裂缝大多数处于受压闭合状态，右岸裂缝有少量处于压剪屈服状态.随着库水位的上升，坝体裂缝受压，处于稳定状态.
| [1] | Wang Weiming, Ding Jianxin, Wang Guojin, et al. Stability analysis of the temperature cracks in Xiaowan arch dam[J]. Science China, 2011, 54(3): 547–555. DOI:10.1007/s11431-010-4280-1 |
| [2] | Fu Shaojun, He Ting, Wang Guojin, et al. Evaluation of cracking potential for concrete arch dam based on simulation feedback analysis[J]. Science China, 2011, 54(3): 565–572. DOI:10.1007/s11431-010-4274-z |
| [3] |
中国水电顾问集团昆明勘测设计研究院. 高拱坝时空特性演化机理及监控体系研究报告[R], 2012.
HydroChina Kunming Engineering Corporation. The research report on time and space evolution mechanism and monitoring system of high arch dam[R], 2012. |
| [4] | Bettess P. More on infinite elements[J]. International Journal for Numerical Methods in Engineering, 1980, 15(11): 1613–1626. DOI:10.1002/(ISSN)1097-0207 |
| [5] | Gerdes K. A review of infinite element methods[J]. Journal of Computational Acoustics, 2000, 8(1): 43–62. DOI:10.1142/S0218396X00000042 |
| [6] |
殷德胜. 岩体渗流与应力分析的数值方法研究[D]. 武汉: 武汉大学, 2010.
Yin Desheng. Study on the numerical method for the seepage and stress analysis of rock mass[D]. Wuhan: Wuhan University, 2010. http://cdmd.cnki.com.cn/Article/CDMD-10486-2010166662.htm |
| [7] | Wang Weiming, Zhang Zhijun, Yin Desheng. Discontinuous finite element method for the stress analysis of dam/foundation system[C]// Proc. of the 6th International Conference on Dam Engineering, 2011. |
| [8] |
朱伯芳.
有限单元法原理与应用[M]. 北京: 中国水利水电出版社, 2009.
Zhu Bofang. The Finite Element Method Theory and Applications[M]. Beijing: China WaterPower Press, 2009. |
| [9] |
陈胜宏, 汪卫明, 徐明毅, 等. 小湾高拱坝坝踵开裂的有限单元法分析[J].
水利学报, 2003(1): 66–71.
Chen Shenghong, Wang Weiming, Xu Mingyi, et al. Analysis of crack propagation in high arch dam heel of Xiaowan project by finite element method[J]. Journal of Hydraulic Engineering, 2003(1): 66–71. |
| [10] |
汪卫明, 徐明毅, 陈胜宏. 三维裂缝扩展的不变网格有限单元分析方法[J].
岩石力学与工程学报, 2001, 20(增2): 1608–1612.
Wang Weiming, Xu Mingyi, Chen Shenghong. A finite element method for 3-D crack propagation without remeshing[J]. Chinese Journal of Rock Mechanics and Engineering, 2001, 20(s2): 1608–1612. |
| [11] | Zhu Bofang. Thermal Stresses and Temperature Control of Mass Concrete[M]. Butterworth-Heinemann Press: 2013. |
| [12] |
张国新. SAPTIS:结构多场仿真与非线性分析软件开发及应用(之一)[J].
水利水电技术, 2013(1): 31–44.
Zhang Guoxin. Development and application of SAPTIS - a software of multi-field simulation and nonlinear analysis of complex structures (Part Ⅰ)[J]. Water Resources and Hydropower Engineering, 2013(1): 31–44. |
| [13] |
邹丽春, 陈胜宏, 王国进, 等. "数值小湾"研究与探索[J].
水利学报, 2011, 42(增): 19–29.
Zou Lichun, Chen Shenghong, Wang Guojin, et al. Study and exploration of Xiaowan arch dam digitized system[J]. Journal of Hydraulic Engineering, 2011, 42(supp): 19–29. |
| [14] |
陆佑楣, 樊启祥, 周绍武, 等. 金沙江溪洛渡高拱坝建设的关键技术[J].
水力发电学报, 2013, 32(1): 187–195.
Lu Youmei, Fan Qixiang, Zhou Shaowu, et al. Key technologies for construction of Xiluodu high arch dam on Jinsha River[J]. Journal of Hydroelectric Engineering, 2013, 32(1): 187–195. |
| [15] |
张洪才.
ANSYS14.0理论解析与工程应用实例[M]. 北京: 机械工业出版社, 2013.
Zhang Hongcai. ANSYS14.0 Theory Analysis and Engineering Application[M]. Beijing: China Machine Press, 2013. |