孫明權，陳姣姣，劉運紅

（華北水利水電學院，河南鄭州 450045）

土石壩歷史悠久，具有就地取材、施工簡單、可以適用多種復雜地質條件等優(yōu)勢，已成為壩工界的主導壩型．隨著壩體高度的增加，除了滲流和穩(wěn)定計算外，其應力和變形分析已經(jīng)成為大型土石壩設計中不可缺少的一部分．有限元法是進行應力和變形分析目前最準確、最有效的方法之一．ANSYS［1］是大型通用的有限元計算軟件，已經(jīng)成為土木建筑行業(yè)分析軟件的主流，但在ANSYS中卻沒有適合于土石壩材料的本構模型．鄧肯－張E－B模型是一種土體的本構模型［2］，被廣泛應用于巖土工程中［3］．筆者通過ANSYS提供的APDL語言二次開發(fā)平臺，實現(xiàn)了在ANSYS中對鄧肯－張E－B模型的模擬．

1 計算原理

1．1 鄧肯張E－B模型

鄧肯（Duncan）和張（Chang）提出的鄧肯－張E－B模型是一種非線性彈性模型，非線性彈性模型是根據(jù)廣義胡克定律建立剛度矩陣D，但考慮到非線性，包含在矩陣D中的彈性常數(shù)E，ν不再視為常量，而是看作隨應力狀態(tài)而改變的變量．

1．1．1 彈性模量的計算

依據(jù)鄧肯－張模型，應當首先判斷單元處于卸荷還是加荷狀態(tài)：當（σ1－σ3）＜（σ1－σ3）0且S＜S0時，單元處于卸荷狀態(tài)，用彈性模量Eur表示;反之，單位處于加荷狀態(tài)，彈性模量用Et表示．（σ1－σ3）0為歷史上曾經(jīng)達到的最大偏應力，S0為歷史上曾經(jīng)達到的最大應力水平．

對于加荷狀態(tài)，切線彈性模量Et為

式中：c為材料凝聚力;φ為材料內摩擦角;pa為標準大氣壓;σ1，σ3分別為單元的大主應力、小主應力;Rf為破壞比;K為彈性模量系數(shù);n為彈性模量指數(shù);Kur為卸載和再加載時的彈性模量系數(shù);nur為卸載和再加載時的彈性模量指數(shù)．1．1．2 體積模量的計算

切線體積模量采用下式計算

式中：Kb為體積模量系數(shù);m為體積模量指數(shù)．

引入切線體積模量Bt后，相當于假定土的泊松比為

只要確定c，φ，Rf，K，n，Kur，nur，Kb，m這幾個參數(shù)，就可以確定鄧肯－張E－B模型，而這些參數(shù)均可由常規(guī)三軸試驗獲得．

1．2 ANSYS的二次開發(fā)原理

1．2．1 APDL語言編寫鄧肯－張E－B模型的宏命令

APDL［4－5］即 ANSYS 參數(shù)化設計語言，用其編寫的鄧肯－張E－B模型的宏命令如下：

1．2．2 ANSYS 模擬土石壩施工分層加載［6－7］

ANSYS模擬土石壩的施工填筑過程，主要利用ANSYS中的“激活”和“殺死”單元．按照填筑順序，首先只激活地基的單元，代表只有地基，荷載是地基的自重;然后在此基礎上激活第一層結構的單元，代表此時施工進行到第一層，荷載是該層結構的自重;進而按照施工順序繼續(xù)激活各層的結構，直至填筑到壩頂．

有限元法利用逐級加荷增量法進行土體的非線性計算．由于計算過程涉及到非線性計算和生死單元，因此應當應用自適應下降關閉的完全牛頓－拉普森選項，即每進行一次平衡迭代，修改剛度矩陣一次．

1．2．3 ANSYS 重啟動分析

在壩體施工階段，為了進行下一級施工，每完成一次施工求解，必須進行重啟動分析，以保證結果的正確．重啟動分析并不保存生死單元的設置，因此在重啟動分析時應當重新設置單元的生死情況．

1．2．4 初始應力狀態(tài)的設置

根據(jù)鄧肯－張E－B模型的公式可以看出，要計算單元的Et和Bt，必須知道單元的σ1和σ3，但每級新填筑層各單元的初始應力狀態(tài)是{σ}=0．如果以此代入鄧肯張公式，則Et=0，無法進行計算．通常采用下面方法確定新填筑層的初始應力狀態(tài)［8－9］：

式中：γ為填土的重度;h為單元形心在土層表面以下的深度;K0為土的靜止側壓力系數(shù);φ為此種材料的內摩擦角．

1．2．5 計算步驟

首先建立土石壩三維模型;然后通過控制單元生死來模擬土石壩的分層施工，在每一層施工完成后通過編制的宏命令來提取各個活單元的最大、最小主應力，執(zhí)行宏修改每個單元的彈性常數(shù);再把當前填筑高度所計算的結果作為下次繼續(xù)計算的初始條件進行重啟動計算．從而可動態(tài)模擬土石壩的施工過程，并在每一層填筑的過程中動態(tài)修改土石壩的彈性常數(shù)，進而可實現(xiàn)鄧肯－張本構模型在土石壩中的應用．

1．2．6 對結果進行后處理

后處理程序的基本設計思想［10］是：將沒有填筑部分的壩體位移歸零．實現(xiàn)方法如下：假設壩體共分n層進行填筑，進行了n步計算，則每一步都包含所有壩體單元的計算結果．對于第一層，第n步的計算結果中的第一層單元的結果即為真實結果;第二層單元的結果需要用第n步的計算結果減去第1步的計算結果得到;同理，第i步的計算結果應當用第n步的計算結果減去第（i－1）步的計算結果得到．假如要求蓄水后正常運行期位移，則可以利用前述步驟求得的結果，加上水壓力單獨作用下壩體產(chǎn)生的位移．

2 工程實例

2．1 計算模型

安寧水電站攔河壩為瀝青混凝土心墻堆石壩，壩頂高程為2 135．00 m，最大壩高62．00 m，壩頂寬度10．00 m，壩頂軸線長度336．00 m．瀝青混凝土心墻厚度為0．80 m，心墻基礎位于混凝土墊座上，墊座寬3．00 m，高3．00 m．壩體河床覆蓋層最大深度約為92．00 m．壩基防滲形式采用混凝土防滲墻（壩基覆蓋層）和帷幕灌漿（壩基基巖），混凝土防滲墻最大深度為87．00 m，厚度為1．00 m．攔河壩筑壩材料分區(qū)從上游到下游分為上游干砌石護坡（厚度為1．00 m）、上游墊層（厚度為 0．50 m）、上游堆石區(qū)、上游2．00 m厚過渡層、瀝青混凝土心墻（厚度為0．80 m）、下游 3．00 m 厚過渡層、下游堆石區(qū)、下游墊層（厚度為0．50 m）、下游干砌石護坡（厚度為1．00 m）．在下游壩基設置厚為1．00 m的過渡料和厚為1．00 m水平反濾層．水庫設計正常蓄水位高程為2 130．00 m．

計算選取具有代表性的土石壩斷面建立二維有限元模型．共剖分2 671個節(jié)點，2 617個單元．分8個荷載步，地基為第1荷載步，壩體填筑分為6個荷載步，蓄水到正常蓄水位為1個荷載步．邊界條件：底部取豎直約束，上、下游取水平約束．計算取定的坐標系：x為順河向，指向下游;y為豎直向，方向向上．

2．2 計算工況及計算參數(shù)

計算工況：竣工期（壩體自重）;蓄水期（壩體自重、浮托力、水壓力），蓄水期上游正常蓄水位57．00 m．計算參數(shù)見表1和表2．

表1 覆蓋層及壩體土石料計算參數(shù)

表2 混凝土及基巖參數(shù)

3 計算分析

竣工期及蓄水期對應壩體及覆蓋層的大小主應力、位移如圖1至圖8所示．圖形只截取中間較重要的一段，自上游壩角向上游延伸132．42 m，自下游壩角向下游延伸138．03 m．圖中應力的符號規(guī)定為：壓應力為正，拉應力為負，單位為MPa;位移的符號規(guī)定為：豎向位移以向下為正，水平位移以向下游方向為正，單位都為cm．

3．1 壩體竣工期

竣工期壩體最大沉降為50．70 cm，占最大壩高

圖1 竣工期壩體及覆蓋層豎直位移（單位：cm）

竣工期大主應力如圖3所示．可以看出在防滲墻以及防滲帷幕與基巖接觸部位附近，有應力集中現(xiàn)象．壩體大主應力的最大值為1．48 MPa，發(fā)生在壩底靠近心墻的位置．竣工期小主應力如圖4所示．

圖3 竣工期壩體及覆蓋層大主應力（單位：MPa）

3．2 壩體初次蓄水

初次蓄水后壩體的最大沉降為56．60 cm，與竣工時相比沉降了5．90 cm，發(fā)生的位置仍然是在上、下游壩體中部與覆蓋層鄰近的位置處，如圖5所示．

竣工期大主應力如圖7所示．可以看出在防滲墻以及防滲帷幕與基巖接觸部位附近，有應力集中現(xiàn)象．壩體大主應力的最大值為1．58 MPa，依然發(fā)生在壩底靠近心墻位置處．竣工期小主應力如圖8所示．可以看出蓄水后，由于水荷載的作用，上游壩殼內的小主應力急劇減少，心墻及下游壩殼內的小主應力有所增大，不再是對稱分布．壩體小主應力的最大值為0．47 MPa，發(fā)生在下游堆石體下方靠近覆（62．00 m）的0．82%，發(fā)生在上、下游壩體中部與覆蓋層鄰近的位置處，如圖1所示．向上游水平位移最大值2．66 cm，占最大壩高的0．04%，發(fā)生在上游覆蓋層內;向下游水平位移最大值7．93 cm，占最大壩高的0．13%，發(fā)生在下游覆蓋層內，如圖2所示．可以看出小主應力關于心墻和防滲墻對稱分布，壩體小主應力的最大值為0．24 MPa，發(fā)生在圍堰下方接近覆蓋層位置處．與竣工期相比，蓄水后壩體及覆蓋層的水平位移變化較為明顯，覆蓋層的向下游水平位移達到了19．87 cm，占最大壩高0．32%，如圖6所示．蓋層處．

圖2 竣工期壩體及覆蓋層水平位移（單位：cm）

圖4 竣工期壩體及覆蓋層小主應力（單位：MPa）

圖6 初次蓄水壩體及覆蓋層水平位移（單位：cm）

圖7 初次蓄水壩體及覆蓋層大主應力圖（單位：MPa）

圖8 初次蓄水壩體及覆蓋層小主應力圖（單位：MPa）

4 結語

1）在ANSYS軟件中，開發(fā)鄧肯－張E－B模型是完全可行的，ANSYS提供的二次開發(fā)功能能夠滿足實際的需要．

2）利用ANSYS二次開發(fā)功能進行的土石壩計算程序理論充足，計算過程由ANSYS自生的非線性求解器實現(xiàn)，保證了計算結果的準確性．

3）通過對安寧水電站瀝青混凝土心墻堆石壩的計算，結果表明，此程序能充分反映土石壩材料的特性和施工過程．

4）通過后處理之后，最終輸出的等值線能夠反映土石壩的豎直沉降和水平位移，符合一般規(guī)律，表明程序是可行的．

［1］張勝民．基于有限元軟件ANSYS7．0的結構分析［M］．北京：清華大學出版社，2003．

［2］錢家歡，殷宗澤．土工原理與計算［M］．北京：中國水利水電出版社，2000．

［3］戴躍華，薛繼樂．ANSYS在土石壩有限元計算中的應用［J］．水利與建筑工程學報，2007，5（4）：74 －77．

［4］龔曙光，謝桂蘭．ANSYS操作命令與參數(shù)化編程［M］．北京：機械工業(yè)出版社，2004．

［5］博弈創(chuàng)作室．APDL參數(shù)化有限元分析技術及其應用實例［M］．北京：中國水利水電出版社，2004．

［6］張愛軍，謝定義．復合地基三維數(shù)值分析［M］．北京：科學出版社，2004．

［7］朱伯芳．有限單元法原理與應用［M］．2版．北京：中國水利水電出版社，1998．

［8］陳慧遠．土石壩有限元分析［M］．南京：河海大學出版社，1988．

［9］范泳賢，劉芳．鄧肯－張E－B模型的ANSYS二次開發(fā)及其應用［OL］．［2009－02－12］．中國科技論文在線，http：//www．paper．edu．cn/releasepaper/content/200902 －568．

［10］彭國倫．FORTRAN 95程序設計［M］．北京：中國電力出版社，2002．