基于統(tǒng)一強(qiáng)度模型的圓形硐室應(yīng)力分析

基于統(tǒng)一強(qiáng)度模型的圓形硐室應(yīng)力分析

1基于abaqus的數(shù)值分析中國力學(xué)專家俞茂宏從滑動位移橫截面的力學(xué)模型出發(fā)，考慮到作用在兩個剪切位移體上的所有應(yīng)力分量對材料的屈服或破壞的不同影響，提出了一套可適用于各種巖石和材料的統(tǒng)一強(qiáng)度理論和統(tǒng)一形式的數(shù)學(xué)方案。Mohr-Coulomb強(qiáng)度理論和雙剪強(qiáng)度理論均為其特例,并且還包含了可以比D-P準(zhǔn)則更合理的新的計算準(zhǔn)則,以及可以描述非凸極限面試驗結(jié)果的新的非凸強(qiáng)度理論。同時統(tǒng)一強(qiáng)度理論可以很好地考慮中間主應(yīng)力對巖體強(qiáng)度的影響,并能與巖土材料的真三軸試驗結(jié)果相吻合,從而在巖土工程領(lǐng)域得到廣泛的應(yīng)用。通用有限元軟件ABAQUS可以分析各種固體力學(xué)、結(jié)構(gòu)力學(xué),特別是能夠處理材料非線性、幾何非線性、接觸非線性以及這三種非線性各種組合的高度非線性問題。ABAQUS本身帶有較多材料非線性彈塑性模型,包括D-P模型、Mohr-Coulomb模型、劍橋模型等,但巖土工程領(lǐng)域中采用的統(tǒng)一彈塑性本構(gòu)模型并沒有在ABAQUS中得以實(shí)現(xiàn),這使得ABAQUS在巖土工程數(shù)值分析中難以充分發(fā)揮,加入了統(tǒng)一彈塑性本構(gòu)模型,則可以充分發(fā)揮ABAQUS軟件的計算能力。為了彌補(bǔ)這一不足,本文利用用戶材料子程序(UMAT)接口,通過Fortran或者VC編程來開發(fā)了統(tǒng)一彈塑性本構(gòu)模型,對單軸壓縮試驗的數(shù)值模擬及圓形硐室問題進(jìn)行了彈塑性分析,并與ABAQUS軟件自帶模型及解析解進(jìn)行比較,以驗證編制的接口程序的正確性。2彈塑性剛度矩陣及點(diǎn)將統(tǒng)一強(qiáng)度理論表示成應(yīng)力不變量的形式為:式中:c0為內(nèi)黏聚力;?為內(nèi)摩擦角;b為中間主應(yīng)力系數(shù)。如圖1所示,統(tǒng)一強(qiáng)度理論包含了多個強(qiáng)度理論:當(dāng)0≤b≤1時為統(tǒng)一強(qiáng)度理論(外凸理論);當(dāng)b<0或b>1時為雙剪非凸強(qiáng)度理論(非凸理論);當(dāng)b=0時統(tǒng)一強(qiáng)度理論退化為Mohr-Coulomb強(qiáng)度理論。本文研究僅考慮0≤b≤1的外凸理論。由屈服函數(shù)F(F′)和塑性勢函數(shù)Q(Q′)決定的彈塑性剛度矩陣可以寫成:式中:A與H′為硬化函數(shù);為等效塑性應(yīng)變。F(F′)是統(tǒng)一強(qiáng)度理論的表達(dá)式。如果采用Q=F(Q′=F′),由式(5)可以得到相關(guān)聯(lián)流動的統(tǒng)一強(qiáng)度彈塑性剛度矩陣。文獻(xiàn)給出了相關(guān)聯(lián)流動的統(tǒng)一彈塑性流動矢量以及本構(gòu)模型中奇異點(diǎn)的數(shù)學(xué)處理。定義統(tǒng)一強(qiáng)度理論的流動矢量為同理,定義屈服面F′流動矢量為如圖1所示,統(tǒng)一彈塑性本構(gòu)模型的屈服面應(yīng)具有單一的硬化參數(shù)和流動方向,但在A、B、C點(diǎn)存在奇異性,定義A、B、C點(diǎn)為奇異點(diǎn)。奇異點(diǎn)是指流動矢量{a}在該點(diǎn)不能唯一確定的點(diǎn)。這里定義兩類奇異點(diǎn),針對不同的奇異點(diǎn),采用不同的處理方法。(1)對于在θ=θb,如圖1中B點(diǎn)產(chǎn)生的奇異性,采用矢量平均的辦法,即(2)當(dāng)b=1時,在點(diǎn)θ=0°和θ=60°處,如圖1中A和C點(diǎn)產(chǎn)生的奇異性,采用數(shù)學(xué)極限的方法,確定流動矢量。當(dāng)θ=0°時,當(dāng)θ=60°時,當(dāng)b≠1時,在點(diǎn)θ=0°和θ=60°處產(chǎn)生的奇異,采用物理的方法,確定流動矢量。當(dāng)θ=0°時,當(dāng)θ=60°時,3比例修正因子確定施加荷載增量后,首先計算試探應(yīng)力:將上述試探應(yīng)力帶入屈服條件(1),如果它不滿足,表明此時材料的行為是彈性的。硬化參數(shù)A保持不變,第r個增量步的應(yīng)力就等于試探應(yīng)力。如果滿足屈服條件,則表示在當(dāng)前荷載步內(nèi),積分點(diǎn)處達(dá)到塑性條件,因此,應(yīng)按照塑性規(guī)律進(jìn)行計算,即該點(diǎn)的應(yīng)力應(yīng)在屈服面上移動。分以下兩種情況:第1種,如圖2所示,應(yīng)力狀態(tài)由A點(diǎn)穿越屈服面達(dá)到B點(diǎn)。此時為保持應(yīng)力狀態(tài)在屈服面上,則C點(diǎn)的位置,可由下式?jīng)Q定:式中:R為比例修正因子第2種,如圖2所示,應(yīng)力狀態(tài)由屈服面A′點(diǎn)達(dá)到B′點(diǎn)。此時如果荷載增量比較大,且應(yīng)力點(diǎn)位于屈服面的大曲率附近,則上述過程仍存在較大的誤差。為了提高計算的精度,可以根據(jù)超出屈服面的應(yīng)力大小動態(tài)設(shè)定n個等分?jǐn)?shù),依次將迭代應(yīng)力拉回到屈服面,最后再采用上述的比例修正法進(jìn)行修正,就可得到較為準(zhǔn)確的結(jié)果。如圖3所示,將超過屈服面應(yīng)力點(diǎn)分成n等分,則經(jīng)過n個循環(huán)后,應(yīng)力點(diǎn)將返回偏離屈服面的E點(diǎn),然后再通過比例修正使其返回屈服面的E′點(diǎn)。4彈塑性本構(gòu)模型求解過程ABAQUS的用戶材料子程序(UMAT)通過與其求解器Standard的接口來實(shí)現(xiàn)數(shù)據(jù)交流。UMAT有自己的書寫格式與一些規(guī)范,與主程序共享的變量必須在子程序開頭予以定義。而主程序通過ABAQUS輸入文件(.inp)中的關(guān)鍵字“USERMATERIAL”來判斷用戶是否使用了自定義材料本構(gòu)模型,從而擴(kuò)展了它的適用性和應(yīng)用空間。統(tǒng)一彈塑性本構(gòu)模型主要求解過程:每一個增量加載步開始時,ABAQUS主程序在單元的積分點(diǎn)上調(diào)用UMAT子程序,并傳入應(yīng)變增量、時間步長及荷載增量,同時也傳入當(dāng)前已知狀態(tài)的應(yīng)力、應(yīng)變及其他與求解過程相關(guān)的變量;UMAT子程序根據(jù)本構(gòu)方程求解應(yīng)力增量并更新應(yīng)力及其他相關(guān)的量,提供Jacobian矩陣給ABAQUS主程序以形成整體剛度矩陣;主程序結(jié)合當(dāng)前荷載增量求解位移增量,繼而進(jìn)行平衡校核;如果不滿足指定的誤差,ABAQUS將進(jìn)行迭代直到認(rèn)為收斂,然后進(jìn)行下一增量步的求解。根據(jù)上述思路,本文的UMAT子程序主要流程如圖4所示。5正確計算格式和程序的一致性5.1關(guān)于統(tǒng)一強(qiáng)度理論及驗證驗證算例為一個四邊形平面應(yīng)變單元,試件長、寬均為100mm,如圖5所示。首先,采用ABAQUS自帶的Mohr-Coulomb材料,進(jìn)行軸向位移壓縮試驗;然后調(diào)用用戶自定義的統(tǒng)一強(qiáng)度理論子程序,令參數(shù)b=0,使其退化為Mohr-Coulomb強(qiáng)度屈服準(zhǔn)則,同樣進(jìn)行單軸壓縮試驗;最后將兩者的計算結(jié)果進(jìn)行對比。ABAQUS軟件中,應(yīng)力受壓為負(fù),受拉為正。材料參數(shù)見表1。圖6為整個計算過程中單元應(yīng)力σ2與軸向位移的關(guān)系曲線。由圖可見,兩者計算結(jié)果十分吻合。表明在ABAQUS中開發(fā)實(shí)現(xiàn)統(tǒng)一強(qiáng)度理論(b=0)退化為Mohr-Coulomb的計算結(jié)果是正確的。對于統(tǒng)一強(qiáng)度理論,當(dāng)b分別取0、0.5、1.0時,單軸試驗的單元應(yīng)力σ2與軸向位移的關(guān)系曲線如圖7所示。從圖中看出,屈服應(yīng)力隨b的增加而增加,這恰與圖1相吻合。當(dāng)考慮硬化參數(shù)A=1.0E7為常數(shù),b分別取0、0.5、1.0時,單軸試驗的單元應(yīng)力σ2與軸向位移的關(guān)系曲線如圖8所示。單元的σ2超過初始屈服應(yīng)力后,隨變形的增加線性增加。5.2圍巖應(yīng)力分析圖9為一圓形硐室的計算網(wǎng)格,隧道半徑為10m,地應(yīng)力p0在x和y方向均為30MPa。本文取1/4模型進(jìn)行計算。統(tǒng)一彈塑性本構(gòu)模型的計算參數(shù)如表2所示。理想彈塑性Mohr-Coulomb材料圍巖內(nèi)應(yīng)力的解析表達(dá)式、塑性區(qū)的徑向應(yīng)力σrp和環(huán)向應(yīng)力σθp分別為式中:a為圓形硐室半徑;r為圍巖內(nèi)一點(diǎn)至圓心的半徑。設(shè)σr1為彈塑性交界面r=r1處的徑向應(yīng)力,則彈性區(qū)域的徑向應(yīng)力σre和環(huán)向應(yīng)力σθe分別為塑性區(qū)半徑r1及彈塑性交界面上的徑向應(yīng)力為當(dāng)b=0時,統(tǒng)一強(qiáng)度理論退化為Mohr-Coulomb強(qiáng)度理論。如圖10所示,計算所得徑向應(yīng)力及環(huán)向應(yīng)力的數(shù)值解與解析解吻合很好。徑向應(yīng)力的最大相對誤差為1.55%,環(huán)向應(yīng)力的最大相對誤差為1.46%,塑性區(qū)半徑最大相對誤差為1.38%。圖11、12為統(tǒng)一強(qiáng)度理論當(dāng)b=0、0.5、1.0時的圍巖內(nèi)部徑向應(yīng)力與環(huán)向應(yīng)力分布曲線。從圖中可以看出,塑性區(qū)半徑外,徑向應(yīng)力和σr環(huán)向應(yīng)力σθ隨b值的增加而增加,但相差不大,且隨著半徑的增加,最終趨于穩(wěn)定值。塑性區(qū)半徑內(nèi),環(huán)向應(yīng)力明顯增加。從圖12可以看出,塑性區(qū)半徑r1隨b值的增加而減小,當(dāng)b=0時,r1=4.59;當(dāng)b=1時,r1=2.29。當(dāng)b=0和1時,考慮硬化參數(shù)A=5e8為常數(shù)與不考慮硬化的徑向應(yīng)力及環(huán)向應(yīng)力的數(shù)值解比較分別見圖13、14所示?？紤]硬化后,塑性區(qū)半徑有所減小,塑性區(qū)內(nèi)徑向應(yīng)力相差不大,環(huán)向應(yīng)力提高較大。6q-us在abaqus軟件中的應(yīng)用巖土材料是一種復(fù)雜的工程介質(zhì),采用現(xiàn)有通用有限元軟件ABAQUS所提供的本構(gòu)模型往往不能滿足工程實(shí)際數(shù)值分析的需要。為此,本文依據(jù)ABAQUS所提供的二次開發(fā)程序接口,結(jié)合統(tǒng)一彈塑性本構(gòu)模型,推導(dǎo)了其在ABAQ-US中的增量迭代格式,并編制了相應(yīng)的接口程序,可在ABAQUS軟件中進(jìn)行應(yīng)用。首先通過使計算參數(shù)b=0,將統(tǒng)一彈塑性本構(gòu)模型退化為Mohr-Coulomb模型,進(jìn)行了單軸試驗?zāi)M,并與軟件自帶的模型進(jìn)行了對比。結(jié)果表明,所開發(fā)的統(tǒng)一模型得到的結(jié)果與軟件自帶模型