混凝土重力坝深层滑动稳定性评估软件模块研发

崔书生，苏怀智，张帅，顾浩钦

（河海大学a.水文水资源与水利工程科学国家重点实验室;

b.水资源高效利用与工程安全国家工程研究中心;c.水利水电学院，南京 210098）

混凝土重力坝深层滑动稳定性评估软件模块研发

崔书生a，b，苏怀智a，c，张帅a，b，顾浩钦a，b

（河海大学a.水文水资源与水利工程科学国家重点实验室;

b.水资源高效利用与工程安全国家工程研究中心;c.水利水电学院，南京 210098）

基于混凝土重力坝深层滑动面上应力分布情况的数值计算成果，结合分项系数极限状态设计方法，研究了一种评估混凝土重力坝深层组合滑动稳定性的改进分项系数法。在此基础上，利用Delphi编程工具，结合AutoCAD二次开发技术对该方法实现了程序化，得到了评估混凝土重力坝深层滑动稳定性的可视化软件模块。该方法和软件模块，充分利用了有限元数值模拟方法在计算应力场和考虑坝基及坝体变形影响等方面的优势，又结合了可靠度设计方法，得出的稳定安全系数能够较好地反映重力坝的抗滑稳定特性，且沿用了设计规范中的Ks表示方法，易于计算和普及。

重力坝;抗滑稳定;Delphi;AutoCAD;二次开发;ActiveX

1 重力坝深层抗滑稳定分析方法

建基面及坝基深层抗滑稳定分析，是重力坝设计中一项非常重要的内容;并且在保证施工、运行的安全性方面起着举足轻重的作用;同时，准确分析坝体的抗滑稳定特性也是重力坝开展其他工作的基本前提［1］。目前重力坝坝基深层抗滑稳定分析方法大致分为3种:刚体极限平衡法、有限单元法和地质力学模型法［2-4］。

随着计算力学和计算机软件及硬件的发展，利用有限单元法可以比较准确地计算出坝体及其坝基的应变、应力和塑性区。依靠计算机，有限单元法能够精确建立模型，并考虑各种材料的非均质、各向异性等特性，因而能较准确地计算接触面以及深层组合滑动面上的阻滑力和滑动力［5-7］。然而仅仅计算出阻滑力和滑动力还不能判断大坝的抗滑稳定性状，需要找到一种计算方法，以较为准确地反映大坝的抗滑稳定特性。本文首先利用有限元法准确计算出重力坝深层滑动面上的阻滑力和滑动力，然后结合分项系数极限状态设计方法，对每一个剖分单元运用抗剪断强度公式，定义一种新的抗滑稳定安全系数Ks的计算公式，得到了一种计算混凝土重力坝深层组合滑动的改进分项系数法。并采用Delphi语言，结合AutoCAD二次开发技术对该方法实现了程序化。该方法充分利用了有限元能够计算出精确应力场和考虑坝基及坝体的变形影响的优势，又考虑了可靠度设计方法，所以得出的Ks能够较好地反映重力坝的抗滑稳定特性，并且本方法沿用了设计规范中的Ks表示方法，易于计算和普及。

2 混凝土重力坝深层组合滑动的改进分项系数法

本方法的基本思想［8］是:引入广义力的概念，将传统的作用效应函数S（·）和结构及构件抗力函数R（·）广义化，结合混凝土重力坝的实际地质情况，分析深层滑移模式，从而决定这种广义力以及该广义力对应的作用效应函数和抗力函数。

下面以平面状2滑面和3滑面（如图1所示）为例，简述混凝土重力坝深层组合滑动的改进分项系数法。根据混凝土承载能力极限状态的基本组合，设计表达式为

式中:γ0为结构重要性系数;ψ为设计状况系数; S（·）为作用效应函数;R（·）为结构及构件抗力函数;Gk为永久作用标准值;γG为永久作用分项系数;Qk为可变作用标准值;γQ为可变作用分项系数;ak为几何参数标准值;γd1为承载能力极限状态基本组合的结构系数;fk为材料性能标准值;γm为材料性能分项系数。

图1 平面状滑面图Fig.1Flat slip surfaces

对于图1中所示的平面状2滑动面和3滑动面，由于各个单元及滑动面上的滑动力和阻滑力不可以直接加减，所以必须将式（1）中的S（·），R（·）进行广义化。事实上，当大坝及坝基沿着组合滑动面出现抗滑失稳的时候，滑移体（由坝体及滑动面以上的岩体组成）受到边界的约束，将绕着滑移体瞬时滑移中心做转动。在此，将阻滑力和抗滑力广义为阻滑力偶和抗滑力偶。

具体计算方法及步骤:对于其中的第1个滑动面，首先根据有限单元法的应力计算成果，利用式（2）和式（3）将滑动面相关单元的应力合成为法向力N1和切向力Q1，同理第2个滑动面为法向力N2和切向力Q2;除了确定力的大小外，还需确定力的作用点。将滑动面上的各单元法向力看成是作用在该滑动面上的平行力系，利用求解平行力系中心的方法可以求得作用点O1和O2的位置。详细求解方法此处不再赘述，见《理论力学》［9］。

过作用点O1，O2画力N1，N2的延长线交于点O。

利用结构计算成果——应力场，得到下面的法向力和切向力的计算公式:

式（2）和式（3）中:n为组合滑动面上划分的单元个数;Ai为第i个单元上的滑动面的面积（二维情况下为长度）;σi为组合滑动面上单元i沿滑动面法向的正应力;τi为组合滑动面上单元i沿滑动面切向的切应力;f'为组合滑动面上材料之间的摩擦系数; c为组合滑动面上材料之间的黏聚力;其余符号意义同前。选取点O作为瞬时转动中心，对中心取矩后，得到广义的S（·），R（·）计算公式:

式（4）、式（5）中:A1，A2分别为组合滑动面中第1、第2滑动面的面积;l1，l2分别为组合滑动面中第1、第2滑动面到交点O的距离;其余符号意义同上。

如图1（b）所示，对于3滑面上的每一个滑动面，按照2滑面的原理得到每个面上的法向力和切向力及其作用点分别为:①N1，Q1，O1;②N2，Q2，O2;③N3，Q3，O3。和2滑面不同的是，经过N1，N2和N3的3条直线正常不会交于一点，所以不能确定瞬心。但是，根据《理论力学》中的知识，可以将任意2条直线的交点作为名义瞬心，将另一边的法向力向该点转化为1个力矩和1个力。然后分别对名义瞬心取矩，方法同2滑面。具体的广义S（·），R（·）计算公式为:

式（6）、式（7）中:d2为瞬时转动中心O点到第2个滑动面上法向合力N2作用线的距离，当N2对O点的矩为逆时针时取正，反之取负。其余符合意义同上。

因此，根据以上对S（·），R（·）的计算，定义混凝土重力坝组合滑面抗滑稳定安全系数Ks，具体定义为

对于3滑面的情况，因为其中可以算出3个不同的Ks，在这里取3个当中最小的作为3滑面的抗滑稳定安全系数。

比较式（1）可得混凝土重力坝深层抗滑稳定分析中的组合滑动面的失稳判据为

除此以外，对于圆弧状单滑面，平面状单滑面和多滑面等组合滑动面，具体计算方法和平面状2滑面或3滑面类似，这里不再赘述。

3 基于ActiveX Automation技术对AutoCAD进行二次开发

AutoCAD作为一种高度开放结构的CAD平台软件，它提供了强大的二次开发环境［10］。其中ActiveX技术是一种完全面向对象的技术，所以许多面向对象化编程的语言和应用程序，可以通过ActiveX与AutoCAD进行通信，并操纵AutoCAD的许多功能［11］。

Delphi利用ActiveX与AutoCAD进行通信，达到操纵AutoCAD的目的。具体做法为:AutoCAD支持Activex技术，其类型库文件中定义了可以被客户程序访问和操纵的对象类型以及这些对象的属性和方法［12］。类型库文件的文件名为ACAD.TLB。通过Delphi的Project菜单中的Import Type Library工具，可以导入类型库文件，将文件转换成Delphi的单元文件格式，这时在Delphi的编辑窗口中会增加标记为AutoCAD2TLB的一页，随后在Delphi的Code Explorer窗口中，用户可以方便地浏览Auto-CAD的ActiveX对象及这些对象的属性和方法，并可以了解调用这些属性和方法时可用的参数及形式。在开发应用程序时，要应用本单元即AutoCAD_ TLB单元［13-14］。

4 计算软件模块研发的关键技术

首先，根据待处理的问题，必须了解以下情况:①在编程之前有哪些已知条件，给了哪些原始资料和数据;②通过编程，要求本程序利用这些资料达到什么目的。

根据以上提出的问题，对于①，已知条件就是混凝土重力坝的网格剖分形式及其应力计算成果;对于问题②，目的就是按照混凝土重力坝深层组合滑动的改进分项系数法计算出给定的组合滑动面的抗滑稳定安全系数Ks。

本程序用户界面如图2。

本程序的总体思路是:

第1步:利用给出的有限元剖分形式，在有限元软件的计算成果中导出每个单元的节点坐标并使之存储在sql数据库中的数据表中，其中有2张表:一张为unit表，用于储存单元信息，包括单元类型、单元号及其包含的节点号;另外一张为coordinate表，用于储存各节点的节点号和三维坐标值。除此之外，程序还设计了unitstress，coslidsurf，relateunit数据表，分别存储单元应力，组合滑动面信息和相关单元信息，各数据表的格式如图3所示。

图2 程序设计图Fig.2Diagram of program design

图3 数据格表Fig.3Formats of data sheets

第2步:利用Delphi调用AutoCAD，然后在CAD图纸上先绘出组合滑动面（这一步骤由使用者完成，这也是人机交互的体现）并用相关技术将其滑动面的参数储存在coslidsurf数据表中，为以后寻找相关单元做准备;之后通过访问sql数据库中的unit表和coordinate表，读出节点坐标，然后利用这些坐标绘出直线，进而绘出有限元剖分网格。

调用AutoCAD的部分代码如下:

从unit表中读取单元号和其中的节点，然后根据节点在coordinate表中查询相应节点号的坐标，作为绘制直线的坐标，从而逐个绘出有限元剖分单元。

第3步:这一步是整个软件的核心，即找出有限元网格中和组合滑动面相交的网格，本文称之为相关单元。当找出所有和组合滑动面相交的相关单元之后，就可以通过这些单元信息，从有限元计算软件的所计算的应力场导出和组合滑动面相关的单元应力信息，最后读取这些应力，通过弹性力学知识，换算成滑动面上的正应力和切应力。然后和单元面积相乘，转化为我们想要知道的滑动力和抗滑力。最后通过改进的组合系数法，求出抗滑稳定安全系数，判断重力坝深层抗滑是否稳定。

寻找相关单元的核心思想是:通过AutoCAD给出的对象模型树中提供的IntersectWith方法，遍历每个单元的每条边，求出组合滑动面和单元的边的交点，确定个数，如果个数不为零，则确定该单元就是相关单元，并且记下滑动面和该单元的交点，求出2交点之间的距离d和滑动面的方向余弦l，m，为之后求解抗滑力和下滑力做准备，具体代码这里不再赘述。

在找出相关单元并且把相关单元的信息，包括单元号、直线号、滑动面和相关单元的交点读入relateunit数据表之后，通过delphi的adoquery控件根据单元号在数据unitstress中查询该单元的应力，即σx，σy，τxy。这里所用单元的应力为单元形心处的应力。

然后根据弹性力学中的知识求出组合滑动面上的应力，即:

式（10）、式（11）中:l，m，σ，τ分别表示滑动面的方向余弦和斜面上的正应力和切应力;然后代入上面式子算出Ks，再利用式（9）判断重力坝的稳定情况。

最后，需要强调的是:本文抗滑稳定安全系数的计算精度与有限元应力结果密切相关，而有限元的应力往往是存在误差的，所以滑动面上的单元应细分，以保证应力求解精度。

5 算例分析

某混凝土重力坝，最大坝高97 m，坝顶全长386.5 m，共由21个坝段组成。其中，5#坝段偏于右岸，是实体重力坝，坝高为60 m，顶宽和底宽分别为6 m和45 m，下游坝坡为1∶0.7。坝基内有倾向下游的与水平面呈18°夹角的软弱结构面，下游抗力体内有倾向上游角度为10°的缓倾角软弱结构面。上述软弱结构面共同形成连续滑动面（见图4），必须验算坝基深层抗滑稳定。

图4 有限元剖分网格及滑动面示意图Fig.4Finite element meshes and the slip surface

应用本文所述的改进分项系数法验算该重力坝深层抗滑稳定。首先建立该坝段的有限元分析数值模型。为了简化实例计算，将该坝段视为平面应变问题，利用平面三角形和四边形等参单元模拟。考虑该坝段及其基础条件的实际情况，选定其有限元计算范围为:上游约1倍坝高，取为60 m;下游约2倍坝高，取为120 m;坝基以下取为60 m。根据上述分析结果，5#坝段有限元网格剖分结果如图4所示，该图为利用有限元软件导出的节点坐标信息，利用本系统在CAD中绘制的网格示意图。该模型共有1 252个单元，1 202个节点。模型中，上游水位为58 m，计算结果见图5。

完成有限元计算之后，将本程序需要的数据按照一定的格式导出，通过相关程序处理之后导入到数据库相应的数据表中（库中各表信息在第4节中已做详细介绍），作为抗滑稳定的已知信息。所需要的数据包括:各节点的坐标信息，各单元的信息以及各单元沿x方向和y方向的应力。

图5 5#坝段最小主应力计算结果Fig.5Calculation result of minimum principal stress for dam segment 5#

运行本程序，完善各计算参数。点击“开始计算”按钮，按系统提示进行操作，最终计算得Ks= 1.24＞1，根据公式（9）判断可得坝体处于稳定状态;而使用传统的刚体极限平衡法可得该坝段深层抗滑稳定安全系数Ks=3.1，同样满足规范规定的相应工况下的安全系数要求。虽然上述采用不同的分析方法，但是所得的结论是一样的。因此，就本算例而言，使用本文提出的改进的分项系数法对混凝土重力坝的深层抗滑稳定性进行分析是可行的，其结果是可信的。但是，该方法仍然需要经过大量的工程实例的验证，以确定其可行性和结果可靠性。

6 结语

（1）随着有限元软件在水利工程中的广泛应用，大量的后续处理工作需要各种软件相互结合，从而要求各软件之间的数据通信。运用ActiveX技术就可以解决数据通信问题。本文利用该技术对AutoCAD进行二次开发，与重力坝深层抗滑稳定分析相结合，收到了很好的效果。结合有限元的计算成果和CAD的强大的绘图功能，从而使水工计算程序化、快速化、自动化、可视化。

（2）基于ActiveX技术，可以对支持该技术的软件，如word，excel等，利用结构化面向对象的编程语言，并结合具体专业问题，进行软件的二次开发，从而对有限元计算结果的后续分析提供一种新的、快捷的、自动化的、实时的处理方法。

（3）结合本专业的具体问题对AutoCAD等支持ActiveX技术的软件进行二次开发的方法，将会给水利工程的设计、分析、监测带来极大的方便，在这方面还有待进一步开发与研究。

［1］李旭.重力坝抗滑稳定问题分析［J］.四川水利，2007，（6）:30-33.（LI Xu.Anti-sliding Stability of Gravity Dam［J］.Sichuan Water Resources，2007，（6）: 30-33.（in Chinese））

［2］苏怀智，胡江，温志萍，等.基于时变风险率的大坝使用寿命评估模型［J］.中国科学E辑:技术科学，2009，52（7）:1966-1973.（SU Huai-zhi，HU Jiang，WEN Zhi-ping，et al.Evaluation Model for Service Life of Dam Based on Time-Varying Risk Probability［J］.Science in China Series E:Technological Sciences，2009，52（7）:1966-1973.（in Chinese））

［3］苏怀智，温志萍，吴中如.智能推理机在大坝健康早期预警中的应用研究［J］.水资源管理，2011，25（6）: 1545-1563.（SU Huai-zhi，WEN Zhi-ping，WU Zhongru.Study on An Intelligent Inference Engine in Early-Warning System of Dam Health［J］.Water Resources Management，2011，25（6）:1545-1563.（in Chinese））

［4］段庆伟，陈祖煜，王杰，等.重力坝抗滑稳定的强度折减法探讨及应用［J］.岩石力学与工程学报，2007，26（增2）:4510-4516.（DUAN Qing-wei，CHEN Zuyu，WANG Jie，et al.Discussion on Strength Reduction Method for Gravity Dam Stability Analysis Against Sliding and Its Application［J］.Chinese Journal of Rock Mechanics and Engineering，2007，26（Sup.2）:4510-4516.（in Chinese））

［5］张军涛，姜礼斌，唐春风.基于有限元法的混凝土重力坝抗滑稳定分析［J］.四川建筑科学研究，2010，36 （6）:132-134.（ZHANG Jun-tao，JIANG Li-bin，TANG Chun-feng.Using FEM Analysis on Slide of Concrete Gravity Dam［J］.Sichuan Building Science，2010，36 （6）:132-134.（in Chinese））

［6］张彬洪，姚激，李泽，等.基于有限元法的混凝土重力坝抗滑稳定分析方法［J］.水利与建筑工程学报，2010，8（1）:62-64.（ZHANG Bin-hong，YAO Ji，LI Ze，et al.Analysis on Stability Against Sliding of Concrete Gravity Dam Based on Finite Element Method［J］. Journal of Water Resources and Architectural Engineering，2010，8（1）:62-64.（in Chinese））

［7］蒋春艳，常晓林，周伟.用于重力坝抗滑稳定分析的分项系数有限元方法［J］.水力发电学报，2006，25 （2）:16-20.（JIANG Chun-yan，CHANG Xiao-lin，ZHOU Wei.Partial Coefficient FEM for Gravity Dam Anti-sliding Stability Analysis［J］.Journal of Hydroelectric Engineering，2006，25（2）:16-20.（in Chinese））

［8］汪亚超.混凝土重力坝抗滑稳定多维时空分析方法研究［D］.南京:河海大学，2011.（WANG Ya-chao. Multi-dimensional Spatio-Temporal Analysis Method for the Anti-sliding Stability of Concrete Gravity Dam［D］. Nanjing:Hohai University，2011.（in Chinese））

［9］陈世民.理论力学简明教程［M］.北京:高等教育出版社，2001.（CHEN Shi-min.Brief Tutorial of Theoretical Mechanics［M］.Beijing:Higher Education Press， 2001.（in Chinese））

［10］戴春来，陈功，陆信.AutoCAD Automation技术的应用［J］.工程设计CAD与智能建筑，1999，（2）:56-58.（DAI Chun-lai，CHEN Gong，LU Xin.Application of AutoCAD Automation［J］.Computer Aided Design and Intelligent Building，1999，（2）:56-58.（in Chinese））

［11］周强.基于ActiveX对象模型树的参数化绘图技术的研究与应用［J］.南京化工大学学报，2001，23（5）: 83-86.（ZHOU Qiang.Research and Application of Graph Parametrizing Technique Based on ActiveX Model Tree［J］.Journal of Nanjing University of Chemical Technology，2001，23（5）:83-86.（in Chinese））

［12］王玉.用VBA开发AutoCAD 2000应用程序［M］.北京:人民邮电出版社，1999.（WANG Yu.Development of AutoCAD 2000 Applications by Using VBA［M］. Beijing:Posts＆Telecom Press，1999.（in Chinese））

［13］赵峰，毕苏萍，周振红.Delphi通过双接口控制AutoCAD绘图［J］.郑州大学学报（工学版），2007，28 （4）:66-69.（ZHAO Feng，BI Su-ping，ZHOU Zhenhong.Delphi Driving AutoCAD to Draw by Way of Dual Interfaces［J］.Journal of Zhengzhou University（Engineering Science），2007，28（4）:66-69.（in Chinese））

［14］李兵，赵锡英.Delphi在CAD系统中的应用［J］.兰州工业高等专科学校学报，2004，11（4）:16-19.（LI Bing，ZHAO Xi-ying.The Application of Delphi in CAD System［J］.Journal of Lanzhou Polytechnic College，2004，11（4）:16-19.（in Chinese））

（编辑:姜小兰）

Software Module to Assess the Deep Sliding Stability of Concrete Gravity Dam

CUI Shu-sheng1，2，SU Huai-zhi1，3，ZHANG Shuai1，2，GU Hao-qing1，2
（1.State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering，Hohai University，Nanjing 210098，China;2.National Engineering Research Center of Water Resources Efficient Utilization and Engineering Safety，Hohai University，Nanjing210098，China;3.College of Water Conservancy and Hydropower，Hohai University，Nanjing210098，China）

To precisely reflect the anti-sliding stability of dam，a visualized software module was researched and developed to assess the deep-seated sliding stability of concrete gravity dam.Firstly，the stress distribution on the deep-seated sliding plane of the dam was numerically calculated.On this basis，by using the design method of the limit state of partial coefficient，the formula of shear strength was applied to each FEM element to define a new formula for the safety coefficient，thus an improved method of partial coefficient for the stability assessment was developed.In subsequence，the method was programmed by programming tool Delphi in association with AutoCAD. Since the advantage of FEM in calculating stress filed and considering dam foundation and body deformation is reflected in the method and the software module，and owing to the reliability degree taken into account as well，the obtained safety coefficient could well reflect the anti-sliding stability.Moreover，Ksexpression in design specification is employed to make the method easy to use and popularize.

gravity dam;anti-sliding stability;Delphi;AutoCAD;secondary development;ActiveX

TV642.3;TP311.56

A

1001-5485（2013）04-0085-06

10.3969/j.issn.1001-5485.2013.04.018

2013，30（04）:85-90

2012-08-27;

2012-10-17

国家自然科学基金（51179066，51139001）;新世纪优秀人才支持计划（NCET-10-0359）;国家重点实验室专项经费资助项目（2009586912，2009586012）;中央高校基本科研业务费专项资金资助（2012B06614）;江苏省第四期“333工程”培养资金资助项目（BRA2011179）;江苏高校优势学科建设工程资助项目（水利工程）（YS11001）

崔书生（1987-），男，安徽六安人，硕士研究生，研究方向为水工结构工程安全监控理论和方法，（电话）13770516712（电子信箱） shushengjiayou@163.com。