鄧肯－張E－B模型的ANSYS二次開發(fā)及應(yīng)用

孫明權(quán)，陳姣姣，劉運(yùn)紅

（華北水利水電學(xué)院，河南鄭州 450045）

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

1 計(jì)算原理

1．1 鄧肯張E－B模型

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

1．1．1 彈性模量的計(jì)算

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

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

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

切線體積模量采用下式計(jì)算

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

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

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

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

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

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

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

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

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

1．2．3 ANSYS 重啟動分析

在壩體施工階段，為了進(jìn)行下一級施工，每完成一次施工求解，必須進(jìn)行重啟動分析，以保證結(jié)果的正確．重啟動分析并不保存生死單元的設(shè)置，因此在重啟動分析時應(yīng)當(dāng)重新設(shè)置單元的生死情況．

1．2．4 初始應(yīng)力狀態(tài)的設(shè)置

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

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

1．2．5 計(jì)算步驟

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

1．2．6 對結(jié)果進(jìn)行后處理

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

2 工程實(shí)例

2．1 計(jì)算模型

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

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

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

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

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

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

3 計(jì)算分析

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

3．1 壩體竣工期

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

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

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

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

3．2 壩體初次蓄水

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

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

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

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

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

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

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

4 結(jié)語

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

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

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

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

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

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

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

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

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

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

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

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

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

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