基于MSC.Marc軟件的面板堆石壩靜力分析

馮龍龍 蘇曉麗 李 星

（1.河海大學 水文水資源與水利工程科學國家重點實驗室，南京 210098；2.河海大學 水資源高效利用與工程安全國家工程研究中心，南京 210098；3.河海大學 水利水電學院，南京 210098）

混凝土面板堆石壩是由防滲面板，防滲接地結構以及壩體堆石體組成的混合結構，由于混凝土面板堆石壩具有工程量小、安全、經濟、施工方便和適應性好的特點［1］，所以它在實際工程中得到廣泛的應用.隨著技術的發展和設計經驗的積累，面板堆石壩的設計高度已經突破了200m，如水布埡水電站的壩高超過了230m.由于壩高的增加，使得壩體的應力和變形變得比較復雜.因此，在設計混凝土面板堆石壩時，對其進行應力和變形分析是十分重要的.本文計算時以某高面板堆石壩為例，利用Marc軟件的非線性求解功能進行應力和變形分析，給出其相應的變化規律.

1 有限元計算原理

非線性有限元按位移模式求解的平衡方程為

式中，［K（u）］為勁度矩陣，｛u｝為節點位移列陣，｛R｝為節點荷載列陣.

本次計算采用中點增量法［2］求解非線性方程組.壩的填筑順序把自重荷載分為若干級，配合荷載分級進行有限元網格劃分.中點增量法把每級自重荷載稱為荷載增量，記為｛ΔR｝，如果第i級為｛ΔR｝i，先根據上一級荷載計算末的｛σ｝i－1確定彈性常數 Ei－1和νi－1，并組成勁度矩陣［K］i－1，在結構中施加本級荷載增量的一半｛ΔR｝i／2，用式（2）求得位移增量：

再計算應力與應變增量，疊加到上次計算的應變與應力上，由此得到本級計算的中點應變｛ε｝i－1／2，然后可以得到彈性常數Ei－1／2和νi－1／2.再用式（3）求本級全荷載發生的位移增量：

最后計算本級的應變與應力增量，累加到上一級終了的應變與應力上，即為本級的應變與應力，壩體上游面的水壓力按照蓄水順序分級加荷.

2 堆石體的本構關系以及在Marc中的實現

2.1 堆石體的本構模型

堆石體的本構關系表達了堆石體應力應變之間的關系，它無疑對應力應變的計算結果起決定性的作用［3］，因而采用恰當的本構模型對面板堆石壩應力變形分析是至關重要的.堆石料是非線性材料，根據三軸實驗，鄧肯－張（E－B）模型能夠較好地反映土體材料應力應變的非線性，因此選擇鄧肯－張（E－B）模型作為堆石料、墊層和過渡料等的本構模型.該本構模型的基本原理［4］為

式中，Rf為破壞比；Pa為大氣壓力；K、Kb、m、n為試驗常數；S為剪應力水平，反映材料強度的發揮程度，S＝，（σ1－σ3）f為破壞時的偏應力，由摩爾－庫倫破壞準則得：

2.2 在Marc軟件中的實現

Marc軟件作為大型有限元商業計算軟件之一，與其它軟件相比，具有能夠快速求解高度非線性問題的優點.但是，它不包含工程中常用的鄧肯－張（E－B）非線性彈性模型，這嚴重影響了它在壩工領域中的使用.本文利用 Marc軟件提供二次開發功能［5］，把通過實驗得到的參數利用Fortran語言編寫子程序，利用Marc子程序來實現鄧肯－張非線性彈性模型的應用，然后利用生死單元技術進行分級加載，以進行壩體施工和蓄水的模擬.

3 面板堆石壩實例分析

3.1 工程概況

某水庫工程位于河南省境內，它的開發任務主要是以防洪、供水為主，同時考慮發電和改善生態環境，并進一步完善黃河下游調水調沙運行條件.該工程為混凝土面板堆石壩，設計壩高124.0m，壩頂高程288.5m，頂部設防浪墻，其高度為1.2m，壩頂總長為481.0m，寬為10.0m，設計時上游采用1∶1.5的坡率，下游的為1∶1.6.壩體主要由主堆石區、次堆石區、混凝土面板和防滲墻等幾部分組成，主堆石區與次堆石區分界線為壩軸線向下游1∶0.6.根據面板壩規范要求，100m以上的壩面板上游面下部設有上游鋪蓋區及蓋重區，上游鋪蓋的頂高程，取為壩體高度的35%，即頂部高程在200.00m，最大高度為35 m；在周邊縫下設小料區（特殊墊層）；下游坡面為大塊石護坡.

3.2 有限元模型及其計算參數

考慮壩體分區、施工程序及加載過程，并考慮到防滲墻的連接型式，對壩體及壩基進行剖分，建立三維有限元模型如圖1所示，總共剖分9 862個單元，11 489個結點.

圖1 有限元計算模型

混凝土面板、趾板、連接板等在達到破壞強度之前線性關系一般較好，故按線彈性材料處理，參數見表1.

表1 面板堆石壩的彈性材料參數

面板堆石壩中主堆石料、次堆石料、粘土夾層和過渡石料等按照非線性彈性材料進行模擬，計算使用的參數見表2.

荷載加載時，利用Marc軟件提供的生死單元功能實現分級加載來模擬面板堆石壩的澆筑過程以及蓄水期的情況.按照面板堆石壩施工進度和蓄水計劃，先后將荷載分為24級模擬.

表2 面板堆石壩的鄧肯－張模型（E－B）參數

3.3 計算結果分析

分析計算結果時，取3個剖面1、2、3（如圖2所示）進行分析，由于篇幅有限，文中只給出2號剖面的等值線圖，其它剖面的極值列入相應的表中.計算時選取的坐標方向：x向為順河向，y方向為壩高方向，向上為正，z軸為壩軸線方向，從左向右；應力按照土力學的習慣，以壓應力為正，拉應力為負.

圖2 剖面示意圖

3.3.1 壩體應力和位移分析

圖3為竣工期堆石體2號斷面的應力分布，圖4為滿蓄期堆石體2號剖面的應力等值線圖，表3分別列出了3個剖面的應力極值.

表3 剖面的應力極值統計表 （單位：MPa）

圖3 竣工期2號剖面應力等值線圖（MPa）

圖4 滿蓄期2號剖面應力等值線（MPa）

由圖3～4和表3可以得到壩體應力分布規律：竣工期壩體主應力等值線與壩坡基本平行，從壩頂向壩基呈逐漸加大的趨勢，堆石體應力極大值出現在2號剖面，第一主應力為2MPa，第三主應力為0.80 MPa；蓄水后，受水荷載作用，堆石應力極值增大，所處的位置進一步向上游主堆石區靠近，主堆石區應力較次堆石區應力有所增加，第一主應力最大值為2 MPa，第三主應力最大值為0.85MPa，極值仍然出現在2號剖面.

圖5～6分別為2號斷面竣工期和滿蓄期堆石體的位移等值線圖，3個斷面的位移極大值列于表4.

圖5 竣工期2號剖面位移等值線圖（cm）

圖6 滿蓄期2號剖面位移等值線圖（cm）

表4 堆石體的位移統計表 （單位：cm）

由圖5～6以及表4可以看出：由于堆石體的泊松效應，使得橫向剖面上水平位移分布規律基本上是上游堆石區位移指向上游，下游堆石區位移指向下游，這符合竣工期面板堆石壩上下游方向位移分布的一般規律.

1號剖面向上游最大位移為3cm，位于上游側剖面1／3壩高處；向下游最大位移為7cm，位于下游側剖面1／3壩高處；最大豎向位移為40cm，位于剖面的1／2壩高處.2號剖面向上游最大位移為9cm，位于上游側基礎覆蓋層以上壩體1／4壩高位置，向下游最大位移為23cm，位于基礎覆蓋層以上壩體1／3壩高位置；最大豎向位移為96cm，位于剖面的1／3壩高處.3號剖面向上游最大位移為3cm，位于上游側剖面1／3壩高處；向下游最大位移為9cm，位于下游側剖面1／3壩高處；豎直向位移極值為42cm，位于斷面的1／2壩高處.

水庫蓄水后，在水荷載的作用下，壩體3個剖面的順河向位移均有較大變化，但仍然有一定的規律.1號剖面向上游最大位移為1.4cm，位于壩底附近；向下游最大位移為9.3cm，位于下游側剖面1／3壩高處；豎直向位移極值為40cm，位于剖面的1／2壩高處.2號剖面向上游最大位移為2cm，位于上游堆石區壩底部位置，向下游最大位移為27cm，位于下游次堆石區壩底部附近；最大豎向位移為98cm，位于剖面的1／3壩高處.3號剖面向上游最大位移為1.4cm，位于壩底部附近；向下游最大位移為11cm，位于下游側剖面1／3壩高處，豎直向位移為43cm，位于剖面的1／2壩高處.

3.3.2 面板應力與變形分析

圖7～8分別是竣工期和滿蓄期面板的應力圖.

圖7 竣工期面板應力分布（MPa）

圖8 滿蓄期面板應力分布（MPa）

由圖可知：竣工期，面板順坡向主要表現為壓縮變形，最大壓應力為3.8MPa，出現在2號剖面182m高程位置.面板軸向應力與順坡向應力相比數值較小，基本呈兩岸受拉、中間受壓分布，最大壓應力為0.8MPa，出現在2號剖面180m高程位置，最大拉應力值為1MPa，出現在3號剖面190m高程位置.

蓄水后，面板順坡向應力仍然以壓為主，最大壓應力為4.1MPa，出現在2號剖面182m高程位置.但靠近河床兩端出現較大的拉應力區，極值為1.8 MPa，出現在面板的最底部.相對于竣工期，軸向應力有所增大，但總體呈兩岸面板受拉，中間受壓分布，符合軸向應力的一般分布規律：2號剖面高程217m位置軸向壓應力最大值為1.5MPa；相對于竣工時面板的拉應力明顯增大，而且在靠近河床的兩端及兩岸周邊處會出現較大的拉應力，最大值為2.5MPa.

圖9為竣工期面板位移分布等值線圖，圖10為滿蓄期面板位移分布等值線圖.

圖9 竣工期面板位移等值線圖（cm）

圖10 滿蓄期面板位移等值線圖（cm）

由圖9～10可以看出：竣工期，由于已經歷了一期面板擋水，因此面板撓度基本指向壩內，位于河床中央的面板底部撓度較大，撓度最大值為7cm；面板在壩軸向上由河谷向中央擠壓［6］，左岸軸向位移最大值為1.6cm（指向右岸）；軸向位移在右岸的最大值為1.8cm，方向指向河床中央.蓄水后，面板變形分布規律較好，面板撓度指向壩內，面板中間區域數值較大，最大值為24cm；面板軸向有進一步由河谷向中央擠壓的趨勢，軸向位移在左岸最大為3.4cm，右岸為3.7 cm，且分布具有一定的對稱性.

4 結 語

本文通過Marc軟件提供的子程序功能，利用Fortran語言編寫子程序以實現非線性彈性E－B模型在Marc中的應用.根據三維非線性有限元分析可以得出：

1）三維情況下，壩體具有較強的空間效應，河床中央剖面的最大主應力值不在壩基面底部，而位于距離底部有一定距離的位置.

2）由于基礎有較深的覆蓋層，導致堆石體的最大沉降發生的位置以及面板最大撓度出現的位置相對于修建在基巖上的面板壩位置偏低一些，而且該工程由于次堆石下部參數較低、變形較大，故造成最大值位置偏向下游次堆石下部附近.

3）面板變形分布規律較好，面板撓度指向壩內，壩軸向上由河谷向中央擠壓；蓄水后面板應力仍然以壓為主，但是在靠近河床兩端出現較大的拉應力，在實際工程中應予以重視.

［1］ 孔憲京，張 宇，鄒德高.高面板堆石壩面板應力分布特性及其規律［J］.水利學報，2013，44（6）：631－639.

［2］ 田 艷，王瑞駿，呂海東.積石峽面板堆石壩應力變形非線性有限元分析［J］.水資源與水工程學報，2007，18（5）：71－74.

［3］ 高 輝，李守義.面板堆石壩非線性有限元應力變形分析［J］.西北水資源與水工程，2000，11（4）：33－37.

［4］ 付國棟，王瑞駿，趙一新.施工期壩面過水面板堆石壩非線性有限元分析［J］.水資源與水工程學報，2009，20（6）：116－119.

［5］ 陳火紅，尹偉奇，薛小香.MSC.Marc二次開發指南［M］.北京：科學出版社，2004.

［6］ 張 峰.深覆蓋層面板堆石壩應力變形分析［J］.工程地質計算機應用，2012，68（4）：30－35.