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

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

（1. 長安大學 公路學院，西安 710064； 2. 中交第一公路勘察設計研究院有限公司，西安 710065）

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

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

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

1 黏彈塑性流變模型

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

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

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

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

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

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

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

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

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

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

3)當σ≥σs且ε≥εc時，五元件模型（H-H/N-S/N），可以描述巖體衰減蠕變階段、穩態蠕變階段和加速蠕變階段的蠕變特性，其相應的狀態方程可以表示為

式中：σs為巖體蠕變的應力閾值（文中將巖體長期強度作為應力閾值）；εc為巖體進入加速蠕變的應變閾值；tc為與之對應的進入加速蠕變的時間；D為損傷變量[17]，其表達式為

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

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

巖體在實際工程中處于復雜的三維應力狀態，其應力張量可分解為球應力張量σm和偏應力張量Sij，其表達式分別為

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

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

根據式（1）~式（12），結合疊加原理，可得巖體三維應力狀態下的蠕變方程

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

式中，J2為第二偏應力不變量。

2 黏彈塑性模型數值程序二次開發

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

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

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

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

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

則3部分的黏塑性應變增量可表示為

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

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

2.2 黏彈塑性模型子程序的二次開發

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

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

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

式中：n是增量步；θ是0~1的積分參數[29]，θ=0對應一個顯式Euler積分算法，θ=1對應一個隱式的Euler積分算法。

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

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

2.3 黏彈塑性UMAT算例驗證

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

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

圖4 三軸壓縮蠕變試驗Fig. 4 Triaxial compression creep test

圖5給出了圓柱體試樣軸向應變蠕變曲線的黏彈塑性五元件模型的UMAT子程序數值計算結果，并與河海模型[10]及其流變試驗[26]結果進行了比較。可以發現，文中提出的黏彈塑性模型可以很好地描述巖體蠕變的3個階段，且基于黏彈塑性模型的UMAT子程序計算結果與河海模型及蠕變試驗結果基本一致，從而驗證了黏彈塑性模型合理性及UMAT子程序的有效性。

圖5 巖體流變模型與試驗結果的比較Fig. 5 Comparison of rock mass rheological model and test results

3 結 論

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

2）通過引入非定常牛頓黏壺改進了Burgers模型，黏壺同時兼顧時間相關性與巖體的微元漸進損傷，使五元件蠕變模型具備了描述減速蠕變、穩態蠕變和加速蠕變的特點。

3）在恒定應力狀態下根據一維黏彈塑性本構方程推導出三維黏彈塑性本構方程，并將其用于ABAQUS有限軟件UMAT子程序的研發，結合河海模型及其蠕變試驗驗證了該模型的合理性與有效性。

4）非定常黏彈塑性五元件模型可滿足巖體蠕變的各個階段，通過有限元軟件可將其用于隧道工程巖體蠕變的數值分析，為重大巖土體工程的結構設計及長期的運營維護提供理論依據，降低工程實踐過程中的潛在風險。