基于ANSYS的拱坝等效应力计算及图形显示

2011-04-28李业盛赖国伟甘海阔

水力发电 2011年4期

李业盛，赖国伟，甘海阔

（武汉大学水资源与水电工程科学国家重点实验室，湖北武汉 430072）

拱坝是高次超静定结构，坝体应力对于边界条件的影响十分敏感，与纯拱法、拱梁分载法等结构力学方法相比，有限元法可以考虑坝体大孔口、不规则外形 (如溢流堰)、拱坝与基岩的相互作用、复杂基础等因素的影响，计算精度较高。但在弹性范围内，有限元计算成果在坝踵、坝趾等角缘部位存在应力集中现象，且网格愈密应力集中程度愈高，不利于确定应力控制指标。为此我国一些学者提出了 “有限元等效应力”。傅作新等[1]将有限元法所求得的应力合成为截面内力，求出对应的线性化应力作为等效应力，但截面内力是拟合出来的，不能满足内力平衡条件[2]。李同春等[3]将坝体分解为拱系和梁系，根据拱和梁的内力平衡条件求解指定截面上的约束内力，进而求解相应截面上的内力和坝体内任一点的等效应力。杨强等[4]将建基面沿高程依次分成若干段曲面，对每段曲面均采用等效矩形进行近似，分别积分求解各等效矩形截面的内力。朱伯芳[5]认为在计算有限元等效应力时，宜直接进行数值积分，而不宜用二次曲线逼近，因拱坝是偏心受压结构，应力分布是一条有拐点的高次曲线。李守义等[6]经过对比分析认为,基于ANSYS的拱坝等效应力计算结果的应力分布规律与拱梁分载法计算结果的分布规律基本一致，结果可靠，计算精度和效率较高。

尽管ANSYS后处理功能强大，但并不能像有限元计算结果一样将等效应力结果在后处理模块中直接输出为工程上常用的、直观的等值线、云图等形式，一般需借助TECPLOT、MATLAB等第三方后处理软件将计算结果重新整理并输出[6]，过程略显繁琐。本文在获得ANSYS拱坝有限元计算结果数据库的基础上，拟利用ANSYS参数化编程语言（ANSYS Parameter Design Language,APDL）直接沿拱坝拱和梁断面进行数值积分得到截面内力，计算上下游坝面各结点等效应力，同时采用单元转化及结点输出列表修改的方法，在ANSYS后处理模块中一次性完成拱坝上下游面各结点等效应力的计算和结果图形化显示。

1 基于ANSY S的有限元等效应力计算

与弹性壳体理论相似，假设三个主要应力分量（σx，σy，σz）沿坝厚为线性分布[1]，根据弹性有限单元法求得的拱坝应力分量，沿梁拱断面直接进行数值积分，得到梁与拱的内力，即可用材料力学方法计算坝面有限元等效应力，据此按规范即可进行拱坝应力安全评价。但基于ANSYS的有限元等效应力计算还有几个关键问题需要解决。

1.1 基于APDL编制计算程序

基于ANSYS有限元计算软件平台，利用APDL编制开发拱坝的有限元等效应力的计算程序。其主要计算步骤为：①根据拱坝体形尺寸，确定有限元等效应力计算结点的局部坐标系与梁、拱截面几何参数；②由有限元法计算的应力分量，分别沿梁、拱截面积分，计算得到梁、拱截面的内力；③用材料力学方法计算上下游坝面各结点有限元等效应力分量；④最后计算上下游坝面各结点有限元等效应力主应力。

1.2 计算坐标系设置

有限元计算模型整体坐标系采用原点为拱冠梁中面坝顶上游面点在海平面的投影、+X轴指向右岸、+Y轴指向顺河流向、+Z铅直向上的直角坐标系。计算等效应力过程中需先后建立两套局部坐标系。在1.1步骤②沿路径积分有限元应力分量计算截面内力时，需在坝面结点i建立如图1所示局部坐标系Ⅰ，其中x轴平行于拱中心线的切线方向，y轴平行于半径方向，指向上游为正，z轴为铅直方向，向上为正，原点在中心线上。在1.1步骤③计算上下游坝面各结点等效应力时，为适应工程上拉负压正的应用习惯，采用拱梁分载法的坐标系统，需将局部坐标系Ⅰ绕Y轴顺时针旋转180°，局部坐标系Ⅱ如图2所示。在ANSYS后处理中可用RSYS命令激活局部坐标系。

图1 内力计算时的局部坐标系Ⅰ

图2 等效应力计算时的局部坐标系Ⅱ

1.3 沿路径积分

在1.1步骤②中局部坐标系Ⅰ下，积分的路径需设定为沿拱圈径向方向。梁的水平截面在拱中心线上取单位宽度，沿厚度方向对梁的应力及其矩进行积分；单位高度拱圈的径向截面，沿厚度方向对拱应力及其矩进行积分，内力积分公式见文献[7]。

以下坝面为例，当计算下坝面上A点时，不难确定出其沿水平拱圈径向对应的上坝面B点，以A、B两点为控制点形成积分路径，如图3所示。在ANSYS中沿路径积分有多种方法，这里介绍一种简单的思路，用PATH命令定义A、B两点确定的积分路径，通过PDEF命令、PVECT命令分别将有限元计算的应力结果、积分点的位置在局部坐标系I下映射到积分路径上，最后用PCALC命令沿路径积分应力分量及其矩，即得到梁与拱的内力。

图3 路径示意

1.4 夹角φ、η求解

利用梁与拱的内力，即可根据材料力学的基本公式和平衡条件，计算坝体表面的有限元等效应力。如图4示，在下游坝面取出一四面微元体ABCO，有6个应力分量作用在相互垂直的3个平面ACO、ABO和BCO上，同时有水压力p作用在坝面ABC上。

水平面上悬臂梁铅直正应力σz、切向水平剪应力τzx及径向铅直平面上拱的水平正应力σx计算公式见文献[7]，根据力的平衡条件，可求出未知的3个应力分量如下

图4 四面体ABCO的应力示意

式中，φ为径向铅直平面内坝面与铅直线的夹角，即AC与铅直向（Z轴）的夹角；η为水平面内坝面与拱中心线切线的夹角，即BC与切向（X轴）的夹角；p为坝面水压力。

夹角φ、η可根据面ABC的法向向量求得，过程为：在局部坐标系Ⅱ下，设下坝面法向向量为（nx， ny， nz），该法向向量可在 ANSYS计算程序中用面的投影方式求解提取。今在此法向向量的基础上任意作一平面，其平面方程为

式中，d为任意常数。

交线AC的方程为

由此并结合四面微元体ABCO上的力平衡推导过程，可得

同理，可得交线BC的方程为

及

由式（6）、（8）求得的夹角 φ、 η 代入式（1）～（3）中即得 τxy、 τyz、 σy三个应力分量，与已求得的 σz、τzx及σx共6个应力分量，可进一步求得下游坝面上各结点的主应力值。上坝面的各点主应力值可类似求得。

2 等效应力结果图形显示

在建立有限元计算模型过程中，以上下游坝面结点集为基础创建不参与有限元计算的坝面MESH200三维表面虚单元集。在后处理模块中完成有限元等效应力计算后，将上下游坝面MESH200表面虚单元集转化为能显示坝面结点位移分量（ux，uy，uz）的SURF65表面效应单元集，再用DNSOL命令将SURF65表面效应单元集的结点位移分量修改为计算得到的结点等效应力主应力值（s1，s2，s3），即可直接输出类似有限元计算结果一样的上下游面的等效应力等值线图或彩色云图。

3 工程实例

3.1 有限元计算模型

某双曲混凝土拱坝，坝顶高程355.30 m，坝顶厚度1.7 m，坝底厚度5.0 m，最大坝高53.3 m，最小中心角40°，坝顶弧长126 m，厚高比0.094，属于薄拱坝。该工程三维有限元网格划分见图5。整体模型模拟了左拱座所在相对较单薄的转向山体和坝下游河谷实际地形，以及地基F3、F4断层。对坝体，模拟了溢流堰、8条临时施工横缝及坝体混凝土、浆砌石不同材料分区。单元采用三维八结点等参单元，沿坝厚方向设置五层单元，整体网格共48 480个结点，49 719个单元，其中坝体结点11 915个，单元9 455个。

图5 某双曲拱坝基础体系三维有限元模型

3.2 有限元应力与等效应力计算结果比较

该坝在基本荷载组合（水库正常蓄水位+相应尾水位+设计正常温降+自重+渗透压力+泥砂压力+波浪压力）下的有限元第一主应力计算结果见图6，图中正值表示拉应力，负值表示压应力。

从图6可以看出：上游坝面最大主拉应力发生在坝踵部位，最大拉应力为3.934 MPa，为应力集中效应所致。由于该坝底拱中心角较小，半中心角仅20°，静水压强又较大，故在坝体下游面靠底部出现1.299 MPa的拉应力。下坝面最大拉应力为1.612 MPa，发生在拱冠梁坝段横缝与基岩交界处，亦为应力集中现象。

有限元等效应力第一主应力计算结果见图7，图中正值表示压应力，负值表示拉应力。

从图7可以看出：上游坝面最大等效应力为1.846 MPa，位于坝踵部位，下游坝面最大等效应力为1.319 MPa。

图6 有限元第一主应力等值线（单位：MPa）

图7 等效应力第一主应力等值线（单位：MPa）

比较图6、7发现，在坝踵、坝趾等应力集中部位有限元等效应力值较有限元应力值显著降低。说明有限元等效应力计算可以有效地对有限元应力计算近基础位置应力过大的结果进行修正。

3.3 有限元应力与等效应力的分布规律差异性分析

上下游坝面有限元应力与等效应力的差值绝对值分布见图8。

由图8可看出，拱坝有限元等效应力与有限元应力主要在坝体建基面向上约3/5倍坝底厚度范围内存在较大的差别，其中在建基面处差值最大，等效应力值较有限元应力值减小53%，沿坝体高度向上，二者差值减小。拱冠梁建基面以上3/10倍坝底厚度层面上，二者差值为1.0 MPa；拱冠梁建基面以上3/5倍坝底厚度层面上，二者差值仅为0.1 MPa。从拱冠梁到左右两岸，差值为0.1 MPa的层高逐渐减小。在复杂形状边界，如溢流堰端部与坝体相交处，等值线密集，差值变化梯度较大。由此可见，在离开坝基面及坝体孔洞边界一定距离的坝体区域内，有限元等效应力与有限元应力十分接近，其绝对差值在0.1 MPa以下；在坝体角缘区等有限元弹性结果失真部位，有限元等效应力较有限元应力显著减小，且不存在应力集中问题。