張 鑫

近年來，諸如無單元法、邊界元法、非連續(xù)變形分析（DDA）等更適用于解決大變形及裂紋擴展等問題的數(shù)值計算方法的快速發(fā)展，為更好地解決水力劈裂等問題提供了有力的數(shù)值計算工具。本文利用損傷力學的相關概念，并基于無單元法建立起應力-滲流-損傷三場耦合數(shù)值計算模型，充分發(fā)揮了無單元法在追蹤裂紋擴展方面的優(yōu)勢。的關系為：巖體處于彈性狀態(tài)時，材料未發(fā)生破壞，損傷變量D=0，滲透系數(shù)和有效應力為負指數(shù)關系。當應力大于材料的抗拉強度時，裂紋開裂，材料發(fā)生破壞，此時 0＜D＜1，引入滲透系數(shù)增大系數(shù)，材料滲透性有所提高，當材料發(fā)生破壞時，D=1，滲透系數(shù)將會大幅增大。

1. 計算方法

根據(jù)損傷力學的相關理論，材料受損傷前后彈性模量是變化的，并認為無損材料的應變與全應力作用于受損材料產生的應變是等價的。根據(jù)這一假定，引入損傷變量，將受損材料的彈性模量定義為

式中E，為產生裂紋材料的彈性模量；E為未產生裂紋材料的彈性模量。D為損傷變量，D=0表示材料處于無損傷狀態(tài)，D=1表示材料處于完全損傷（斷裂或者破壞）狀態(tài)，0＜D＜1對應不同的損傷程度。

巖體等材料的拉應力大于其抗拉強度時，材料的損傷變量為

式中，σcr為材料損傷后的殘余強度；E意義同上。

巖體材料的滲透系數(shù)和其應力狀態(tài)存在一定的關聯(lián)，主要表現(xiàn)在隨著巖體應力的增大，當其裂尖材料的應力大于抗拉強度時，巖體開始破壞，裂紋開始擴展，裂紋的開度隨之增加，裂隙間材料的滲透性增加，相應結點的滲透系數(shù)將會增大；相反，若裂紋閉合或其未擴展，材料未損傷，相應區(qū)域材料的滲透系數(shù)不增大，其數(shù)值和有效應力仍為負指數(shù)關系。根據(jù)以上原理，建立應力和滲透系數(shù)

式中，ξ為滲透系數(shù)增大系數(shù)。

基于以上原理建立應力-滲流-損傷耦合計算模型，并編制相關計算程序。

2. 注意事項

滲流場、應力場以及損傷場3場耦合關系十分復雜，需要各自考慮的邊界條件也不盡相同，得到解析解是十分困難的。本文分別計算3場，并利用3場之間的相互耦合關系進行迭代循環(huán)，以達到它們的動態(tài)平衡，以滿足計算精度作為迭代計算終止條件。

3. 算例分析

如圖1所示，建立巖體二維平面應力模型，邊界條件為：水平方向兩側圍壓4MPa，豎直方向頂端水壓力 2.3MPa（向下），底端水壓力 3.8MPa（向上）；巖體的彈性模量10GPa，巖體滲透系數(shù)為 1.0×10-6m/s，抗拉強度 10MPa，殘余抗拉強度0.1MPa，抗壓強度 100MPa，殘余抗壓強度10MPa。巖體二維計算模型和無單元模型高斯點分布如圖1所示，其中初始裂紋與底邊的夾角為58°。分別對不考慮損傷和考慮損傷2種情況進行對比計算分析。

圖1平面應力模型示意圖

圖2不考慮損傷變量流速矢量分布圖

圖3考慮損傷變量流速矢量分布圖

通過對比圖2和圖3考慮及不考慮損傷時流速矢量分布圖可以看出，考慮損傷作用后，在裂紋擴展的部分區(qū)域的損傷變量達到了材料的抗拉和抗壓強度，巖體發(fā)生了破壞，裂紋體內部裂隙的連通性有所提高，使得裂紋計算區(qū)域的滲透系數(shù)有了較大幅度的增加，相關區(qū)域內的滲透流速矢量較之未考慮損傷作用時亦有所增大。此外，在考慮損傷作用后，沿裂紋走向以及裂尖的部分區(qū)域的應力值較未考慮損傷作用時應力值偏小。這主要是由于在迭代計算過程中，裂紋計算區(qū)域的應力達到了材料的損傷闕值，材料產生了不同程度的損傷，相關高斯計算點附近的材料強度按照彈脆性損傷本構方程式進行了剛度退化處理，故其應力有所降低。

在模擬水力劈裂的數(shù)值模型中加入損傷變量，考慮其對于巖體滲透系數(shù)和巖體剛度的影響后得出的結論。

無單元法作為新興的計算方法，擺脫了網(wǎng)格劃分的局限，在水力劈裂數(shù)值模擬計算過程中可以隨著裂紋的擴展在裂尖位置隨意增加結點，較傳統(tǒng)的有限元法其在處理大變形及裂紋擴展等問題方面具有明顯的優(yōu)勢。無單元法模擬分析了含有裂縫的巖體的破裂過程。通過對比計算表明，考慮了損傷變量對于水力劈裂過程中的滲流場特性和應力分布影響的計算結果與文獻的實驗結果更為接近。因此，耦合模型可以用來研究巖體的破裂或水力劈裂，具有理論意義和應用價值。