坝基面受力状态的有限元计算误差影响分析

程帅+李守义+司政+杜占科+李萌+郝晓飞

摘要：有限元数值计算在水利工程中运用广泛，其计算精度和计算误差是需要考虑和控制的重要问题，对其进行研究可为有限元法在水利设计和复核中提供技术支撑。以重力坝为研究载体，根据坝体荷载合力值与有限元计算的坝基面内力值（坝基面应力积分值）的平衡条件，分析有限元应力计算误差的影响因素，提出减小误差的措施。研究结果表明：在同一受力状态下，提取应力的对象不同计算误差不同，需针对所研究的问题及误差控制要求合理选择相应提取应力的对象；将紧接坝体的薄层地基保持与坝体相同的模型尺寸和约束条件，可充分减小坝基面有限元应力计算误差；合理选择模型网格尺寸、地基模拟范围、应力提取路径的条数、提取点间距，可在不消耗过量的计算资源和时间下满足精度要求。研究成果可为控制类似结构的有限元计算误差提供参考。

关键词：水工结构工程；有限元法；误差分析；坝基面；受力状态

中图分类号：TV314文献标志码：A文章编号：1672-1683（2017）01-0179-07

Abstract：Finite element numerical calculation is widely used in hydraulic engineering.Its accuracy and error is one of the important problems that need to be considered and controlled.Studies on this issue can provide technical support for the FEM in the design and check of hydraulic engineering.Gravity dam was used as the object of research.The factors influencing the error in finite element stress calculation were analyzed by comparing the resultant load value of the gravity dam with the internal force value of the dam base surface from finite element calculation （integral value of stress of dam base surface）.The results showed that extracting the stress from different objects results in different finite element calculation error，even under the same stress state.It is necessary to select the proper object to extract the stress from according to the research goal and error control requirements.The error can be greatly reduced by making the thin-layer foundation adjacent to the dam the same size and in the same constraint condition as the dam body.Proper determination on mesh size，simulation range of foundation，number of paths to extract the stress，and spacing of extracting points can meet the required precision without overly consuming resources or time.The results can provide reference for controlling the finite element calculation error in similar projects.

Key words：hydraulic structure engineering；finite element method；error analysis；dam base surface；forced state

隨着计算机技术及数值分析方法的迅速发展，有限元数值计算在各个领域得到了越来越广泛的应用[1-2]。近几年，有限元法在水利工程中正发挥着重要作用，特别是大中型水利工程中有限元的应用已经相当广泛[3-5]。目前，对于复杂水利工程中的有限元计算结果，在没有试验数据和标准答案的情况下，如何判断其正确性，如何用所掌握的信息对设计进行完善和修正，如何控制其计算精度和计算误差，这些问题没有得到很好的解决。这也是为何至今利用有限元软件进行有限元分析的水工建筑物或水力模型仅作为水利工程设计参考和复核，并未给出设计或检验标准的原因之一。

影响有限元计算精度的因素很多，Jin[4]，Mirzabozorg[5]，Kennedy[6]，Natarajan[7]等都针对提高有限元计算精度开展了很多研究。国内已有学者对其中的部分因素进行了探讨。虞皓影[8]、张国新[9]等在利用有限元分析结构稳定性中指出边界计算范围、应力积分方法等因素对计算结果影响较为敏感，并各自提出了减小有限元计算误差的方法。邹超英等[10]通过对拱坝进行应力分析，研究了有限元单元类型、自由度数量、单元阶次等因素与能量误差的关系。网格密度也是一个影响计算精度的因素，李涛等[11]、章春亮等[12]等指出，有限元求解模型中所有关键性区域的网格划分水平决定整个有限元模型的计算结果精度，理论上单元数量的增多，结果的精确度就会提高，并收敛于真实解。单元类型及网格划分方法对有限元计算精度也有一定影响，姬贺炯等[13]、Cavin等[14]、张丽媛等[15]分别对这些因素进行了研究，为相应结构的有限元计算误差提供了理论依据。然而，对于大体积水工结构的模拟，网格划分过密，会消耗大量的资源和时间，需确定相应网格尺寸来兼顾计算效率和精度。此外，运用有限元分析结构受力状态、稳定性等问题时，通常需要选择整体结构中的不同对象来提取计算截面的应力，选择不同的提取应力对象，同一受力位置的有限元计算误差不同，而在既定的应力提取对象下应力路径的范围及疏密程度也会影响其计算精度。目前，在大体积水工结构的有限元计算中，针对上述几点影响因素及误差影响分析方面的研究甚少。

坝基面作为坝体的薄弱控制面，该部位的受力状态对分析坝体稳定性等问题有很大影响。本文以重力坝为研究载体，根据力学平衡条件得出坝基面有限元计算结果的误差，对同一位置不同提取应力对象下有限元计算值的偏差及其原因进行了分析，提出了减小误差的措施，考察了网格尺寸大小、地基模拟范围、应力路径条数、应力点间距等因素对坝基面有限元计算误差的影响。研究结果可为分析导墙、闸墩、厂房、挡土墙等类似工程结构的有限元计算误差提供依据。

1 研究资料

1.1 研究对象概况

本文研究主要依托某水利枢纽工程中一重力坝坝段开展。该工程属Ⅱ等大（2）型工程，拦河坝为混凝土重力坝，最大坝高75.97 m，坝顶宽度7.0 m，上游坝面铅直，下游坝坡1∶0.8，最大坝底宽度60.74 m。正常蓄水位450 m，死水位440 m，下游正常尾水位405 m，坝前淤沙高程389.3 m。

1.2 荷载施加

本文在进行有限元计算时，考虑的荷载有：自重、上下游静水压力、扬压力、浪压力、泥沙压力。

1.3 材料参数

不同材料的计算参数取值见表1。

2.1 有限元模型

本文采用ANSYS有限元软件对某一重力坝坝段进行分析，坝段宽18 m，整体有限元计算模型包括坝体与基础岩体。基岩的模拟范围取矩形区域，基岩边界值：深度方向均取152 m，约为2倍坝高，左、右岸方向与坝段同宽，上下游方向取2倍坝高（在研究地基范围影响时分类考虑）。模型整体坐标系的原点设在坝踵右侧底部，沿水流方向指向下游为X轴正方向；沿高度方向铅直向上为Y轴正方向；垂直水流方向指向右岸为Z轴正方向。整体模型网格剖分基本采用8结点六面体实体单元，坝基紧接坝体10 m范围内网格尺寸与坝体相同，坝基其它区域网格尺寸较大，中间10 m厚度的槽型区域通过四面体单元进行过渡。整体有限元计算模型见图1。

2.2 计算假定

根据重力坝、厂房等类似结构的特点及其相关研究成果[16-18]，本文在三维有限元静力计算中，做如下假定。

（1）坝段间设置有横缝，因此各方案下坝段独立承受荷载，坝段间无相互作用。

（2）充分考虑坝体的材料分区，在静力计算中假定坝体混凝土是线弹性的，地基为均匀无质量弹性体，二者皆满足弹性力学基本假定。

（3）由于横缝作用，在计算模型中坝体两侧无约束，地基岩体上下游面及左右侧面考虑法向约束，基岩底面考虑三向约束（在研究边界约束影响时分类考虑）。

2.3 研究方法

分析截面受力状态，需对应力进行曲面积分以得到建基面截面内力[19-21]。由以上计算模型和资料，施加荷载，得到有限元应力计算结果。选取上部坝体结构或下部地基结构作为应力提取对象，沿着坝轴线方向取N条路径，并选定每条路径的提取点间距，提取建基面上的法向应力和切向应力，按式（1）沿建基面进行积分后得到法向和切向力的合力。

式中：σy，τxy为提取的各点法向和切向应力；Fy，Fs为法向和切向力的合力值。

根据力学平衡原理计算出外荷载合力值与有限元计算得出的内力值作对比，按式（2）分析其相对误差。

式中：F0为外荷载合力值；Fi为有限元计算的内力值。

采取一系列方案得到有限元计算结果，根据各方案相对误差分析坝基面受力状态的有限元计算误差影响因素，研究减小误差的措施。

3 研究内容与结果分析

3.1 提取应力对象的影响

分别选取上部坝体结构和下部地基结构为应力提取对象，提取坝基面上的法向和切向应力，积分得合力值与外荷载合力值对比，计算其相对误差见表2。

由表2可知，有限元计算结果中，提取应力的对象不同，法向力和切向力都有所变化，且误差不同，法向力误差对于应力提取对象的敏感度小于切向力误差。选择下部地基结构作为应力提取对象计算误差偏大，法向力误差超过5%，切向力误差甚至达到12.3%。在实际分析中，需合理选择提取应力的对象，才能保证计算结果的精确性。由于坝基面部位对分析坝体稳定性等问题极为关键，且某些问题中必须以下部地基结构作为提取应力的对象，以下分析选择下部地基作为提力对象时误差过大的原因。

由地基与坝体交接面的特点可知，这种上下偏差较大的现象可能是由于材料屬性变化、模型尺寸突变、边界约束条件的不同等因素造成的。现针对以上几种可能的影响因素做以下六种方案进行研究：方案①：原始计算模型，材料、尺寸、边界约束都不变；方案②：将地基的材料属性变为与坝体相一致，其它条件不变；方案③：将地基范围的X方向长度取为与坝体相等，其它条件不变；方案④：地基范围不变，在紧接坝体部位一个单元厚度处地基与坝体X方向长度相等，且约束不变；方案⑤：在地基的原约束基础上，将建基面向下一个单元厚度的地基约束去掉，使其与坝体的约束条件相一致，地基其他边界约束不变；方案⑥：地基范围不变，在紧接坝体部位一个单元厚度处地基与坝体X方向长度相等，且无约束，其它条件不变。各方案的模型、材料及边界约束示意图见图2。

根据图2中各种方案，建立相应的有限元计算模型，施加荷载，分别以上部坝体结构和下部地基结构作为提取应力的对象，坝基面上的法向和切向力的相对误差计算结果见图3。

由图3可知，选择下部地基结构误差较大，其中法向力误差在2.8%～6.5%之间，切向力误差在0.3%～27.9%之间；选取上部坝体结构计算相对较精确，其中法向力误差在1.7%～3.5%之间，切向力误差在0.7%～1.5%之间。切向力误差变幅较法向力更大，其中下部切向力误差对方案的选取最为敏感，变幅最大。

由方案①和②可知，虽然将坝体与地基的材料属性取为一致，但其上下切向应力值相差依然很大，下部切向力误差都在10%以上，故造成这种误差的原因并非材料属性的变化所致。由方案③和④可知，两种不同地基范围下，在坝基面处模型尺寸均无突变，但上下切向应力相差依然很大，下部切向力误差达到28%，可知造成这种误差的原因并非模型尺寸突变所致。方案⑤和⑥验证是否为边界约束条件的影响，方案⑤中解除了建基面向下一个单元厚度的地基约束，误差值并未减小，下部切向应力依然有12%，各项误差与方案①相近。方案⑥中紧接坝体部位一个单元厚度处地基与坝体X方向长度相等，且无约束，使此处一小部分地基的约束条件和模型尺寸皆与坝体相一致，这种情况下，地基面上、下部的法向力和切向力误差都达到最小，上部和下部提取的力几乎相等，法向误差减小至1.7%，切向误差减小至0.3%。

上述六种方案充分说明：造成坝基上部与下部的所提取力的较大偏差的原因，不是材料属性变化所致，而是模型尺寸突变和边界约束条件变化的双重原因所致。故为充分减小有限元计算误差，建议将紧接坝体很小厚度范围（一个单元尺寸）内的薄层地基保持与坝体同宽且不做任何约束，即保持与坝体具有相同的模型尺寸和约束条件。

3.2 网格尺寸的影响

为研究坝体网格尺寸与地基网格尺寸对计算精度的影响，现设置A、B两组方案见表3。其中A组方案中地基网格尺寸固定为10 m，研究坝体网格尺寸由10 m减小至1 m时的计算误差；B组方案中坝体网格尺寸固定为1 m，研究地基网格尺寸由10 m减小至1 m时的计算误差。计算结果见图4。

由图4可知，网格尺寸的大小会影响计算精度。坝体网格尺寸的变化对计算误差的影响更为敏感，结果与文献[11]中关键区域的网格划分水平决定整个有限元计算结果精度的观点相吻合。当坝体与地基的网格尺寸同为10 m时，法向力误差达到7.45%，切向力误差达到2.61%，随着坝体网格尺寸的减小，计算误差也随之减小，计算精度越精确。当网格尺寸减小到一定程度以后，对计算结果准确性的提高就很小，且会导致计算耗时大幅度增加。当坝体网格尺寸为2 m时，法向误差降至2.41%，切向误差降至0.63%。地基的网格尺寸对计算结果影响甚微，当坝体网格尺寸固定为1 m时，地基网格尺寸由10 m减小至1 m，计算误差都维持在一个相对稳定的水平。综合考虑计算机耗时和计算精度的影响，在水利工程中大体积结构模拟中，地基网格尺寸取为10 m，结构网格尺寸取为2 m，即可满足分析要求。

3.3 地基范围的影响

为分析地基模拟范围对计算精度的影响，分C、D两组方案进行研究，各方案取值如表4所示（表中H为坝高）。其中C组方案中地基深度方向（Y向）保持2倍坝高不变，研究上下游方向（X向）地基范围由坝踵和坝趾分别向上、下游延伸不同长度时坝基面合力的计算误差；D组方案中X向地基范围分别向上、下游延伸2倍坝高不变，研究Y向地基不同深度范围下的坝基面合力值的计算误差。计算结果见图5。

由图5可知：地基范围的不同会导致计算精度的不同。X方向地基范围的变化对计算误差的影响更为敏感。Y方向地基深度保持2倍坝高不变，随着X方向地基范围的增大，计算误差整体上呈减小趋势，法向误差在X方向地基范围增大时，也会出现略微的增大现象，当X方向地基范围由坝踵和坝趾分别向上游和下游延伸2倍坝高及以上时，法向误差会控制在2.8%以下，切向误差会控制在0.85%以下。地基Y方向的范围对计算结果的误差影响较小，当X方向地基范围由坝踵和坝趾分别向上游和下游延伸2倍坝高不变时，地基深度范围由0.2倍坝高增加至1倍坝高时，随深度的增加误差略有减小，当深度大于1倍坝高时，计算误差虽有所偏差，但浮动较小，法向误差在2.8%左右，切向误差在0.85%左右。

考虑到地基范围过大会造成模型单元与节点数目的成倍增加，需要更多的计算资源和时间。故建议在类似工程的模拟中，X方向地基范围由坝踵和坝趾分别向上游和下游延伸2倍坝高，Y方向地基深度取为2倍坝高。

3.4 应力路径条数及应力点间距的影响

以上是针对有限元应力计算的误差影响分析。在确定的应力计算结果下，坝基面的合力计算值也会有不同的精度，这是由于计算合力时应力积分路径条数与应力点间距的选取不同而造成的，现通过E、F两组方案分析这两个因素对坝基面合力值计算的误差影响，各方案设置见表5（表中B为坝段宽度，L为路径总长度）。其中，E组方案在相同的应力点间距下（0.01L），沿着坝轴线方向取不同的应力路径条数，分析应力路径条数对计算精度的影响；F组方案固定路径条数（B/2），分析不同应力点间距下计算误差的影响。各方案计算结果见图6。

由图6可知，应力提取路径的条数，以及每条路径提取点的间距，对计算精度有一定的影响。同一提取点间距，随着路径条数的增多，计算误差会减小，当条数大于B/2时，法向力误差可控制在2.6%以内，切向力误差可控制在1.0%以内，且继续增多路径精度的提高效果不明显。在同一提取路径条数（B/2）下，当提取应力点间距过大时，为负偏差（标准值小于計算值），间距过小时，为正偏差（标准值大于计算值）。随着应力提取点间距的减小，计算误差会由负到正逐渐增大，最后会处于一个稳定状态，当间距小于0.01 L时，法向力误差可控制在2.0%以内，切向力误差可控制在0.9%以内。

由以上结果可知，应力路径条数取B/2即可满足分析要求。综合考虑计算精度和计算工作量，提取应力点间距取为0.006 25L～0.02L为宜。

4 结论

（1） 在不同提取应力的对象下，坝基面受力的有限元计算值不同，且计算误差不同。选择上部坝体结构作为研究对象计算误差较小，地基面上的法向应力误差对于提力对象的敏感度小于切向应力误差。在类似工程中，需针对所研究的结构及误差控制要求合理选择相应提取应力的对象。

（2） 在有限元计算模型中，将紧接坝体很小厚度范围（一个单元尺寸）内的薄层地基保持与坝体相同的模型尺寸和约束条件，可充分减小有限元计算误差。其中法向误差可控制在2.0%以内，切向误差可控制在0.5%以内。

（3）对类似大体积水工结构进行模拟并分析其受力状态时，综合考虑计算精度、计算资源和时间、工作效率等因素，结构网格尺寸取为2 m，地基网格尺寸取为10 m，地基范围在上下游方向向两侧延伸2倍建筑物高度，深度方向取为2倍建筑物高，即可在不消耗过量的计算资源和时间下满足精度要求；应力路径条数取B/2即可满足分析要求，应力点间距取为0.006 25L～0.02L为宜。

考虑到影响有限元计算精度的因素很多，如模型本身的结构误差、阶次误差，荷载分布的简化，材料本构模型选取等，各个因素是相互关联，共同作用的。此外，亦可通过其他方法评价有限元计算结果精度，如对比试验数据、观察应力等值线图、核算关键节点位移等。这些规律还有待于进一步研究。本文以经典弹性力学为基础，视混凝土和地基材料为均质连续体，通过对比有限元应力计算值与外荷载合力值，论述了部分影响坝基面受力状态计算误差的因素，所得结论适用于水利工程中的大体积混凝土结构，可满足实际工程的分析需求，对类似工程具有一定的指导意义。

参考文献（References）：

[1] 曾攀.工程有限元方法[M].北京：科学出版社，2010.（ZENG Pan.The finite element method in engineering[M].Beijing：Science Press，2010.（in Chinese））

[2] 何本国.ANSYS土木工程应用实例[M].北京：中国水利水电出版社，2011.（HE Ben-guo.Application examples of ANSYS civil engineering[M].Beijing：China Water Power Press，2011.（in Chinese））

[3] 王树仿，李向明.三维有限元静力分析在小型引水式电站厂房中的应用[J].浙江水利科技，2014，42（4）：49-53.（WANG Shu-fang，LI Xiang-ming.A small diversion-type power station plant three-dimensional finite element static analysis[J].Zhejiang Hydrotechnics，2014，42（4）：49-53.（in Chinese））DOI：10.13641/j.cnki.33-1162/tv.2014.04.020

[4] Jin B，Lazarov R，Pasciak J，et al.Error analysis of a finite element method for the space-fractional parabolic equation[J].SIAM Journal on Numerical Analysis，2014，52（5）：2272-2294.DOI：10.1137/13093933X

[5] Mirzabozorg H，Hariri-Ardebili M A，Heshmati M，et al.Structural safety evaluation of Karun III Dam and calibration of its finite element model using instrumentation and site observation[J].Case Studies in Structural Engineering，2014，1：6-12.DOI：10.1016/j.csse.2014.02.001

[6] Kennedy F E，Colin F，Floquet A，et al.Improved techniques for finite element analysis of sliding surface temperatures[J].Developments in Numerical and Experimental Methods Applied to Tribology，2014：138-150.DOI：10.1016/B978-0-408-22164-1.50025-9

[7] Natarajan S，Ooi E T，Chiong I，et al.Convergence and accuracy of displacement based finite element formulations over arbitrary polygons：Laplace interpolants，strain smoothing and scaled boundary polygon formulation[J].Finite Elements in Analysis and Design，2014，85：101-122.DOI：10.1016/j.finel.2014.03.006

[8] 虞皓影，何建平.影響有限元法分析滑坡稳定性的因素研究[J].河南工程学院学报：自然科学版，2009，21（1）：33-35.（YU Hao-ying，HE Jian-ping.Research on factors of element method for slope stability analysis[J].Journal of Henan Institute of Engineering，2009，21（1）：33-35.（in Chinese））

[9] 张国新，刘毅.坝基稳定分析的有限元直接反力法[J].水力发电，2006，32（12）：30-32.（GHANG Guo-xin，LIU Yi.Finite element direct force method in stability analysis of dam foundation [J].Water Power，2006，32（12）：30-32.（in Chinese））

[10] 邹超英，苏志敏，曾海军，等.拱坝有限元分析网格剖分方案研究[J].南水北调与水利科技，2011，9（1）：47-49.（ZOU Chao-ying，SU Zhi-min，ZENG Hai-jun，et al.Mesh generation of finite element method for arch dam[J].South-to-North Water Diversion and Water Science&Technology，2011，9（1）：47-49.（in Chinese））DOI：10.3724/SP.J.1201.2001.01047

[11] 李涛，左正兴，廖日东.结构仿真高精度有限元网格划分方法[J].机械工程学报，2009，45（6）：304-308.（LI Tao，ZUO Zheng-xing，LIAO Ri-dong.Meshing method of high precision FEM in structural simulations[J].Chinese Journal of Mechanical Engineering，2009，45（6）：304-308.（in Chinese））DOI：10.3901/JME.2009.06.304

[12] 章春亮.有限元分析中的单元划分质量与计算精度[J].轻工机械，2002 （1）：27-30.（ZHANG Chun-liang.Meshing quality and calculation precision of finite element analysis[J].Light Industry Machinery，2002 （1）：27-30.（in Chinese））

[13] 姬贺炯，白长青，韩省亮.输流管道动力有限元建模及实验研究[J].应用力学学报，2013，30（3）：422-428.（JI He-jiong，BAI Chang-qing，HAN Sheng-liang.Dynamic finite element modeling and experimental research of the fluid-filled pipeline[J].Chinese Journal of Applied Mechanics，2013，30（3）：422-428.（in Chinese））DOI：10.11776/cjam.30.03.C059

[14] Cavin P，Gravouil A，Lubrecht A A，et al.Efficient FEM calculation with predefined precision through automatic grid refinement[J].Finite Elements in Analysis and Design，2005，41（11）：1043-1055.DOI：10.1016/j.finel.2004.09.004

[15] 张丽媛，胡昱，李庆斌，等.基于平面布尔运算的重力坝3D有限元建模[J].水力发电学报，2012，31（5）：209-215.（ZHANG Li-yuan，HU Yu，LI Qing-bin，et al.Hex-meshing method for gravity dam based on planar Boolean operations[J].Journal of Hydroelectric Engineering，2012，31（5）：209-215.（in Chinese））

[16] 袁自立，马福恒，焦延涛.石漫滩碾压混凝土重力坝温度效应分析[J].南水北调与水利科技，2013（5）：61-64.（YUAN Zi-li，MA Fu-heng，JIAO Yan-tao.Analysis of temperature effects of Shimantan RCC gravity dam[J].South-to-North Water Transfers and Water Science & Technology，2013（5）：61-64.（in Chinese））DOI：10.3724/SP.J.1201.2013.05061

[17] 童偉，张连明，范书立.官地碾压混凝土重力坝的抗震分析[C]中国碾压混凝土筑坝技术.2010.（TONG Wei，ZHANG Lian-ming，FAN Shu-li.Antiseismic analysis for Guandi RCC gravity dam[C].RCC Damming Technology of China，2010.（in Chinese））

[18] 李守义，肖静.河床式水电站厂房横缝止水布置形式研究[J].水力发电学报，2014，33（5）：165-168.（LI Shou-yi，XIAO Jing.Study on waterstop arrangements of transverse joints of hydropower house in river channel[J].Journal of Hydroelectric Engineering，2014，33（5）：165-168.（in Chinese））

[19] 朱伯芳，高季章，陈祖煜，等.拱坝设计与研究[M].北京：中国水利水电出版社，2002.（ZHU Bo-fang，GAO Ji-zhang，CHEN Zu-yu，et al.Design and research of arch dam[M].Beijing：China Water & Power Press，2002.（in Chinese））

[20] 杨强，刘福深，周维垣.基于直接内力法的拱坝建基面等效应力分析[J].水力发电学报，2006，25（1）：19-23.（YANG Qiang，LIU Fu-shen，ZHOU Wei-yuan.Equivalent stress analysis for the base surfaces of arch dams by direct internal force method[J].Journal of Hydroelectric Engineering，2006，25（1）：19-23.（in Chinese））

[21] 覃海艺，马宁.薄壁弯曲结构有限元计算精度研究[J].应用力学学报，2015，32（1）：139-145.（QIN Hai-yi，MA Ning.Finite element calculating precision of thin-wall bending structure[J].Chinese Journal of Applied Mechanics，2015，32（1）：139-145.（in Chinese））DOI：10.11776/cjam.32.01.B109