基于互層巖體的非定常黏彈塑性蠕變模型及有限元分析

胡濤濤，康志斌，陳建勛，胡雄，王棟

（1. 長安大學(xué) 公路學(xué)院，西安 710064； 2. 中交第一公路勘察設(shè)計(jì)研究院有限公司，西安 710065）

巖體的蠕變效應(yīng)是隧道工程、采礦工程及邊坡工程發(fā)生時(shí)效變形乃至失穩(wěn)破壞的重要原因之一[1-3]。對于軟弱巖體及含有夾層破碎帶的松散巖體，其蠕變現(xiàn)象尤為顯著[4]；在高地應(yīng)力狀態(tài)下，即便是中等強(qiáng)度的巖體或存在節(jié)理發(fā)育的硬巖也存在一定程度的流變效應(yīng)[5-6]；由此可見，巖體的蠕變現(xiàn)象在自然界是普遍存在的。巖體的蠕變特性對隧道及邊坡結(jié)構(gòu)的安全性和長期穩(wěn)定性有重大影響，建立能夠描述巖體全程蠕變特性的蠕變模型至關(guān)重要。然而，軟硬互層巖體的蠕變特性更為復(fù)雜，其在相同的應(yīng)力條件下軟巖與硬巖具有不同的變形特性（軟巖較硬巖更容易出現(xiàn)加速蠕變現(xiàn)象）[7]，為重大巖土體工程實(shí)踐帶來了嚴(yán)峻的挑戰(zhàn)，這就需要對巖體蠕變模型展開進(jìn)一步的研究。

針對巖體的蠕變特性，Griggs[8]最先對灰?guī)r、頁巖等軟巖進(jìn)行蠕變試驗(yàn)研究，發(fā)現(xiàn)當(dāng)加載于軟巖時(shí)的荷載達(dá)到破壞荷載的12.5%~80%，就會產(chǎn)生一定程度的蠕變現(xiàn)象；Fujii等[9]對花崗巖進(jìn)行了三軸蠕變試驗(yàn)，得到了軸向應(yīng)變、環(huán)向應(yīng)變和體積應(yīng)變等多種蠕變曲線，提出環(huán)向應(yīng)變可以作為判斷巖體損傷的重要指標(biāo)。近年來，隨著損傷斷裂力學(xué)的等新理論的引入，巖體加速蠕變模型成為了巖體流變力學(xué)研究的熱點(diǎn)。徐衛(wèi)亞等[10]將非線性黏塑性體與五元件黏彈性模型串聯(lián)建立了河海模型，充分反映了巖體的加速流變特性；王閆超等[11]以巴東組泥巖為研究對象，通過引入應(yīng)變閾值的非線性牛頓體元件，建立了一個(gè)新的八元件非線性蠕變本構(gòu)模型；張亮亮等[12]在經(jīng)典元件模型的基礎(chǔ)上引入非線性元件和蠕變損傷的方法，描述了巖體整個(gè)蠕變過程中的非線性蠕變特性；此外，大型數(shù)值仿真軟件在工程實(shí)踐中的應(yīng)用越來越廣泛，多被用來解決應(yīng)力條件復(fù)雜的特殊巖土體問題，有學(xué)者將復(fù)雜的流變模型通過二次開發(fā)集成于有限元軟件或有限差分軟件中[13-15]，為數(shù)值研究提供了更為豐富的基礎(chǔ)資料，對復(fù)雜地質(zhì)條件的工程實(shí)踐提供了指導(dǎo)作用。時(shí)至今日，國內(nèi)外學(xué)者已經(jīng)在巖石蠕變特性的理論研究及工程實(shí)踐中取得重大成果[16-22]，但仍然存在一些不足：其一，大量的研究仍然聚焦于完整巖體蠕變特性的研究，對于軟硬互層巖體蠕變特性的研究極少[23-24]；其二，巖體的非線性蠕變特性與巖體的應(yīng)力狀態(tài)和應(yīng)力作用時(shí)間密切相關(guān)，而忽略了巖體的應(yīng)變狀態(tài)對加速蠕變特性的影響，其合理性有待商榷。

為此，筆者在Burgers模型的基礎(chǔ)上提出了一種考慮應(yīng)力、應(yīng)變閾值的黏彈塑性蠕變模型，用以描述互層巖體的蠕變特性；并結(jié)合ABAQUS有限元軟件的UMAT接口子程序完成了該數(shù)值程序的二次開發(fā)，根據(jù)文獻(xiàn)[25-26]中的單軸壓縮試驗(yàn)驗(yàn)證了該程序的正確性與有效性。

1 黏彈塑性流變模型

1.1 黏彈塑性流變模型的特性

完整巖體尤其是軟弱巖體在給定的應(yīng)力條件下具有時(shí)效變形的特性，即巖體的流變特性。具有流變特性的巖體的蠕變分為減速蠕變Ⅰ（ab）、穩(wěn)態(tài)蠕變Ⅱ(bc）和加速蠕變Ⅲ(cd）3個(gè)階段，如圖1所示。丁秀麗等[7]表明軟硬互層巖體的流變特性明顯不同，軟巖蠕變變形包括減速蠕變、穩(wěn)態(tài)蠕變和加速蠕變3個(gè)階段，符合黏彈塑性本構(gòu)模型的蠕變特性；硬巖的蠕變變形僅存在減速蠕變一個(gè)階段，符合廣義開爾文模型的蠕變特性。此外，通過室內(nèi)試驗(yàn)研究發(fā)現(xiàn)[5,27]，在不同的應(yīng)力狀態(tài)下，巖樣由穩(wěn)態(tài)蠕變階段進(jìn)入加速蠕變階段是巖土體材料內(nèi)部微元應(yīng)變損傷累積的結(jié)果；其宏觀表現(xiàn)為軸向應(yīng)變達(dá)到某一數(shù)值εc時(shí)，巖樣進(jìn)入加速蠕變階段。相較于單一的完整巖體，軟硬互層巖體對元件組合模型的要求更高。

圖1 巖體蠕變特性曲線Fig. 1 Creep characteristic curve of rock mass

對此，文中提出了一種考慮應(yīng)力與應(yīng)變閾值的黏彈塑性非定常元件組合模型，如圖2所示。

圖2 五元件黏彈塑性模型Fig. 2 Five-element viscoelastoplastic model

圖2中模型可以充分反映軟硬互層巖體蠕變變形的相應(yīng)階段，且該組合模型應(yīng)滿足如下3個(gè)條件。

1）當(dāng)σ<σs時(shí)，1、2部分所有元件參與蠕變工作，組合模型為三元件模型（H-H/N），可以描述巖體的衰減蠕變階段的特性，與之相應(yīng)的狀態(tài)方程[1]可以表示為

根據(jù)式（1）得三元件模型的一維蠕變方程為

2）當(dāng)σ≥σs且ε<εc時(shí)，1、2、3部分所有元件參與蠕變工作，組合模型為五元件模型（H-H/N-S/N），可以描述巖體衰減蠕變階段和穩(wěn)態(tài)蠕變階段的蠕變特性，其相應(yīng)的狀態(tài)方程可以表示為

根據(jù)式（3）得五元件模型的一維蠕變方程為

3)當(dāng)σ≥σs且ε≥εc時(shí)，五元件模型（H-H/N-S/N），可以描述巖體衰減蠕變階段、穩(wěn)態(tài)蠕變階段和加速蠕變階段的蠕變特性，其相應(yīng)的狀態(tài)方程可以表示為

式中：σs為巖體蠕變的應(yīng)力閾值（文中將巖體長期強(qiáng)度作為應(yīng)力閾值）；εc為巖體進(jìn)入加速蠕變的應(yīng)變閾值；tc為與之對應(yīng)的進(jìn)入加速蠕變的時(shí)間；D為損傷變量[17]，其表達(dá)式為

根據(jù)式（3）得五元件模型的一維蠕變方程為

1.2 黏彈塑性模型蠕變的三維形式

巖體在實(shí)際工程中處于復(fù)雜的三維應(yīng)力狀態(tài)，其應(yīng)力張量可分解為球應(yīng)力張量σm和偏應(yīng)力張量Sij，其表達(dá)式分別為

球應(yīng)力張量σm能引起物體體積的改變，但不能改變其形狀；而偏應(yīng)力張量Sij只能引起形狀改變，但不能改變其體積。相應(yīng)地，可將巖體的應(yīng)變張量分解為球應(yīng)變張量εm和偏應(yīng)變張量εij，其表達(dá)式分別為

一維蠕變方程中巖體的彈性模量E和泊松比ν在三維蠕變方程中應(yīng)采用與之相應(yīng)的體積模量G與剪切模量K，可分別表示為

根據(jù)式（1）~式（12），結(jié)合疊加原理，可得巖體三維應(yīng)力狀態(tài)下的蠕變方程

式中：σs為巖體穩(wěn)態(tài)蠕變的長期強(qiáng)度，由室內(nèi)試驗(yàn)確定；β為非線性牛頓體的元件參數(shù)，均由試驗(yàn)確定；H(σ)為單位躍階函數(shù)；F為巖土體材料的屈服函數(shù)[28]；Q為塑性勢函數(shù)，且采用相關(guān)聯(lián)流動法則，其表達(dá)式分別為

式中，J2為第二偏應(yīng)力不變量。

2 黏彈塑性模型數(shù)值程序二次開發(fā)

2.1 黏彈塑性模型的有限元分析方法

在有限元法的計(jì)算過程中，根據(jù)前文給出的一維及三維形式的黏彈塑性蠕變方程，需將其表示為相應(yīng)的增量方程形式?？倯?yīng)變增量Δε包括彈性應(yīng)變增量Δεe、黏彈性應(yīng)變增量Δεev和黏塑性應(yīng)變增量Δεvp三部分為

對圖2中1部分彈性元件有：

對圖2中2部分黏彈性元件有：

對圖2中3部分黏塑性元件有：

則3部分的黏塑性應(yīng)變增量可表示為

記蠕變增量為Δεc，則有：

式中：D為彈性剛度矩陣；Δσ為應(yīng)力增量矩陣；Δε應(yīng)變增量矩陣。

2.2 黏彈塑性模型子程序的二次開發(fā)

ABAQUS有限元軟件為用戶提供了子程序二次開發(fā)接口，用戶可以使用Fortran語言編寫程序，自定義材料子程序（user material subroutine, UMS），通過Standard求解器的接口實(shí)現(xiàn)與主程序之間的數(shù)據(jù)交流。按照ABAQUS非線性增量加載和Newton平衡迭代算法（增量迭代法），軟件在每一個(gè)增量步對每一個(gè)計(jì)算單元均調(diào)用UMAT，獲得Jacobian矩陣DDSDDE（即剛度系數(shù)矩陣D），然后通過Standard接口傳遞給ABAQUS主程序，在主程序完成平衡校核；如果平衡校核不滿足誤差要求，ABAQUS將繼續(xù)進(jìn)行迭代，直至滿足誤差要求，再進(jìn)入以下增量步的求解?？梢奤MAT會被頻繁調(diào)用，因此要充分考慮子程序代碼的質(zhì)量，提高計(jì)算效率。

筆者提出的黏彈塑性模型用于解決與時(shí)間相關(guān)的非線性問題，其剛度矩陣隨時(shí)間不斷發(fā)生變化，每個(gè)增量步調(diào)用UMAT均需要重新計(jì)算剛度矩陣，計(jì)算效率大大降低。為此，采用常剛度法進(jìn)行計(jì)算，在時(shí)間增量Δtn=tn+1-tn內(nèi)通過應(yīng)變增量更新子程序的應(yīng)力增量和狀態(tài)變量，根據(jù)式（18）~（26）可得：

為提高計(jì)算精度，Δεcn采用差分法計(jì)算，即：

式中：n是增量步；θ是0~1的積分參數(shù)[29]，θ=0對應(yīng)一個(gè)顯式Euler積分算法，θ=1對應(yīng)一個(gè)隱式的Euler積分算法。

黏彈塑性本構(gòu)模型通過UMAT在ABAQUS中實(shí)現(xiàn)，每一個(gè)增量步開始ABAQUS 主程序在單元的積分點(diǎn)傳入應(yīng)變增量、時(shí)間增量步，同時(shí)也傳入當(dāng)前已知狀態(tài)的應(yīng)力、應(yīng)變及其他與求解過程相關(guān)的變量；UMAT根據(jù)本構(gòu)方程計(jì)算應(yīng)力增量并更新應(yīng)力及狀態(tài)變量，提供 Jacobian 矩陣DDSDDE給 ABAQUS 主程序以形成整體剛度矩陣；主程序根據(jù)當(dāng)前荷載增量求解位移增量，完成平衡校核；如果平衡校核不滿足誤差要求，ABAQUS 將進(jìn)行迭代直至收斂（筆者判定收斂的標(biāo)準(zhǔn)為平衡校核相對誤差小于10-6），然后進(jìn)行下一增量步。綜上，該模型的UMAT主要流程如圖3所示。

圖3 黏彈塑性模型UMAT子程序流程圖Fig. 3 Sticky elastoplastic model UMAT subprogram flowchart

2.3 黏彈塑性UMAT算例驗(yàn)證

為了驗(yàn)證文中提出的黏彈塑性五元件模型的正確性與合理性，本算例建立了直徑35.5 mm，高70 mm的圓柱體試樣模型，且與文獻(xiàn)[25-26]三軸流變試驗(yàn)條件保持一致，圓柱體試樣底部固定（Ux=Uy=Uz=0），試樣頂部施加100 MPa的均布荷載，試樣周圍施加15 MPa的圍壓，如圖4（a）所示。通過文中開發(fā)的黏彈塑性UMAT子程序，在ABAQUS中完成數(shù)值計(jì)算，軸向應(yīng)變云圖，如圖4（b）所示，巖體流變力學(xué)參數(shù)見表1所示。相比于河海模型（七元件模型），文中模型（五元件模型）組合元件更少，使得2種模型的個(gè)別參數(shù)存在較大的差異。

表1 蠕變模型參數(shù)[10, 26]Table 1 Creep model parameters

圖4 三軸壓縮蠕變試驗(yàn)Fig. 4 Triaxial compression creep test

圖5給出了圓柱體試樣軸向應(yīng)變?nèi)渥兦€的黏彈塑性五元件模型的UMAT子程序數(shù)值計(jì)算結(jié)果，并與河海模型[10]及其流變試驗(yàn)[26]結(jié)果進(jìn)行了比較?？梢园l(fā)現(xiàn)，文中提出的黏彈塑性模型可以很好地描述巖體蠕變的3個(gè)階段，且基于黏彈塑性模型的UMAT子程序計(jì)算結(jié)果與河海模型及蠕變試驗(yàn)結(jié)果基本一致，從而驗(yàn)證了黏彈塑性模型合理性及UMAT子程序的有效性。

圖5 巖體流變模型與試驗(yàn)結(jié)果的比較Fig. 5 Comparison of rock mass rheological model and test results

3 結(jié) 論

1）根據(jù)軟巖互層巖體的蠕變特性，將應(yīng)力、應(yīng)變閾值引入黏彈塑性蠕變模型，建立了五元件非線性蠕變方程，使其能夠同時(shí)描述硬巖的線性黏塑性蠕變特性與軟巖的非線性黏塑性蠕變特性。

2）通過引入非定常牛頓黏壺改進(jìn)了Burgers模型，黏壺同時(shí)兼顧時(shí)間相關(guān)性與巖體的微元漸進(jìn)損傷，使五元件蠕變模型具備了描述減速蠕變、穩(wěn)態(tài)蠕變和加速蠕變的特點(diǎn)。

3）在恒定應(yīng)力狀態(tài)下根據(jù)一維黏彈塑性本構(gòu)方程推導(dǎo)出三維黏彈塑性本構(gòu)方程，并將其用于ABAQUS有限軟件UMAT子程序的研發(fā)，結(jié)合河海模型及其蠕變試驗(yàn)驗(yàn)證了該模型的合理性與有效性。

4）非定常黏彈塑性五元件模型可滿足巖體蠕變的各個(gè)階段，通過有限元軟件可將其用于隧道工程巖體蠕變的數(shù)值分析，為重大巖土體工程的結(jié)構(gòu)設(shè)計(jì)及長期的運(yùn)營維護(hù)提供理論依據(jù)，降低工程實(shí)踐過程中的潛在風(fēng)險(xiǎn)。