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

馮龍龍 蘇曉麗 李 星

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

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

1 有限元計算原理

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

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

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

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

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

2 堆石體的本構(gòu)關(guān)系以及在Marc中的實現(xiàn)

2.1 堆石體的本構(gòu)模型

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

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

2.2 在Marc軟件中的實現(xiàn)

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

3 面板堆石壩實例分析

3.1 工程概況

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

3.2 有限元模型及其計算參數(shù)

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

圖1 有限元計算模型

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

表1 面板堆石壩的彈性材料參數(shù)

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

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

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

3.3 計算結(jié)果分析

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

圖2 剖面示意圖

3.3.1 壩體應(yīng)力和位移分析

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

表3 剖面的應(yīng)力極值統(tǒng)計表 （單位：MPa）

圖3 竣工期2號剖面應(yīng)力等值線圖（MPa）

圖4 滿蓄期2號剖面應(yīng)力等值線（MPa）

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

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

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

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

表4 堆石體的位移統(tǒng)計表 （單位：cm）

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

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

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

3.3.2 面板應(yīng)力與變形分析

圖7～8分別是竣工期和滿蓄期面板的應(yīng)力圖.

圖7 竣工期面板應(yīng)力分布（MPa）

圖8 滿蓄期面板應(yīng)力分布（MPa）

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

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

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

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

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

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

4 結(jié) 語

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

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

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

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

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

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

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

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

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

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