基于ANSYS的拱壩等效應力計算及圖形顯示

李業(yè)盛，賴國偉，甘海闊

（武漢大學水資源與水電工程科學國家重點實驗室，湖北 武漢 430072）

拱壩是高次超靜定結(jié)構(gòu)，壩體應力對于邊界條件的影響十分敏感，與純拱法、拱梁分載法等結(jié)構(gòu)力學方法相比，有限元法可以考慮壩體大孔口、不規(guī)則外形 (如溢流堰)、拱壩與基巖的相互作用、復雜基礎等因素的影響，計算精度較高。但在彈性范圍內(nèi)，有限元計算成果在壩踵、壩趾等角緣部位存在應力集中現(xiàn)象，且網(wǎng)格愈密應力集中程度愈高，不利于確定應力控制指標。為此我國一些學者提出了 “有限元等效應力”。傅作新等[1]將有限元法所求得的應力合成為截面內(nèi)力，求出對應的線性化應力作為等效應力，但截面內(nèi)力是擬合出來的，不能滿足內(nèi)力平衡條件[2]。李同春等[3]將壩體分解為拱系和梁系，根據(jù)拱和梁的內(nèi)力平衡條件求解指定截面上的約束內(nèi)力，進而求解相應截面上的內(nèi)力和壩體內(nèi)任一點的等效應力。楊強等[4]將建基面沿高程依次分成若干段曲面，對每段曲面均采用等效矩形進行近似，分別積分求解各等效矩形截面的內(nèi)力。朱伯芳[5]認為在計算有限元等效應力時，宜直接進行數(shù)值積分，而不宜用二次曲線逼近，因拱壩是偏心受壓結(jié)構(gòu)，應力分布是一條有拐點的高次曲線。李守義等[6]經(jīng)過對比分析認為,基于ANSYS的拱壩等效應力計算結(jié)果的應力分布規(guī)律與拱梁分載法計算結(jié)果的分布規(guī)律基本一致，結(jié)果可靠，計算精度和效率較高。

盡管ANSYS后處理功能強大，但并不能像有限元計算結(jié)果一樣將等效應力結(jié)果在后處理模塊中直接輸出為工程上常用的、直觀的等值線、云圖等形式，一般需借助TECPLOT、MATLAB等第三方后處理軟件將計算結(jié)果重新整理并輸出[6]，過程略顯繁瑣。本文在獲得ANSYS拱壩有限元計算結(jié)果數(shù)據(jù)庫的基礎上，擬利用ANSYS參數(shù)化編程語言（ANSYS Parameter Design Language,APDL）直接沿拱壩拱和梁斷面進行數(shù)值積分得到截面內(nèi)力，計算上下游壩面各結(jié)點等效應力，同時采用單元轉(zhuǎn)化及結(jié)點輸出列表修改的方法，在ANSYS后處理模塊中一次性完成拱壩上下游面各結(jié)點等效應力的計算和結(jié)果圖形化顯示。

1 基于ANSY S的有限元等效應力計算

與彈性殼體理論相似，假設三個主要應力分量（σx，σy，σz） 沿壩厚為線性分布[1]， 根據(jù)彈性有限單元法求得的拱壩應力分量，沿梁拱斷面直接進行數(shù)值積分，得到梁與拱的內(nèi)力，即可用材料力學方法計算壩面有限元等效應力，據(jù)此按規(guī)范即可進行拱壩應力安全評價。但基于ANSYS的有限元等效應力計算還有幾個關鍵問題需要解決。

1.1 基于APDL編制計算程序

基于ANSYS有限元計算軟件平臺，利用APDL編制開發(fā)拱壩的有限元等效應力的計算程序。其主要計算步驟為：①根據(jù)拱壩體形尺寸，確定有限元等效應力計算結(jié)點的局部坐標系與梁、拱截面幾何參數(shù)；②由有限元法計算的應力分量，分別沿梁、拱截面積分，計算得到梁、拱截面的內(nèi)力；③用材料力學方法計算上下游壩面各結(jié)點有限元等效應力分量；④最后計算上下游壩面各結(jié)點有限元等效應力主應力。

1.2 計算坐標系設置

有限元計算模型整體坐標系采用原點為拱冠梁中面壩頂上游面點在海平面的投影、+X軸指向右岸、+Y軸指向順河流向、+Z鉛直向上的直角坐標系。計算等效應力過程中需先后建立兩套局部坐標系。在1.1步驟②沿路徑積分有限元應力分量計算截面內(nèi)力時，需在壩面結(jié)點i建立如圖1所示局部坐標系Ⅰ，其中x軸平行于拱中心線的切線方向，y軸平行于半徑方向，指向上游為正，z軸為鉛直方向，向上為正，原點在中心線上。在1.1步驟③計算上下游壩面各結(jié)點等效應力時，為適應工程上拉負壓正的應用習慣，采用拱梁分載法的坐標系統(tǒng)，需將局部坐標系Ⅰ繞Y軸順時針旋轉(zhuǎn)180°，局部坐標系Ⅱ如圖2所示。在ANSYS后處理中可用RSYS命令激活局部坐標系。

圖1 內(nèi)力計算時的局部坐標系Ⅰ

圖2 等效應力計算時的局部坐標系Ⅱ

1.3 沿路徑積分

在1.1步驟②中局部坐標系Ⅰ下，積分的路徑需設定為沿拱圈徑向方向。梁的水平截面在拱中心線上取單位寬度，沿厚度方向?qū)α旱膽捌渚剡M行積分；單位高度拱圈的徑向截面，沿厚度方向?qū)皯捌渚剡M行積分，內(nèi)力積分公式見文獻[7]。

以下壩面為例，當計算下壩面上A點時，不難確定出其沿水平拱圈徑向?qū)纳蠅蚊鍮點，以A、B兩點為控制點形成積分路徑，如圖3所示。在ANSYS中沿路徑積分有多種方法，這里介紹一種簡單的思路，用PATH命令定義A、B兩點確定的積分路徑，通過PDEF命令、PVECT命令分別將有限元計算的應力結(jié)果、積分點的位置在局部坐標系I下映射到積分路徑上，最后用PCALC命令沿路徑積分應力分量及其矩，即得到梁與拱的內(nèi)力。

圖3 路徑示意

1.4 夾角φ、η求解

利用梁與拱的內(nèi)力，即可根據(jù)材料力學的基本公式和平衡條件，計算壩體表面的有限元等效應力。如圖4示，在下游壩面取出一四面微元體ABCO，有6個應力分量作用在相互垂直的3個平面ACO、ABO和BCO上，同時有水壓力p作用在壩面ABC上。

水平面上懸臂梁鉛直正應力σz、切向水平剪應力τzx及徑向鉛直平面上拱的水平正應力σx計算公式見文獻[7]，根據(jù)力的平衡條件，可求出未知的3個應力分量如下

圖4 四面體ABCO的應力示意

式中，φ為徑向鉛直平面內(nèi)壩面與鉛直線的夾角，即AC與鉛直向 （Z軸）的夾角；η為水平面內(nèi)壩面與拱中心線切線的夾角，即BC與切向 （X軸）的夾角；p為壩面水壓力。

夾角φ、η可根據(jù)面ABC的法向向量求得，過程為：在局部坐標系Ⅱ下，設下壩面法向向量為（nx， ny， nz）， 該法向向量可在 ANSYS計算程序中用面的投影方式求解提取。今在此法向向量的基礎上任意作一平面，其平面方程為

式中，d為任意常數(shù)。

交線AC的方程為

由此并結(jié)合四面微元體ABCO上的力平衡推導過程，可得

同理，可得交線BC的方程為

及

由式（6）、 （8）求得的夾角 φ、 η 代入式（1）～（3）中即得 τxy、 τyz、 σy三個應力分量， 與已求得的 σz、τzx及σx共6個應力分量，可進一步求得下游壩面上各結(jié)點的主應力值。上壩面的各點主應力值可類似求得。

2 等效應力結(jié)果圖形顯示

在建立有限元計算模型過程中，以上下游壩面結(jié)點集為基礎創(chuàng)建不參與有限元計算的壩面MESH200三維表面虛單元集。在后處理模塊中完成有限元等效應力計算后，將上下游壩面MESH200表面虛單元集轉(zhuǎn)化為能顯示壩面結(jié)點位移分量 （ux，uy，uz） 的SURF65表面效應單元集，再用DNSOL命令將SURF65表面效應單元集的結(jié)點位移分量修改為計算得到的結(jié)點等效應力主應力值 （s1，s2，s3）， 即可直接輸出類似有限元計算結(jié)果一樣的上下游面的等效應力等值線圖或彩色云圖。

3 工程實例

3.1 有限元計算模型

某雙曲混凝土拱壩，壩頂高程355.30 m，壩頂厚度1.7 m，壩底厚度5.0 m，最大壩高53.3 m，最小中心角40°，壩頂弧長126 m，厚高比0.094，屬于薄拱壩。該工程三維有限元網(wǎng)格劃分見圖5。整體模型模擬了左拱座所在相對較單薄的轉(zhuǎn)向山體和壩下游河谷實際地形，以及地基F3、F4斷層。對壩體，模擬了溢流堰、8條臨時施工橫縫及壩體混凝土、漿砌石不同材料分區(qū)。單元采用三維八結(jié)點等參單元，沿壩厚方向設置五層單元，整體網(wǎng)格共48 480個結(jié)點，49 719個單元，其中壩體結(jié)點11 915個，單元9 455個。

圖5 某雙曲拱壩基礎體系三維有限元模型

3.2 有限元應力與等效應力計算結(jié)果比較

該壩在基本荷載組合 （水庫正常蓄水位+相應尾水位+設計正常溫降+自重+滲透壓力+泥砂壓力+波浪壓力）下的有限元第一主應力計算結(jié)果見圖6，圖中正值表示拉應力，負值表示壓應力。

從圖6可以看出：上游壩面最大主拉應力發(fā)生在壩踵部位，最大拉應力為3.934 MPa，為應力集中效應所致。由于該壩底拱中心角較小，半中心角僅20°，靜水壓強又較大，故在壩體下游面靠底部出現(xiàn)1.299 MPa的拉應力。下壩面最大拉應力為1.612 MPa，發(fā)生在拱冠梁壩段橫縫與基巖交界處，亦為應力集中現(xiàn)象。

有限元等效應力第一主應力計算結(jié)果見圖7，圖中正值表示壓應力，負值表示拉應力。

從圖7可以看出：上游壩面最大等效應力為1.846 MPa，位于壩踵部位，下游壩面最大等效應力為1.319 MPa。

圖6 有限元第一主應力等值線 （單位：MPa）

圖7 等效應力第一主應力等值線 （單位：MPa）

比較圖6、7發(fā)現(xiàn)，在壩踵、壩趾等應力集中部位有限元等效應力值較有限元應力值顯著降低。說明有限元等效應力計算可以有效地對有限元應力計算近基礎位置應力過大的結(jié)果進行修正。

3.3 有限元應力與等效應力的分布規(guī)律差異性分析

上下游壩面有限元應力與等效應力的差值絕對值分布見圖8。

由圖8可看出，拱壩有限元等效應力與有限元應力主要在壩體建基面向上約3/5倍壩底厚度范圍內(nèi)存在較大的差別，其中在建基面處差值最大，等效應力值較有限元應力值減小53%，沿壩體高度向上，二者差值減小。拱冠梁建基面以上3/10倍壩底厚度層面上，二者差值為1.0 MPa；拱冠梁建基面以上3/5倍壩底厚度層面上，二者差值僅為0.1 MPa。從拱冠梁到左右兩岸，差值為0.1 MPa的層高逐漸減小。在復雜形狀邊界，如溢流堰端部與壩體相交處，等值線密集，差值變化梯度較大。由此可見，在離開壩基面及壩體孔洞邊界一定距離的壩體區(qū)域內(nèi)，有限元等效應力與有限元應力十分接近，其絕對差值在0.1 MPa以下；在壩體角緣區(qū)等有限元彈性結(jié)果失真部位，有限元等效應力較有限元應力顯著減小，且不存在應力集中問題。

圖8 第一主應力差值絕對值分布 （單位：MPa）

4 結(jié)語

（1）結(jié)合ANSYS參數(shù)化設計語言 （APDL），闡述了拱壩有限元等效應力計算在ANSYS中實現(xiàn)的幾個關鍵問題及解決方法。

（2）通過采用單元轉(zhuǎn)化及結(jié)點輸出列表修改的方法，無需借助第三方軟件，可以將等效應力計算結(jié)果的圖形化顯示與等效應力計算過程在ANSYS后處理模塊中一體化實現(xiàn)，大大提高了計算分析的效率。

（3）經(jīng)實例計算分析、比較發(fā)現(xiàn)，在拱壩壩體的大部分區(qū)域，有限元等效應力與有限元應力十分接近，其差值絕對值在0.1MPa以下；在壩體角緣區(qū)等應力集中部位，有限元等效應力較有限元應力顯著減小，且不存在應力集中問題。

（4）本方法可為工程設計人員以智能化的手段完成拱壩有限元等效應力計算、分析全過程提供一條有效途徑，方法簡單實用。

［1］ 傅作新，錢向東.有限單元法在拱壩設計中的應用［J］.河海大學學報， 1991， 19（2）:8-15.

［2］ 李同春，溫召旺.拱壩應力分析中的有限元內(nèi)力法［J］.水力發(fā)電學報， 2002， 79（4）:18-24.

［3］ 李同春，章杭惠.改進的拱壩等效應力分析方法［J］.河海大學學報 （自然科學版）， 2004， 32（1）:104-107.

［4］ 楊強，劉福深，周維垣.基于直接內(nèi)力法的拱壩建基面等效應力分析［J］.水力發(fā)電學報， 2006， 25（1）:19-23.

［5］ 朱伯芳.拱壩的有限元等效應力及復雜應力下的強度儲備［J］.水利水電技術， 2005， 36（1）:43-47.

［6］ 李守義，周偉，蘇禮邦，等.基于ANSYS的拱壩等效應力研究［J］.水力發(fā)電學報， 2007， 26（5）:38-41.

［7］ 朱伯芳，高季章，陳祖煜，等.拱壩設計與研究 ［M］.北京：中國水利水電出版社，2002.