王俊奇,陸 峰
統(tǒng)一強度理論模型嵌入ABAQUS軟件及在隧道工程中的應(yīng)用
王俊奇1,陸 峰2
(1.華北電力大學(xué)可再生能源學(xué)院水利水電工程系,北京 102206;2.中國水利水電科學(xué)研究院 科研規(guī)劃處,北京 100044)
利用非線性有限元軟件ABAQUS提供的二次開發(fā)功能,實現(xiàn)了統(tǒng)一強度理論本構(gòu)模型的嵌入,完成了有理論解的承受內(nèi)壓作用的厚壁圓筒三維算例數(shù)值測試,以及采用該模型進行隧道開挖三維數(shù)值分析。結(jié)果表明:在ABAQUS中增加統(tǒng)一強度理論本構(gòu)模型,豐富了材料單元庫,提高了計算精度和效率,而且,通過算例驗證和隧道開挖模擬,說明在巖土工程中,考慮材料的中主應(yīng)力效應(yīng),可以充分利用材料強度,指導(dǎo)工程實踐,節(jié)省造價。
統(tǒng)一強度理論;中主應(yīng)力;隧道;ABAQUS有限元
ABAQUS軟件屬國際上最先進的大型通用非線性有限元分析軟件之一,在幾何、材料和接觸非線性問題方面的分析能力居世界領(lǐng)先水平,其中包括了巖土工程中常用的 Drucker-Prager,Mohr-Coulomb和劍橋等眾多材料本構(gòu)模型[1,2]。在巖土工程實踐中廣泛采用的Mohr-Coulomb理論認為,抗剪強度和中主應(yīng)力無關(guān)。不考慮中主應(yīng)力的Mohr-Coulomb準(zhǔn)則對材料強度的反映與試驗結(jié)果有一定出入[3];考慮中主應(yīng)力效應(yīng),可以提高巖石的強度20%~30%[4]。因此,發(fā)展了統(tǒng)一強度理論,通過取用不同參數(shù),考慮中主應(yīng)力的影響[4-6]。
雖然ABAQUS有許多本構(gòu)模型可供選擇,國內(nèi)也開發(fā)了土工分析中常用的 Duncan-Chang模型[7,8],但 Duncan-Chang模型實際屬于非線性彈性模型,處理相對簡單。為彌補ABAQUS中尚無可以考慮中主應(yīng)力的統(tǒng)一強度理論本構(gòu)模型這一不足,本文基于ABAQUS提供的二次開發(fā)用戶子程序,完成了統(tǒng)一強度理論彈塑性本構(gòu)模型的接口開發(fā)工作,與ABAQUS強大的解決大型復(fù)雜土工問題的非線性數(shù)值分析能力結(jié)合,探討考慮中主應(yīng)力對隧道開挖計算塑性區(qū)范圍的影響。
統(tǒng)一強度理論是西安交通大學(xué)俞茂宏教授在1961年提出的雙剪屈服準(zhǔn)則、1983年的雙剪強度理_論和1987年的雙剪角隅模型基礎(chǔ)上,概括提出的工程強度理論體系,在特定條件下,可以退化為現(xiàn)今各種強度理論。考慮材料的拉壓強度差異效應(yīng),同時考慮中主應(yīng)力的影響[4-6]。
統(tǒng)一強度理論可以表示為主應(yīng)力形式,也可表達為用應(yīng)力不變量和應(yīng)力角I1,J2,θ與巖土工程中常用的2個強度參數(shù)c和φ表示的形式[4-6]。
交接處的角度θb可由F=F'的條件求得
加載可以理解為由小屈服面F到大屈服面Φ。加載函數(shù)和加載面不僅與應(yīng)力狀態(tài){σ}有關(guān),而且還取決于塑性應(yīng)變{εp}和反映加載歷史的強化參數(shù)к。用Φ表示加載函數(shù),則加載條件或加載面(后繼屈服面)為 Φ({σ},{εp},κ)=0。
根據(jù)流動法則推導(dǎo)得出彈塑性增量理論本構(gòu)關(guān)系的表達式[9-11]為
其中
ABAQUS為了方便和鼓勵用戶開發(fā)自己研究或感興趣的本構(gòu)模型,采用FORTRAN語言接口方式,提供了若干用戶子程序(User Subroutines)和在編程時可以調(diào)用的實用工具(Utility Routines)。與開發(fā)用戶材料力學(xué)本構(gòu)模型直接相關(guān)的子程序是UMAT,按照約定,開發(fā)者需利用UMAT子程序定義其單元材料積分點的Jacobian矩陣,即材料本構(gòu)關(guān)系的剛度系數(shù)矩陣
用戶材料子程序(user-defined material mechanical behavior,簡稱UMAT)通過與ABAQUS主求解程序的接口實現(xiàn)與ABAQUS的數(shù)據(jù)交流。在輸入文件中,使用關(guān)鍵字“user material”來定義用戶材料屬性。
根據(jù)ABAQUS的約定,用戶子程序體結(jié)構(gòu)應(yīng)至少包括6部分,分別是:ABAQUS約定的子程序題名說明、ABAQUS定義的參數(shù)說明表、開發(fā)者定義的局部變量說明表、開發(fā)者編寫的程序代碼段和子程序返回與結(jié)束語句等,編程代碼目的是得到式(6)的具體表達。
以下就程序設(shè)計中彈塑性增量關(guān)系的一些問題分4小節(jié)作出說明。
彈塑性狀態(tài)確定的一般步驟如下(僅討論各向同性硬化情況):
(1)根據(jù)當(dāng)前增量步或迭代步的位移增量計算應(yīng)變增量
(2)按彈性應(yīng)力應(yīng)變關(guān)系計算應(yīng)力增量的預(yù)測值及應(yīng)力的預(yù)測值(試探性應(yīng)力):
其中tσ是上一增量步或迭代結(jié)束時的應(yīng)力值。
(3)按單元內(nèi)各個積分點計算本增量步或迭代結(jié)束時的t+Δtσ,t+Δtεp等狀態(tài)量。
① 計算屈服函數(shù)值 F(t+Δt珟σ,t珔εp),然后區(qū)分 3種情況,即
若 F(t+Δt珟σ,t珔εp)≤0,則該積分點為彈性加載,或由塑性按彈性卸載,這時均有
若 F(t+Δt珟σ,tεp)>0,且 F(t珟σ,t珔εp)<0,則該積分點為由彈性進入塑性的過渡情況,應(yīng)由
來計算彈性因子m。計算m是為了確定應(yīng)力到達屈服面的時刻。試探應(yīng)力到達屈服面的比例因子m,可以采用線性插值和屈服函數(shù)Taylor展開修正的方法得到[9,10],也可以采用二分的方法獲得,本文采用前者。
若 F(t+Δt珟σ,t珔εp)>0,且 F(t珟σ,t珔εp)=0,則該積分點為塑性繼續(xù)加載,這時令m=0。
②對于①中前兩種情況,均有對應(yīng)于彈塑性部分的應(yīng)變增量 Δε',即
③ 計算彈塑性部分應(yīng)力增量Δσ',即
一般情況下用數(shù)值積分方法進行此積分,在積分過程中可以同時得到Δ珔εp。
④ 計算本增量步或迭代結(jié)束時刻的t+Δtσ,t+Δtεp,即
對率形式的本構(gòu)方程進行積分的算法稱為應(yīng)力更新算法,在彈性預(yù)測階段,塑性應(yīng)變和內(nèi)變量保持固定,而在塑性修正階段,總體應(yīng)變保持不變。
本文采用的是基于顯式積分的切向預(yù)測徑向返回的子增量法[9-12]。
2.5.1 切向預(yù)測徑向返回
所謂切向預(yù)測就是將歐拉方法用于式(13),得到應(yīng)力增量的預(yù)測值,即
進一步得到應(yīng)力的預(yù)測值
同時還可以得到 Δt珔εp。
因為式(16)所表達的算法是顯式的歐拉方法,其中的 Dep(t+Δtσ,t珔εp)是起點切線剛度,所以 Δ珟σ是在加載曲面的切線方向。同時由于加載曲面是外凸的,因此t+Δt珟σ總是在加載面之外。但是屈服準(zhǔn)則要求應(yīng)力t+Δtσ只能在加載曲面之上或者之內(nèi),所以常再采用徑向返回的方法以求得滿足屈服條件的t+Δtσ。具體做法是令
其中r是比例因子,它由下式
得到。
雖然經(jīng)過校正以后得到的t+Δtσ是位于屈服曲面上的,但因為假設(shè)應(yīng)變增量Δε和等效塑性應(yīng)變t+Δt珔εp均保持不變,所以這樣的彈塑性狀態(tài)并不是完全一致的。這種不一致性隨增量步長的增加而增加。為了減小由于這種不一致引起的誤差,可將上述方法和子增量法相結(jié)合。
本文在徑向返回過程中,做了適當(dāng)調(diào)整,采用文獻[9,10]的映射方法直接求得應(yīng)力增量的修正值。
2.5.2 子增量法[12-14]
子增量法就是將總的彈塑性應(yīng)變增量分成若干個子增量,對于每個子增量,利用上述切線預(yù)測徑向返回的方法確定應(yīng)力增量,將每一個子增量結(jié)束時的彈塑性狀態(tài)作為下一個子增量的初始狀態(tài)。子增量有助于提高應(yīng)力增量的計算精度,從而加快迭代的收斂速度。
本文采用的子增量數(shù)N由等效應(yīng)力確定,即
式中:σ0y為初始屈服時對應(yīng)的等效應(yīng)力;珚σy為按試應(yīng)力得到的等效應(yīng)力;σ'y為t+Δtσ=tσ+mΔ珟σ對應(yīng)的等效應(yīng)力。
隱式后向歐拉積分方法中多采用一致切線剛度矩陣。本文采用向前歐拉顯式積分方法,由前面的公式(4)可見,這里形成的彈塑性剛度矩陣是標(biāo)準(zhǔn)的切線模量矩陣,即非一致切線模量矩陣。
統(tǒng)一強度理論在空間屈服面有交線,在π平面上并不光滑,存在奇異點,需特殊處理,一般有2種處理方式,一是Koiter提出的塑性應(yīng)變增量求和法,二是 Owen等建議的回代求導(dǎo)法[4-6,14,15]。本文采用的方法分別如下:
(1)對于θ=θb點產(chǎn)生的奇異性,用矢量和的辦法處理。
(2)對于θ=0°和θ=60°產(chǎn)生的奇異,用回代求導(dǎo)法。
選取一個成熟的問題作為實例,進行對比和驗證。長厚壁圓筒承受內(nèi)壓p=160 MPa作用,內(nèi)壁半徑a=0.1 m,外壁半徑b=0.2 m,如圖1所示(取四分之一),材料彈性模量 E=2.1×105MPa,泊松比ν=0.3,單軸屈服應(yīng)力 σs=240 MPa。
圖1 有限元計算網(wǎng)格Fig.1 Finite element calculation meshes
根據(jù)塑性力學(xué),若采用Tresca屈服準(zhǔn)則,彈性極限壓力和塑性極限壓力分別為[16]
相應(yīng)于所給材料屈服條件,得pe=90 MPa,ps=166.35 MPa,所受內(nèi)壓介于彈性極限與塑性極限之間,按照理論解,圓筒處于彈塑性狀態(tài)。根據(jù)與理論解的比較,驗證數(shù)值試驗結(jié)果的正確性。
首先比較商業(yè)軟件ABAQUS就所給模型的平面和三維計算結(jié)果,采用ABAQUS中Mohr-Coulomb準(zhǔn)則分別計算。長圓筒作為平面應(yīng)變問題,劃分8結(jié)點等參單元單元12個,結(jié)點51個。三維模型沿厚度方向劃為四層單元,采用20結(jié)點等參塊體單元12個,共有結(jié)點122個。網(wǎng)格輪廓見圖1。
Mohr-Coulomb準(zhǔn)則內(nèi)摩擦角取 0°時退化為Tresca準(zhǔn)則,但 ABAQUS軟件中的 Mohr-Coulomb準(zhǔn)則不允許內(nèi)摩擦角取0°,計算時取一小值0.3°,剪脹角取0°,以模擬Tresca準(zhǔn)則。計算得到二維和三維的結(jié)果完全一樣,沿徑向軸的位移和應(yīng)力見表1。
3.3.1 應(yīng)力和位移
按照統(tǒng)一強度理論,當(dāng)b=0時,退化為Tresca準(zhǔn)則。應(yīng)用自編的UMAT程序,就三維模型計算得到沿徑向軸的位移和應(yīng)力見表1。比較表1中的數(shù)據(jù),二者吻合很好,同時說明所編模型程序的正確。
由圖可見,ABAQUS軟件Mohr-Coulomb準(zhǔn)則和b=0時UMAT程序計算得到圓筒切向應(yīng)力基本一致,徑向位移在圖示的結(jié)點上重合。
b增大時,圓筒周邊切向應(yīng)力最大值向筒內(nèi)壁靠近,徑向應(yīng)力變化不大,位移則可以看出隨b的增大而減小。
3.3.2 塑性區(qū)
圖2 圓筒應(yīng)力及位移(ABAQUS程序M-C準(zhǔn)則和UMAT程序b=0,b=0.366,b=1)Fig.2 Stresses and displacements of the cylinder obtained by the M-C criterion of ABAQUSprogram and b=0,b=0.366 and b=1 in UMAT program respectively
表1 三維模型沿水平徑向軸的位移和應(yīng)力Table 1 Displacements and stresses of 3-D model along horizontal radial direction
圖3 圓筒塑性區(qū)分布(等效塑性應(yīng)變大于0.000 1)Fig.3 Distribution of the cylinder plastic zone with the equivalent plastic strain being more than 0.000 1
3.3.3 結(jié)果分析
將統(tǒng)一強度理論本構(gòu)關(guān)系嵌入有限元程序,通過厚壁圓筒承受內(nèi)壓的算例與有關(guān)文獻二維有限元分析結(jié)果對比,在ABAQUS中二次開發(fā)彈塑性程序獲得成功。
某工程隧洞開挖洞徑D=10.2 m,埋深620 m。擬以TBM法施工為主,鉆爆法施工為輔。所選斷面巖石為II級圍巖,主要物理力學(xué)參數(shù)見表2。
表2 材料參數(shù)Table 2 Material parameters
計算時依照對稱原則,取結(jié)構(gòu)的一半進行分析。有限元網(wǎng)格劃分如圖4。為了簡化分析,自地面向下570 m的巖體部分其自重用施加等值面力來代替,約束條件,底部全約束,垂直邊界法向約束。計算時不考慮分步開挖和支護,一次開挖成毛洞,考察位移和塑性區(qū)。
圖4 有限元計算網(wǎng)格Fig.4 Finite element calculation meshes
4.3.1 位 移
統(tǒng)一強度理論中參數(shù)b可以考慮中主應(yīng)力效應(yīng),當(dāng)b=0時,退化為Mohr-Coulomb模型。本工程通過改變b值,考察中主應(yīng)力效應(yīng)對隧道工程的影響。將b=0,b=0.366和 b=0.8時隧道開挖后的水平和豎直位移場作于圖5和圖6,圖中可見,位移場形狀一致,只是數(shù)值有別。將洞周側(cè)壁水平方向及洞頂和洞底豎直方向的位移作于圖7??梢钥闯?,隨b值的增大,3個方向的洞周位移逐漸減小。
圖5 水平方向位移場Fig.5 Displacement fields in the horizontal direction
圖6 垂直方向位移場Fig.6 Displacement fields in the vertical direction
圖7 洞周位移比較Fig.7 Comparison for the peripheral displacements
4.3.2 塑性區(qū)
塑性區(qū)在一定程度上等同于洞室的破壞區(qū),將b=0,b=0.366和 b=0.8時隧道開挖后的塑性區(qū)作于圖8,當(dāng)b=0時,即Mohr-Coulomb法則計算得到洞周塑性區(qū)出現(xiàn)范圍為三層網(wǎng)格,b=0.366時外層網(wǎng)格減小,b=0.8時,兩層網(wǎng)格,說明隨b值的增大,塑性區(qū)范圍在逐漸縮小。
圖8 洞周塑性區(qū)范圍比較(等效塑性應(yīng)變大于0.0001)Fig.8 Comparison of plastic zones of the cylinder with equivalent plastic strain being more than 0.000 1
4.3.3 結(jié)果分析
在現(xiàn)有文獻中也有采用統(tǒng)一強度理論研究洞室開挖的實例[4],但僅限于二維模型,本工程采用三維分析,獲得了與二維模型相似的結(jié)果,即考慮中主應(yīng)力對洞室開挖破壞和位移是有不同的影響的。由工程算例可見,對應(yīng) b=0,b=0.366,和 b=0.8三種條件,洞周位移隨b增大而減小,塑性區(qū)范圍也逐步減小??梢?,隨反映中間主切應(yīng)力及相應(yīng)面上正應(yīng)力對材料破壞影響程度的系數(shù)b值的增大,可使材料強度在工程應(yīng)用中得以充分發(fā)揮,對于指導(dǎo)工程支護措施和節(jié)約工程造價都有很強的實用上的意義。
研究開發(fā)工作表明:依據(jù)ABAQUS為用戶提供的高起點的二次開發(fā)平臺,在ABAQUS中嵌入統(tǒng)一強度理論本構(gòu)模型是可行的。本文給出的具有理論解的典型算例結(jié)果顯示了理想的計算精度。
同時,通過巖土工程中隧道開挖的工程實例計算表明,改變參數(shù)b實現(xiàn)中主應(yīng)力效應(yīng),與不考慮中主應(yīng)力的Mohr-Coulomb破壞準(zhǔn)則比較,顯示增大b值后,洞周位移及塑性區(qū)都有減小趨勢,對于指導(dǎo)工程支護措施和節(jié)約工程造價有積極意義。
[1] ABAQUS/Standard有限元軟件入門指南[M].莊 茁,譯.北京:清華大學(xué)出版社,1998.(Instruction of ABAQUS/Standard[M].ZHUANG Zhuo,translated.Beijing:Tsinghua University Press,1988.(in Chinese))
[2] Hibbitte,Karlsson,SorensonINC(2002).ABAQUS/StandardUser's Manual[S].
[3] 吳天行.土力學(xué)(第二版)[M].馮國棟,譯 .成都:成都科技大學(xué)出版社,1982.(WU Tian-xing.Soil Mechanics[M].FENGGuo-dong,translated.Chengdu:Chengdu U-niversity of Science and Technology Press,1982.(in Chinese))
[4] 俞茂宏 .雙剪理論及其應(yīng)用[M].北京:科學(xué)出版社,1998.(YU Mao-hong.Twin Shear Theory and Application[M].Beijing:Science Press,1998.(in Chinese))
[5] 俞茂宏 .工程強度理論[M].北京:高等教育出版社,1999.(YU Mao-hong.Engineering Strength Theory[M].Beijing:Advanced Education Press,1999.(in Chinese))
[6] 趙均海.強度理論及其工程應(yīng)用[M].北京:科學(xué)出版社,2003.(ZHAO Jun-hai.Strength Theory and Application[M].Beijing:Advanced Education Press,2003.(in Chinese))
[7] 徐遠杰,王觀琪.在ABAQUS中開發(fā)實現(xiàn) Duncan-Chang本構(gòu)模型[J].巖土力學(xué),2004,25(7):1032-1036.(XU Yuan-jie,WANG Guan-qi.Development and implementation of Duncan-Chang constitutive model in ABAQUS[J],Rock and Soil Mechanics,2004,25(7):1032-1036.(in Chinese))
[8] 張 欣,丁秀麗,李術(shù)才.ABAQUS有限元分析軟件中Duncan-Chang模型的二次開發(fā)[J].長江科學(xué)院院報,2005,22(4):45-47.(ZHANG Xin,DING Xiu-li,LI Shu-cai.Secondary development of Duncan-Chang Model in ABAQUSsoftware[J].Journal of Yangtze River Scientific Research Institute,2005,22(4):45-47.(in Chinese))
[9] 朱伯芳.有限單元法原理與應(yīng)用(第二版)[M].北京:中國水利水電出版社,1998.(Zhu Bo-fang.Finite Element Method Principle and Application[M].Beijing:China Hydraulic and Hydroelectric Press,1998.(in Chinese))
[10]王 仁,熊祝華,黃文彬 .塑性力學(xué)基礎(chǔ)[M].北京:科學(xué)出版社,1982.(WANGRen,XIONGZhu-hua,HUANG Wen-bing.Elementary of Plastic Mechanics[M].Beijing:Science Press,1982.(in Chinese))
[11]蔣友諒.非線性有限元法[M].北京:北京工業(yè)學(xué)院出版社,1998.(JIANG You-liang.Nonlinear finite Element Method[M].Beijing:Beijing Institute of Technology Press,1988.(in Chinese))
[12]王勖成.有限單元法[M].北京:清華大學(xué)出版社,
2003.(WANG Xu-cheng.Finite Element Method[M].Beijing:Tsinghua University Press,2003.(in Chinese))
[13]SMITH I M,GRIFFITHS D V.有限元方法編程,王 崧,周堅鑫,王 來,等譯.北京:電子工業(yè)出版社,2003.(SMITH I M,GRIFFITHS D V.Programming the Finite Element Method[M].Translated by WANG Song,ZHOU Jian-xin,WANG Lai,et al.Bejing:Publishing House of Electronic Industry,2003.(in Chinese))
[14]OWEN D R J,HINTON E.塑性力學(xué)有限元-理論與應(yīng)用[M].曾國平,劉 忠,徐家禮,等譯.北京:兵器工業(yè)出版社,1989.(OWEN D R J,HINTON E.Finite Elements in Plasticity[M].ZENGGuo-ping,LIU Zhong,XU Jia-li.trans.Beijing:Publishing Company of Weapon Industry,1989.(in Chinese))
[15]胡其志,周輝,楊雪強,屈服面角點奇異性的非平滑處理對塑性應(yīng)變的影響[J],巖土工程學(xué)報,2009,31(1):66-71.(HU Qi-zhi,ZHOU Hui,YANG Xue-qiang.Effect on plastic strain through the non-smoothness management of corner singularity,Chinese Geotechnical Journal of Geotechnical Engineering,2009,31(1):66-71.(in Chinese))
[16]王 仁,黃文彬,黃筑平 .塑性力學(xué)引論[M].北京:北京大學(xué)出版社,1992.(WANG Ren,HUANG Wen-bin,HUANG Zhu-ping. Introduction of Plastic Mechanics[M].Beijing:Beijing University Press,1992.(in Chinese) )
Unified Strength Theory Constitutive Model Embedded Software ABAQUS and Its Application in Tunnel Engineering
WANG Jun-qi1,LU Feng2
(1.Department of Hydraulic and Hydropower Engineering,North China Electric Power University,Beijing 102206,China;2.Institute of Water Resources and Hydropower Research,Beijing 100044,China)
The unified strength theory constitutive model can be embedded nonlinear finite element software ABAQUSfor providing the user subroutine second development function.A numerical test using thick wall cylinder as an illustration is completed,which gives the theoretical resolution and three dimension analysis of tunnel excavation.All numerical results indicate that the approach of adding the model to the ABAQUSenriches the material element library and enhances the computation efficiency and precision,meanwhile,by the example verification,the tunnel excavation simulation demonstrates that considering influences of middle main stress in geotechnical engineering,the method can fully utilize the material strength enough,guide the engineering practice and reduce investment outlay.
unified strength theory constitutive model;middle main stress;tunnel;ABAQUSfinite element
TU311.41
A
1001-5485(2010)02-0068-07
2009-04-24
中國水利水電科學(xué)研究院開放研究基金資助項目(IWHRO2009005)
王俊奇(1971-),男,山西朔州人,博士,副教授,從事水工結(jié)構(gòu)、巖土工程研究,(電話)15801452959(電子信箱)jqwangmail@163.com。
(編輯:陳紹選)