土石料雙屈服面彈塑性模型的二次開發(fā)算法與應(yīng)用

岑威鈞, 陳司寧, 鄧同春, 熊 堃

(1.河海大學(xué)水利水電學(xué)院, 江蘇 南京 210098; 2. 水利部巖土力學(xué)與工程重點(diǎn)實(shí)驗(yàn)室, 湖北 武漢 430010; 3. 長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究有限責(zé)任公司, 湖北 武漢 430010)

ADINA、ABAQUS和ANSYS等大型商用有限元軟件前后處理方便、計(jì)算求解能力高、擴(kuò)展開發(fā)性強(qiáng),已在巖土工程數(shù)值分析領(lǐng)域中得到了較為廣泛的應(yīng)用[1-2].對(duì)土石壩結(jié)構(gòu)分析而言,目前商用軟件均能方便地實(shí)現(xiàn)壩基分期開挖、壩體分期填筑、壩料分區(qū)設(shè)置及各種接觸關(guān)系的仿真模擬.然而,這些商用軟件的本構(gòu)庫(kù)中僅有Mohr-Coulomb模型、Drucker-Prager模型和(修正的)CAM黏土模型等少數(shù)巖土材料模型,缺乏我國(guó)土石壩工程界廣泛應(yīng)用的Duncan雙曲線非線性彈性模型、沈珠江雙屈服面彈塑性模型(簡(jiǎn)稱沈珠江模型)等實(shí)用土石料本構(gòu)模型,不適合精細(xì)求解土石壩的應(yīng)力變形問題,特別在高土石壩相關(guān)問題研究中受到了較大的應(yīng)用限制.

所幸這些商用軟件都預(yù)設(shè)了二次開發(fā)平臺(tái),用戶可以根據(jù)實(shí)際需要添加所需的本構(gòu)模型:Duncan非線性彈性模型已被開發(fā)應(yīng)用于ABAQUS、ANSYS及ADINA等大型商用軟件中[3-5];“南水”模型已被添加到ABAQUS的本構(gòu)庫(kù)中[6-7];修正劍橋模型也在ANSYS中進(jìn)行了二次開發(fā)[8].相對(duì)于非線性彈性模型的開發(fā),雙曲面彈塑性模型的二次開發(fā)要復(fù)雜許多,涉及到應(yīng)力的精確更新問題,即本構(gòu)積分算法,這是彈塑性本構(gòu)模型二次開發(fā)的重點(diǎn)和難點(diǎn).Wissman[9]和Sloan[10]等較早地提出了次階應(yīng)力積分算法,將應(yīng)變?cè)隽考?xì)分成很多次階,在每個(gè)次階中使用歐拉法、改進(jìn)歐拉法或龍格-庫(kù)塔法進(jìn)行積分得到應(yīng)力增量,然后進(jìn)行疊加;Borija等[11]提出了隱式回歸算法,通過迭代算法使應(yīng)力狀態(tài)返回到屈服面上.回歸算法又包括常彈性模量回歸和變彈性模量回歸,Potts等[12]對(duì)上述兩種應(yīng)力積分算法做了比較,認(rèn)為兩種算法都能得到準(zhǔn)確的結(jié)果,但在多屈服面本構(gòu)開發(fā)以及彈性非線性很強(qiáng)的本構(gòu)模型上,次階算法更為占優(yōu).

為了充分利用商用有限元軟件的優(yōu)勢(shì),拓展其在巖土工程領(lǐng)域中的應(yīng)用范圍,本文詳細(xì)介紹了沈珠江雙屈服面彈塑性模型[13]的二次開發(fā)流程,對(duì)關(guān)鍵的本構(gòu)積分算法進(jìn)行了改進(jìn)和驗(yàn)證,最后進(jìn)行了工程實(shí)例應(yīng)用.由于ADINA、ABAQUS和ANSYS等軟件均用FORTRAN語(yǔ)言進(jìn)行本構(gòu)模型的二次開發(fā),且這些軟件開發(fā)本構(gòu)模型的思想是相同的,因此文中所述二次開發(fā)流程和關(guān)鍵算法同樣可用于上述商用有限元軟件中開發(fā)雙屈服面或多屈服面彈塑性本構(gòu)模型時(shí)提供參考.

1 沈珠江雙屈服面彈塑性模型

沈珠江雙屈服面彈塑性模型[13-14]的屈服面由橢圓曲線和冪曲線組成,其定義的屈服函數(shù)為

(1)

式中:

p和τ8分別為八面體正應(yīng)力和剪應(yīng)力;

t和s分別為橢圓和冪函數(shù)的屈服面系數(shù),不同的土石料t和s取值不同,對(duì)細(xì)粒土建議t=2,s=3,對(duì)堆石料等粗粒土建議t=s=2.

根據(jù)經(jīng)典彈塑性理論,總應(yīng)變可分解為彈性和塑性兩部分.對(duì)于雙屈服面模型,塑性應(yīng)變?cè)隽坑蓛刹糠纸M成,各自對(duì)應(yīng)一個(gè)屈服面.采用正交流動(dòng)法則,則應(yīng)變?cè)隽靠蓪憺?/p>

(2)

式中:

A1和A2分別為對(duì)應(yīng)于屈服面函數(shù)f1和f2的塑性系數(shù),由常規(guī)三軸試驗(yàn)確定;

σij為應(yīng)力張量分量;

Δεe,ij為彈性應(yīng)變?cè)隽?

Δf1和Δf2分別為f1和f2的增量.

按照正交流動(dòng)法則,沈珠江模型的塑性系數(shù)A1和A2反映了兩個(gè)屈服面各自產(chǎn)生的塑性應(yīng)變大小,即塑性勢(shì)面法向增量,均應(yīng)非負(fù)值.在用常規(guī)三軸試驗(yàn)確定A1和A2時(shí),其應(yīng)力比值始終位于兩屈服面法向量之間,可確保A1和A2的非負(fù)性[13].在土石壩的施工和蓄水過程中,目前的研究大多認(rèn)為壩內(nèi)應(yīng)力更符合等應(yīng)力比路徑.有限元模擬計(jì)算發(fā)現(xiàn)壩體某些部位土體單元的應(yīng)力路徑會(huì)在兩屈服面法向量所圍區(qū)域之外,導(dǎo)致A1和(或)A2出現(xiàn)負(fù)數(shù)的不合理現(xiàn)象.為此,在計(jì)算程序中需人為限制塑性系數(shù)A1和A2的非負(fù)性.

設(shè)屈服函數(shù)f1和f2的歷史最大值分別為f1 max和f2 max,則加載準(zhǔn)則按如下定義:(1) 當(dāng)f1>f1 max且f2>f2 max時(shí),為全加載,此時(shí)A1>0、A2>0;(2)f1≤f1 max且f2≤f2 max時(shí),為全卸載,此時(shí)A1=0、A2=0;(3) 當(dāng)f1≤f1 max、f2>f2 max或f2≤f2 max、f1>f1 max時(shí),為部分加載,相應(yīng)有A1=0、A2>0或A1>0、A2=0.沈珠江模型共有黏聚力c、內(nèi)摩擦角φ、破壞比Rf、彈性模量基數(shù)KE、回彈模量基數(shù)Kur、模量指數(shù)nk、圍壓為一個(gè)大氣壓時(shí)的最大體應(yīng)變cd、最大體應(yīng)變隨圍壓變化的冪次nd和最大體應(yīng)變發(fā)生時(shí)的應(yīng)力差與偏應(yīng)力漸近值的比值Rd等 9個(gè)材料參數(shù).

2 二次開發(fā)子程序關(guān)鍵算法

開發(fā)彈塑性本構(gòu)模型時(shí)的關(guān)鍵問題是選擇合適的應(yīng)力積分算法.本文對(duì)常用的基本和中點(diǎn)剛度應(yīng)力積分算法中應(yīng)力比例因子ri的求解進(jìn)行修正改進(jìn),同時(shí)還采用了Sloan提出的帶誤差控制的修正向后Euler返回算法[10]進(jìn)行比較分析.

(1) 計(jì)算彈性預(yù)測(cè)應(yīng)力.在第n步應(yīng)力σn已知的情況下,按虎克定律彈性預(yù)測(cè)第n+1步的試探應(yīng)力σtrial,n+1.

σtrial,n+1=σn+DeΔεn+1,

(3)

式中:

De為彈性剛度矩陣;

Δεn+1為第n+1步的應(yīng)變?cè)隽?

(2) 判斷屈服條件.若彈性試探應(yīng)力σtrial,n+1均不滿足兩個(gè)屈服條件ftrial,i(σtrial,n+1)>fi,max,說明應(yīng)力處于彈性狀態(tài),無需修正,轉(zhuǎn)到第(6)步直接進(jìn)行下一荷載步計(jì)算;除試探應(yīng)力σtrial,n+1在兩個(gè)屈服面內(nèi)的情況以外,其余3種情況表明試探應(yīng)力至少“穿越”了一個(gè)屈服面,需要進(jìn)行應(yīng)力修正.

(3) 計(jì)算修正應(yīng)力.當(dāng)試探應(yīng)力“穿越”屈服面時(shí),需要進(jìn)行應(yīng)力修正,先按線性假設(shè)計(jì)算各屈服面的彈性應(yīng)力比例因子ri為

(4)

式中:

fi(σn)為時(shí)刻n所對(duì)應(yīng)的屈服函數(shù)值,若σn位于該時(shí)刻的屈服面上則為0,反之為負(fù)數(shù);

ftrial,i為彈性試探應(yīng)力對(duì)應(yīng)的屈服函數(shù)值,i=1,2分別對(duì)應(yīng)雙屈服面模型的體積屈服面和剪切屈服面.

對(duì)于雙屈服面彈塑性模型,按線性假設(shè)求得彈性應(yīng)力比例因子ri可能會(huì)有較大誤差,因此需要按式(5)重新計(jì)算,直至相鄰rk,i和rk+1,i之差小于允許值為止.

rk+1,i=rk,i-

(5)

式中:

Δσ為應(yīng)力增量;

r0,i=0,r1,i=1.

彈性比例因子ri做上述迭代修正后發(fā)現(xiàn)算法更準(zhǔn)確有效.

令α=rk+1,i,則應(yīng)力增量αΔσtrail,n+1為彈性分量,應(yīng)力增量(1-α)Δσtrail,n+1為塑性分量.修正應(yīng)力Δσ0只根據(jù)塑性階段計(jì)算.

Δσ0=(1-α)Δσtrail,n+1-Dep(1-α)Δε=

(1-α)(Δσtrail,n+1-DepΔε).

(6)

由式(3)和式(6)可得

Δσn+1=Δσtrail,n+1-Δσ0=

αΔσtrail,n+1+Dep(1-α)Δε=

(De-(1-α)Dp)Δε,

(7)

式(6)、(7)中:

Δε為應(yīng)變?cè)隽?

α為對(duì)應(yīng)“穿越”屈服面的修正系數(shù);

Dep為σn對(duì)應(yīng)的彈塑性矩陣;

Dp為塑性矩陣.單屈服面穿越的情況下,另一個(gè)屈服面的塑性系數(shù)為0.

由式(7)可知,在求第n+1級(jí)應(yīng)力增量時(shí),應(yīng)力應(yīng)變矩陣相當(dāng)于在彈塑性矩陣Dep計(jì)算中將塑性部分乘以“塑性比例因子”(1-α).當(dāng)試探應(yīng)力同時(shí)“穿越”兩個(gè)屈服面,只需令α=(rk+1,1+rk+1,2)/2即可.

(4) 更新應(yīng)力.若采用基本剛度法時(shí),則更新后的第n+1級(jí)應(yīng)力σn+1為

σn+1=σn+Δσn+1.

(8)

為了增加計(jì)算精度,可采用中點(diǎn)剛度法對(duì)第n+1級(jí)應(yīng)力σn+1進(jìn)行2次計(jì)算,“中點(diǎn)”應(yīng)力可按式(9)計(jì)算,取塑性應(yīng)力增量的中點(diǎn)值為

(9)

上述兩種常規(guī)應(yīng)力積分算法在程序代碼編寫時(shí)不進(jìn)行誤差控制,隨著荷載步的增加其應(yīng)力應(yīng)變會(huì)越發(fā)偏離真值,只有在荷載步足夠小時(shí)滿足計(jì)算精度要求.對(duì)于荷載步較大情況,建議使用帶誤差控制的改進(jìn)Euler積分算法.該法是一種子步應(yīng)力積分算法,其核心思想就是按照應(yīng)力誤差對(duì)第n+1級(jí)應(yīng)變?cè)隽坎介L(zhǎng)進(jìn)行動(dòng)態(tài)調(diào)整,即上述計(jì)算中Δεn+1用βkΔεn+1替代,其中調(diào)整系數(shù)βk滿足0<βk≤1和∑βk=1,k為子步數(shù).按式(10)計(jì)算第n+1級(jí)應(yīng)力計(jì)算時(shí)的迭代誤差.

(10)

式中:

Δσk-1,n+1、Δσk,n+1分別為時(shí)刻n+1第k-1和k子步的應(yīng)力增量.

如果計(jì)算得到的應(yīng)力誤差Err大于控制值Eps(一般可取10-3),則采用Sloan的建議對(duì)βk按式(11)進(jìn)行更新.

(11)

繼續(xù)迭代計(jì)算至滿足誤差結(jié)束為止.最終要求在第n+1級(jí)下滿足應(yīng)力誤差的各次迭代計(jì)算的βk之和為1,這樣完成了該級(jí)的應(yīng)力增量更新計(jì)算,轉(zhuǎn)入到下一步.

(5) 計(jì)算彈塑性矩陣.根據(jù)更新后的應(yīng)力σn+1和式(3)計(jì)算第n+1級(jí)彈塑性矩陣Dep.

(6) 重復(fù)上述過程,直至荷載施加完成.

3 算法驗(yàn)證與應(yīng)用

3.1 算例1——三軸試驗(yàn)

采用雅礱江水電站粗粒土的三軸試驗(yàn)資料[7],進(jìn)行固結(jié)排水剪切試驗(yàn)的數(shù)值模擬,試驗(yàn)采用應(yīng)力控制式,圍壓分別按0.5、1.0、2.0 MPa三級(jí)進(jìn)行.

圖1為試驗(yàn)數(shù)據(jù)及數(shù)值模擬結(jié)果,其中:σ1-σ3為主應(yīng)力差;εa為軸向應(yīng)變;εV為體積應(yīng)變.限于篇幅,文中僅給出了圍壓1 MPa和2 MPa 時(shí)3種不同應(yīng)力更新算法的計(jì)算曲線.由圖1可見,不同應(yīng)力積分算法得到的數(shù)值模擬曲線與原試驗(yàn)曲線均吻合程度有所不同,其中以帶誤差控制的改進(jìn)Euler積分算法最優(yōu).

(a) (σ1-σ3)-εa曲線(b) εV-εa曲線(c) (σ1-σ3)-εa曲線(d) εV-εa曲線圖1 三軸試驗(yàn)數(shù)值模擬結(jié)果Fig.1 Numerical simulation results of triaxial test

3.2 算例2——模型面板壩的應(yīng)力變形

本算例為100 m高的模型面板堆石壩,上下游坡坡比均為 1∶1.4,大壩堆石料計(jì)算參數(shù)見文獻(xiàn)[14].利用本文算法在ADINA中二次開發(fā)程序計(jì)算得到的大壩應(yīng)力變形分布規(guī)律與文獻(xiàn)[14]中相應(yīng)結(jié)果一致,限于篇幅,未給出等值線圖,其中各物理量極值對(duì)比見表1.

除壩體向下游水平位移和面板拉應(yīng)力極值與文獻(xiàn)結(jié)果有一定出入外,壩體沉降、壩體大主應(yīng)力和面板撓度極值兩者均很接近,且各物理量極值所在位置完全一致.由此驗(yàn)證本文所用算法及編程過程是精確可靠的.

表1 蓄水期大壩各物理量極值對(duì)比Tab.1 Extreme-value comparison at water-impounding stage

3.3 工程應(yīng)用

某混凝土面板堆石壩的最大壩高為132.5 m,上游壩坡坡比為 1∶1.4,下游平均壩坡坡比為 1∶1.57,混凝土趾板置于弱風(fēng)化基巖上,基巖采用水泥灌漿固結(jié),基礎(chǔ)防滲采用灌漿帷幕.大壩典型剖面圖見圖2,有限元計(jì)算網(wǎng)格見圖3,其中結(jié)點(diǎn)數(shù)為 8 126個(gè),單元數(shù)為7 199個(gè).壩體及覆蓋層計(jì)算參數(shù)見表2,計(jì)算中趾板及面板采用C25混凝土,彈性模量E=28 GPa, 泊松比μ=0.167.

根據(jù)圖2所示的大壩設(shè)計(jì)斷面,結(jié)合壩體填筑及面板分期澆筑施工進(jìn)程,大壩計(jì)算分20級(jí)加載模擬.第1級(jí)模擬壩基、覆蓋層及趾板澆筑,第2～9級(jí)為1期壩體填筑,第10級(jí)模擬1期混凝土面板澆筑,第11～18級(jí)模擬2期壩體填筑,第19級(jí)為2期面板澆筑,第20級(jí)模擬水庫(kù)蓄水,蓄水高程為136.04 m.

圖2 大壩典型剖面圖Fig.2 Typical dam section

圖3 有限元計(jì)算模型Fig.3 Finite element model

竣工期壩體向上、下游的水平位移極值分別為23.1 cm和19.1 cm,沉降極值為90.9 cm,沉降率為0.686%.圖4為蓄水期壩體及覆蓋層變形等值線圖,大壩變形等值線主要在上游發(fā)生變化,向上游的水平位移極值減小至14.1 cm,向下游增至20.6 cm,大壩沉降極值增至91.8 cm,這些數(shù)值符合采用沈珠江模型得到的計(jì)算結(jié)果的一般規(guī)律.

表2 壩體及覆蓋層沈珠江模型計(jì)算參數(shù)Tab.2 Constitutive material parameters

竣工期壩體大、小主應(yīng)力極值分別為2.10 MPa和1.12 MPa,蓄水期壩體大、小主應(yīng)力極值增至2.21 MPa和 1.17 MPa.圖5為蓄水期混凝土面板變形等值線分布圖,面板在壩體變形和自身重力的影響下由兩岸朝河床中央變形.竣工期面板壩軸向變形主要集中在1期面板上,左右岸變形極值分別為1.46 cm和1.24 cm;沉降和撓度極值位于以1、2期面板交接線的河床中心線附近,分別為7.55 cm 和8.55 cm(凸向上游).蓄水后,順河向位移、沉降和撓度極值分別增大,軸向變形極值向壩頂方向移動(dòng),大致位于兩岸壩坡1、2期面板交接處,左右岸變形極值分別為2.79 cm和2.64 cm;沉降極值增至20.2 cm;在庫(kù)水的作用下,撓度極值向趾板方向移動(dòng),大致位于2期面板中心處,其值為21.0 cm(凹向下游).整個(gè)面板無論是在竣工期還是在蓄水期,其變形情況均較好地反映了混凝土面板堆石壩面板的工作特性[15].

(a) 順河向

(b) 豎向圖4 蓄水期壩體及覆蓋層變形等值線(單位:cm)Fig.4 Deformation isoline of dam body and overburden at water-impounding stage (unit: cm)

(a) 壩軸向

(b) 撓度圖5 蓄水期面板變形等值線(單位:cm)Fig.5 Deformation isoline of concrete slab at water-impounding stage (unit: cm)

圖6為蓄水期混凝土面板應(yīng)力等值線分布圖.蓄水后,面板順坡向應(yīng)力基本以受壓為主,極值為5.38 MPa,右岸2期面板頂部小范圍內(nèi)出現(xiàn)微小拉應(yīng)力,其極值為0.33 MPa.面板壩軸向應(yīng)力呈河床中央部位受壓、兩岸部位受拉的特點(diǎn),中部應(yīng)力極值為3.81 MPa,兩岸局部區(qū)域(大致位于左右岸1、2期面板交界處)出現(xiàn)拉應(yīng)力,極值為 1.79 MPa.無論是順坡向還是壩軸向,其面板拉、壓應(yīng)力都分別小于C25混凝土的抗拉和抗壓強(qiáng)度.由于面板垂直縫采用了薄層單元模擬,面板應(yīng)力與變形等值線的光滑程度有一定影響.

面板周邊縫和垂直縫的三向變形規(guī)律及極值符合面板壩變形基本規(guī)律,限于篇幅不再給出接縫三向變形的分布規(guī)律圖.總體而言,壩體和面板及接縫的受力變形均符合彈塑性壩體材料下的基本特性,各物理量的分布規(guī)律和極值也較合理,說明經(jīng)上述驗(yàn)證的二次開發(fā)代碼用于預(yù)測(cè)實(shí)際土石壩工程的應(yīng)力變形是可行和可靠的.

(a) 順坡向

(b) 壩軸向圖6 蓄水期混凝土面板應(yīng)力等值線分布(單位:MPa)Fig.6 Stress isoline of concrete slab at water-impounding stage (unit: MPa)

4 結(jié)束語(yǔ)

土石料本構(gòu)模型的合理選擇很大程度上決定了巖土工程應(yīng)力變形數(shù)值模擬過程的準(zhǔn)確性.本文采用應(yīng)力積分中的基本增量法、中點(diǎn)增量法以及帶誤差控制的修正向后Euler返回算法,詳細(xì)介紹了將沈珠江雙屈服面模型集成到大型通用有限元軟件的本構(gòu)庫(kù)中的具體開發(fā)流程,拓寬了ADINA、ABAQUS和ANSYS等大型商用有限元軟件的應(yīng)用范圍.文中所述二次開發(fā)算法同樣可供用戶開發(fā)多屈服面彈塑性本構(gòu)模型時(shí)參考,主體內(nèi)容和核心代碼是相似的,僅在屈服區(qū)域判斷上稍有不同.

致謝:長(zhǎng)江科學(xué)院開放研究基金資助項(xiàng)目(CKWV2016376/KY); 水利部土石壩破壞機(jī)理與防控技術(shù)重點(diǎn)實(shí)驗(yàn)室開放研究基金(YK914019).

[1] 熊堃,岑威鈞,胡清義,等. 多途徑綜合開發(fā)商業(yè)軟件精細(xì)求解土石壩結(jié)構(gòu)靜動(dòng)力反應(yīng)[J]. 巖石力學(xué)與工程學(xué)報(bào),2013,32(1): 117-125.

XIONG Kun, CEN Weijun, HU Qingyi, et al. Static and dynamic response analysis of earth rockfill dams with united multipath development of commercial softwares[J]. Chinese Journal of Rock Mechanics and Engineering, 2013, 32(1): 117-125.

[2] LIU Changhong, LIU Xiaoling, CHEN Qiu. Interval analysis of the finite element method for stochastic structures[J]. Journal of Southwest Jiaotong University, 2001, 12(1): 46-48.

[3] 熊玉春,房營(yíng)光. ADINA有限元軟件中材料本構(gòu)的二次開發(fā)[J]. 巖土力學(xué),2008,29(8): 2221-2225.

XIONG Yuchun, FANG Yingguang. Secondary development of material constitutive model in ADINA software[J]. Rock and Soil Mechanics, 2008, 29(8): 2221-2225.

[4] 江守燕,謝慶明,杜成斌. 基于ABAQUS平臺(tái)的鄧肯-張E-B和E-v模型程序開發(fā)[J]. 河海大學(xué)學(xué)報(bào):自然科學(xué)版,2011,39(1): 61-65.

JIANG Shouyan, XIE Qingming, DU Chengbin. Development of program of Duncan-Chang E-B and E-v models based on ABAQUS[J]. Journal of Hohai University: Natural Sciences, 2011, 39(1): 61-65.

[5] 孫明權(quán),陳姣姣,劉運(yùn)紅. 鄧肯-張E-B模型的ANSYS二次開發(fā)及應(yīng)用[J]. 華北水利水電學(xué)院學(xué)報(bào),2013,34(2): 30-34.

SUN Mingquan, CHEN Jiaojiao, LIU Yunhong. The second development of Duncan-Chang E-B model in ANSYS and its application[J]. Journal of North China Institute of Water Conservancy and Hydroelectric Power, 2013, 34(2): 30-34.

[6] 岑威鈞,朱岳明. 基于 ABAQUS 的土石料本構(gòu)模型二次開發(fā)及其應(yīng)用[J]. 水利水電科技進(jìn)展,2005,25(6): 78-81.

CEN Weijun, ZHU Yueming. ABAQUS-based secondary development of constitutive model for earth rockfill materials and its application[J]. Advances in Science and Technology of Water Resources, 2005, 25(6): 78-81.

[7] 司海寶,化西婷. 南水模型在ABAQUS中的實(shí)現(xiàn)及在工程中的應(yīng)用[J]. 南水北調(diào)與水利科技,2010,8(1): 52-55.

SI Haibao, HUA Xiting. Development of NHRI constitutive model in ABAQUS and application in engineering[J]. South-to-North Water Transfers and Water Science & Technology, 2010, 8(1): 52-55.

[8] 關(guān)云飛,高 峰,趙維炳,等. ANSYS軟件中修正劍橋模型的二次開發(fā)[J]. 巖土力學(xué),2010,31(3): 976-980.

GUAN Yunfei, GAO Feng, ZHAO Weibing, et al. Secondary development of modified Cambridge model in ANSYS software[J]. Rock and Soil Mechanics, 2010, 31(3): 976-980.

[9] WISSMAN J W, HAUCK C. Efficient elasto-plastic finite element analysis with higher order stress point algorithms[J]. Computers & Structures, 1983, 17(1): 89-95.

[10] SLOAN S W, ABBO A J. Biot consolidation analysis with automatic time stepping and error control. part Ⅰ: theory[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1999, 23: 467-492.

[11] BORIJA R I, LEE S R. Cam clay plasticity, part Ⅰ: implicit integration of constitutive relations[J]. Computer Methods in Applied Mechanics and Engineering, 1990, 78(1): 49-72.

[12] POTTS D, GANENDRA D. An evaluation of substepping and implicit stress point algorithms[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 119(3): 341-354.

[13] 沈珠江. 土體應(yīng)力應(yīng)變分析的一種新模型[C]∥第五屆全國(guó)土力學(xué)及基礎(chǔ)工程學(xué)術(shù)會(huì)議論文選集.北京:中國(guó)建筑工業(yè)出版社,1990.

[14] 冀春樓. 深厚覆蓋層上高堆石壩靜、動(dòng)力分析方法的研究-冶勒堆石壩靜、動(dòng)工作性態(tài)研究[D]. 南京:河海大學(xué),1995.

[15] CEN Weijun, WEN Langshen, ZHANG Ziqi, et al. Numerical simulation on seismic damage and cracking of concrete slab for high concrete face rockfill dams[J]. Water Science and Engineering, 2016, 9(3): 205-211.