歐陽君 林 飛

（湖南省水利水電勘測設計研究總院 長沙市 410007）

引 言

ABAQUS軟件具有強大的計算功能和廣泛模擬性能，特別是對強非線性和不連續(xù)非線性問題的求解具有明顯的優(yōu)勢[1]。利用ABAQUS軟件進行土石壩地震反應分析時，要涉及到土體的動力本構(gòu)模型的選擇問題。ABAQUS軟件中材料的本構(gòu)模型都是建立在傳統(tǒng)的塑性力學基礎上的，然而傳統(tǒng)的彈塑性理論不能充分反映土料的動力變形機制[2]，同時，軟件中對材料的阻尼的定義采用Rayleigh法和直接模態(tài)阻尼法，按上述方法確定的材料阻尼隨著材料變形的變大而減小，而實際土體材料的阻尼比隨著其變形的變大而增大，且變化幅度約為0.2～30[3]。

本文針對ABAQUS軟件在土石壩地震動力反應分析中的不足，利用其提供的二次開發(fā)用戶子程序接口，完成了反映土體動力變形特性的Hardin-Drnevich[4]本構(gòu)模型的開發(fā)工作。選擇了新疆某土石壩作為研究對象，用鄧肯-張E-B模型完成壩體在穩(wěn)定滲流期的應力與變形計算[5]，再輸入1985年烏恰7.4地震中喀什水電站6.8級余震的實測地震加速度時程曲線[6]進行計算。計算結(jié)果驗證了本文模型的準確性和實用性，為土石壩動力反應分析材料模型的選取提供有益的參考。

1 Hardin-Drnevich本構(gòu)模型

目前用于土石壩地震反應分析的本構(gòu)模型大致可分為兩大類：一類是等效線性粘彈性模型，另一類是真非線性粘彈塑性模型。真非線性粘彈塑性模型根據(jù)加載歷史、屈服條件、硬化規(guī)律、流動法則等建立，在理論上較為合理。由于真非線性粘彈塑性模型比較復雜，雖然相比等效粘彈性模型有一些優(yōu)點[7]，但因其計算參數(shù)較難確定，方程建立及求解工作量大、所需時間長，受應力路徑的影響較大，且目前缺乏合理的計算模型，故工程上的實際應用較少。等效線性粘彈性模型盡管存在一些缺陷，但概念明確、參數(shù)易確定、應用方便、且計算結(jié)果有一定的精度，補充一些相關(guān)的計算模式，包括殘余變形模式，就能夠全面分析地震反應，而且在參數(shù)的確定和應用方面積累了較豐富的試驗資料和工程經(jīng)驗，故被工程界所接受，實用性強。

等效線性粘彈性模型把土視為粘彈性體，采用等效剪切模量Geq和等效阻尼比λeq這兩個參數(shù)來反映土的動應力動應變關(guān)系的兩個基本特征：非線性與滯后性，并且將剪切模量與阻尼比均表示為動應變幅的函數(shù)，即 Geq=G（γ）和 λeq=λ（γ）。 這種模型的關(guān)鍵是要確定最大動剪切模量Gmax與平均有效主應力σ0′的關(guān)系以及剪切模量和阻尼比與動剪應變幅之間的關(guān)系式，國內(nèi)外許多學者對此進行了研究，其中以Hardin-Drnevich模型[4]最有代表性、應用最廣。

Hardin-Drnevich模型假定應力-應變關(guān)系的骨架曲線為一條雙曲線，當動剪應變γ=∞時，曲線以靜力極限剪應力τmax為漸近線。當γ=0時，曲線的切線斜率為最大剪切模量Gmax，見圖1。將此曲線的縱坐標轉(zhuǎn)變?yōu)?γ/τ，繪制成 γ/τ～γ 的關(guān)系線，見圖 2。 此線為一條直線，在縱軸的截距為 1/Gmax，斜率為 1/τmax，有

等效線性剪切模量為：

式中γ為動剪應變幅；γr為參考剪應變，定義γr=τmax/Gmax；Gmax為最大動剪切模量；τmax為最大剪應變。根據(jù)土石料動力試驗資料，最大剪切模量Gmax可以表示為：

式中 σ0′為平均有效主應力，可由 σ′0=[σ′1+（1+K0）σ′3]/3估算。K和n為實驗參數(shù)。τmax由下式可以求得：

式中 σ′1、σ′3分別是最大和最小有效主應力，c′、φ′為有效內(nèi)聚力和內(nèi)摩擦角。

圖1 Hardin-Drnevich模型的主干線

圖2 Hardin-Drnevich模型中的γ/τ～γ關(guān)系線

圖3為Hardin-Drnevich模型的滯回環(huán)，陰影部分的面積為滯回曲線面積的一半,主干線在原點的切線斜率為Gmax，開始卸載時曲線的切線坡度亦為Gmax。假定滯回環(huán)的面積AL為三角形面積Aabc的百分數(shù)，有：

能量損耗系數(shù)的定義為 η=△w/（2πw），△w 為應變一周內(nèi)所損耗的能量，即滯回環(huán)的面積，w為應力應變一周內(nèi)物體內(nèi)部積累的最大彈性變形能，即三角形adc的面積。根據(jù)結(jié)構(gòu)動力學知識η=2λ，加上圖示的幾何關(guān)系和上述假定，得到等效阻尼比為：

當γ接近無窮大時，Geq=0，定義λmax=2K1/π為最大阻尼比，同時根據(jù)式（2）和式（6）可以得出等效阻尼比對動剪應變幅的關(guān)系式：

圖3 Hardin-Drnevich滯回環(huán)模型

由于各種土的性質(zhì)差異很大，試驗表明各類土的應力應變關(guān)系并不完全符合雙曲線的形狀，故需要加以適當?shù)男拚訦ardin[8]（1972年）建議將雙曲線模型修正為如下形式：

其中

式中γh為剪應變，a、b為試驗參數(shù)，與土的振動次數(shù)、振動頻率、土單元的平均有效主應力有關(guān)。

2 Hardin-Drnevich模型的程序?qū)崿F(xiàn)

在進行土石壩的動力反應時，土料的應力-應變關(guān)系通常用八面體上的應力-應變關(guān)系表示[9]，設八面體上剪應力為τ，八面體上的剪應變?yōu)棣?，則剪切變形模量為：

八面體上的剪應變?yōu)椋?/p>

式中J′為應變偏量的不變量，用下式計算：

式中：εx，εy，εz，γxy，γyz和 γzx分別為各單元的正應變和剪應變。

根據(jù)彈性力學原理，應力應變關(guān)系為：

其中：

式中δij稱為克羅納克爾，ν為土的泊松比。

將式（13）寫成增量形式為：

采用隱式積分求解非線性動力問題時，雅克比（Jacobi）矩陣的定義將直接影響到積分格式的收斂速度，上述模型的雅克比矩陣在二維條件下定義如下：

在ABAQUS軟件中的隱式積分模塊Standard和顯式積分模塊Explicit中[10]，分別采用外掛子程序的辦法實現(xiàn)Hardin-Drnevich模型的二次開發(fā)。

3 算例分析

3.1 計算模型

以新疆某位于基巖上的均質(zhì)壩為例進行數(shù)值計算。壩高44.0 m，壩頂寬8.0 m，上游壩坡坡比為1∶2.5，下游壩坡坡比為1:2.75。為獲得動力計算所需的初始應力場，靜力計算采用Duncan-chang（鄧肯-張）E-B模型，模型參數(shù)如表1所示。本文模型進行動力時計算參數(shù)見表2?；鶐r面運動輸入1985年烏恰7.4地震中喀什水電站6.8級強余震的實測地震加速度時程曲線，歷時17.80 s，步長0.01 s。地震動水平向加速度幅值為2.43 m/s2，豎向加速度幅值為2.07 m/s2。水平和豎向地震波時程曲線見圖4。

表1 靜力模型計算參數(shù)

表2 動力模型計算參數(shù)

圖4 喀什地震加速度-時程曲線

3.2 計算結(jié)果與分析

壩體絕對加速度響應是大壩對地震動力作用最直接的反應，為了驗證本文子程序的實用性和準確性，選取壩頂節(jié)點，分析其峰值加速度和放大倍數(shù)，并繪制其絕對加速度響應時程時程曲線見圖5。

圖5 壩頂節(jié)點絕對加速度響應時程曲線

由圖4和圖5可見壩頂節(jié)點水平向和豎向加速度幅值與地震加速度幅值出現(xiàn)的時間相比有一些滯后。計算所得地震加速度反應峰值和壩頂加速度放大倍數(shù)列于表3，由此可見，考慮豎向和水平向地震加速度時的水平向和豎向放大倍數(shù)分別為3.11和1.38，這和現(xiàn)行《水工建筑抗震設計規(guī)范》（SL 203-97）中的5.1.3條所推測設定值相吻合：設計烈度7度時壩頂加速度放大倍數(shù)為3.11與規(guī)范所規(guī)定值3.0相近。以上結(jié)論表明本文基于ABAQUS軟件的二次開發(fā)平臺，利用其子程序UMAT開發(fā)的Hardin-Drnevich動力本構(gòu)模型在土石壩地震反應分析中的研究是可行的。

表3 地震加速度反應峰值及放大倍數(shù)

4 結(jié) 語

基于ABAQUS有限元軟件的二次開發(fā)平臺，編制了Hardin-Drnevich模型的外接子程序，對該子程序的實用性和準確性進行了試驗和理論驗證，說明編制的子程序能較好的模擬土體的動力本構(gòu)特征；選取新疆某土石壩作為研究對象，基于ABAQUS軟件強大的非線性有限元分析環(huán)境，對其進行地震反應分析，給出了壩頂節(jié)點加速度時程曲線；計算結(jié)果表明：基于ABAQUS軟件二次開發(fā)的Hardin-Drnevich動力本構(gòu)模型是可行的，對土石壩地震反應的計算結(jié)果也符合相應規(guī)律；同時，利用ABAQUS軟件強大的非線性求解平臺，進行土石壩的動力反應分析，并利用其用戶子程序UMAT，修改不同類型的Hardin-Drnevich本構(gòu)模型也易于實現(xiàn)。

[1]徐遠杰，王光琪.在ABAQUS中開發(fā)實現(xiàn)Duncan-Chang本構(gòu)模型[J].巖土力學，2004，25（7）：1032-1036.

[2]鄭穎人，沈珠江，龔曉南.巖土塑性力學原理[M].北京:中國建筑工業(yè)出版社，2002.

[3]陳國興，謝君斐，張克緒.土的動模量和阻尼比的經(jīng)驗估計[J].地震工程與工程振動，1995，15（1）：73-84.

[4]顧淦臣.土石壩地震工程 [M].南京:河海大學出版社，1989.

[5]歐陽君,徐千軍,嚴新軍,等.基于ABAQUS軟件的土石壩應力應變分析[J].水資源與水工程學報，2009，20（6）：112-115.

[6]新疆維吾爾地震局.1985年新疆烏恰7.4級地震強余震觀測報告[M].北京:地震出版社，1987.

[7]趙劍明,汪聞韶,常亞屏,陳寧.高面板壩三維真非線性地震反應分析方法及模型試驗驗證[J].水利學報，2003，（9）:12-18.

[8]Hardin B O,Drnevich V P.Shear modulus and damping in soils design equations and curves [J].Journal of Soil Mechanics and Foundation,ASCE,1972,98（7）:603-642.

[9]朱伯芳.有限元法原理與應用[M].北京:中國水利水電出版社，1998.

[10]莊茁,張帆,岑松，等.ABAQUS非線性有限元分析與實例[M].北京:科學出版社，2005.