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

（武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072）

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

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

1 基于ANSY S的有限元等效應(yīng)力計(jì)算

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

1.1 基于APDL編制計(jì)算程序

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

1.2 計(jì)算坐標(biāo)系設(shè)置

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

圖1 內(nèi)力計(jì)算時(shí)的局部坐標(biāo)系Ⅰ

圖2 等效應(yīng)力計(jì)算時(shí)的局部坐標(biāo)系Ⅱ

1.3 沿路徑積分

在1.1步驟②中局部坐標(biāo)系Ⅰ下，積分的路徑需設(shè)定為沿拱圈徑向方向。梁的水平截面在拱中心線上取單位寬度，沿厚度方向?qū)α旱膽?yīng)力及其矩進(jìn)行積分；單位高度拱圈的徑向截面，沿厚度方向?qū)皯?yīng)力及其矩進(jìn)行積分，內(nèi)力積分公式見(jiàn)文獻(xiàn)[7]。

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

圖3 路徑示意

1.4 夾角φ、η求解

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

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

圖4 四面體ABCO的應(yīng)力示意

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

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

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

交線AC的方程為

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

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

及

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

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

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

3 工程實(shí)例

3.1 有限元計(jì)算模型

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

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

3.2 有限元應(yīng)力與等效應(yīng)力計(jì)算結(jié)果比較

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

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

有限元等效應(yīng)力第一主應(yīng)力計(jì)算結(jié)果見(jiàn)圖7，圖中正值表示壓應(yīng)力，負(fù)值表示拉應(yīng)力。

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

圖6 有限元第一主應(yīng)力等值線 （單位：MPa）

圖7 等效應(yīng)力第一主應(yīng)力等值線 （單位：MPa）

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

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

上下游壩面有限元應(yīng)力與等效應(yīng)力的差值絕對(duì)值分布見(jiàn)圖8。

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

圖8 第一主應(yīng)力差值絕對(duì)值分布 （單位：MPa）

4 結(jié)語(yǔ)

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

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

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

（4）本方法可為工程設(shè)計(jì)人員以智能化的手段完成拱壩有限元等效應(yīng)力計(jì)算、分析全過(guò)程提供一條有效途徑，方法簡(jiǎn)單實(shí)用。

［1］ 傅作新，錢(qián)向東.有限單元法在拱壩設(shè)計(jì)中的應(yīng)用［J］.河海大學(xué)學(xué)報(bào)， 1991， 19（2）:8-15.

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

［3］ 李同春，章杭惠.改進(jìn)的拱壩等效應(yīng)力分析方法［J］.河海大學(xué)學(xué)報(bào) （自然科學(xué)版）， 2004， 32（1）:104-107.

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

［5］ 朱伯芳.拱壩的有限元等效應(yīng)力及復(fù)雜應(yīng)力下的強(qiáng)度儲(chǔ)備［J］.水利水電技術(shù)， 2005， 36（1）:43-47.

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

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