土石料雙屈服面彈塑性模型的二次開發算法與應用

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

(1.河海大學水利水電學院, 江蘇 南京 210098; 2. 水利部巖土力學與工程重點實驗室, 湖北 武漢 430010; 3. 長江勘測規劃設計研究有限責任公司, 湖北 武漢 430010)

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

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

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

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

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

(1)

式中:

p和τ8分別為八面體正應力和剪應力;

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

根據經典彈塑性理論,總應變可分解為彈性和塑性兩部分.對于雙屈服面模型,塑性應變增量由兩部分組成,各自對應一個屈服面.采用正交流動法則,則應變增量可寫為

(2)

式中:

A1和A2分別為對應于屈服面函數f1和f2的塑性系數,由常規三軸試驗確定;

σij為應力張量分量;

Δεe,ij為彈性應變增量;

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

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

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

2 二次開發子程序關鍵算法

開發彈塑性本構模型時的關鍵問題是選擇合適的應力積分算法.本文對常用的基本和中點剛度應力積分算法中應力比例因子ri的求解進行修正改進,同時還采用了Sloan提出的帶誤差控制的修正向后Euler返回算法[10]進行比較分析.

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

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

(3)

式中:

De為彈性剛度矩陣;

Δεn+1為第n+1步的應變增量.

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

(3) 計算修正應力.當試探應力“穿越”屈服面時,需要進行應力修正,先按線性假設計算各屈服面的彈性應力比例因子ri為

(4)

式中:

fi(σn)為時刻n所對應的屈服函數值,若σn位于該時刻的屈服面上則為0,反之為負數;

ftrial,i為彈性試探應力對應的屈服函數值,i=1,2分別對應雙屈服面模型的體積屈服面和剪切屈服面.

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

rk+1,i=rk,i-

(5)

式中:

Δσ為應力增量;

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

彈性比例因子ri做上述迭代修正后發現算法更準確有效.

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

Δσ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)中:

Δε為應變增量;

α為對應“穿越”屈服面的修正系數;

Dep為σn對應的彈塑性矩陣;

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

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

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

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

(8)

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

(9)

上述兩種常規應力積分算法在程序代碼編寫時不進行誤差控制,隨著荷載步的增加其應力應變會越發偏離真值,只有在荷載步足夠小時滿足計算精度要求.對于荷載步較大情況,建議使用帶誤差控制的改進Euler積分算法.該法是一種子步應力積分算法,其核心思想就是按照應力誤差對第n+1級應變增量步長進行動態調整,即上述計算中Δεn+1用βkΔεn+1替代,其中調整系數βk滿足0<βk≤1和∑βk=1,k為子步數.按式(10)計算第n+1級應力計算時的迭代誤差.

(10)

式中:

Δσk-1,n+1、Δσk,n+1分別為時刻n+1第k-1和k子步的應力增量.

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

(11)

繼續迭代計算至滿足誤差結束為止.最終要求在第n+1級下滿足應力誤差的各次迭代計算的βk之和為1,這樣完成了該級的應力增量更新計算,轉入到下一步.

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

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

3 算法驗證與應用

3.1 算例1——三軸試驗

采用雅礱江水電站粗粒土的三軸試驗資料[7],進行固結排水剪切試驗的數值模擬,試驗采用應力控制式,圍壓分別按0.5、1.0、2.0 MPa三級進行.

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

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

3.2 算例2——模型面板壩的應力變形

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

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

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

3.3 工程應用

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

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

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

圖3 有限元計算模型Fig.3 Finite element model

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

表2 壩體及覆蓋層沈珠江模型計算參數Tab.2 Constitutive material parameters

竣工期壩體大、小主應力極值分別為2.10 MPa和1.12 MPa,蓄水期壩體大、小主應力極值增至2.21 MPa和 1.17 MPa.圖5為蓄水期混凝土面板變形等值線分布圖,面板在壩體變形和自身重力的影響下由兩岸朝河床中央變形.竣工期面板壩軸向變形主要集中在1期面板上,左右岸變形極值分別為1.46 cm和1.24 cm;沉降和撓度極值位于以1、2期面板交接線的河床中心線附近,分別為7.55 cm 和8.55 cm(凸向上游).蓄水后,順河向位移、沉降和撓度極值分別增大,軸向變形極值向壩頂方向移動,大致位于兩岸壩坡1、2期面板交接處,左右岸變形極值分別為2.79 cm和2.64 cm;沉降極值增至20.2 cm;在庫水的作用下,撓度極值向趾板方向移動,大致位于2期面板中心處,其值為21.0 cm(凹向下游).整個面板無論是在竣工期還是在蓄水期,其變形情況均較好地反映了混凝土面板堆石壩面板的工作特性[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為蓄水期混凝土面板應力等值線分布圖.蓄水后,面板順坡向應力基本以受壓為主,極值為5.38 MPa,右岸2期面板頂部小范圍內出現微小拉應力,其極值為0.33 MPa.面板壩軸向應力呈河床中央部位受壓、兩岸部位受拉的特點,中部應力極值為3.81 MPa,兩岸局部區域(大致位于左右岸1、2期面板交界處)出現拉應力,極值為 1.79 MPa.無論是順坡向還是壩軸向,其面板拉、壓應力都分別小于C25混凝土的抗拉和抗壓強度.由于面板垂直縫采用了薄層單元模擬,面板應力與變形等值線的光滑程度有一定影響.

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

(a) 順坡向

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

4 結束語

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

致謝:長江科學院開放研究基金資助項目(CKWV2016376/KY); 水利部土石壩破壞機理與防控技術重點實驗室開放研究基金(YK914019).

[1] 熊堃,岑威鈞,胡清義,等. 多途徑綜合開發商業軟件精細求解土石壩結構靜動力反應[J]. 巖石力學與工程學報,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] 熊玉春,房營光. ADINA有限元軟件中材料本構的二次開發[J]. 巖土力學,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平臺的鄧肯-張E-B和E-v模型程序開發[J]. 河海大學學報:自然科學版,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] 孫明權,陳姣姣,劉運紅. 鄧肯-張E-B模型的ANSYS二次開發及應用[J]. 華北水利水電學院學報,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 的土石料本構模型二次開發及其應用[J]. 水利水電科技進展,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中的實現及在工程中的應用[J]. 南水北調與水利科技,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] 關云飛,高 峰,趙維炳,等. ANSYS軟件中修正劍橋模型的二次開發[J]. 巖土力學,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] 沈珠江. 土體應力應變分析的一種新模型[C]∥第五屆全國土力學及基礎工程學術會議論文選集.北京:中國建筑工業出版社,1990.

[14] 冀春樓. 深厚覆蓋層上高堆石壩靜、動力分析方法的研究-冶勒堆石壩靜、動工作性態研究[D]. 南京:河海大學,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.