基于無單元法的水力劈裂數(shù)值模擬計(jì)算

張 鑫

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

1. 計(jì)算方法

根據(jù)損傷力學(xué)的相關(guān)理論，材料受損傷前后彈性模量是變化的，并認(rèn)為無損材料的應(yīng)變與全應(yīng)力作用于受損材料產(chǎn)生的應(yīng)變是等價(jià)的。根據(jù)這一假定，引入損傷變量，將受損材料的彈性模量定義為

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

巖體等材料的拉應(yīng)力大于其抗拉強(qiáng)度時(shí)，材料的損傷變量為

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

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

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

基于以上原理建立應(yīng)力-滲流-損傷耦合計(jì)算模型，并編制相關(guān)計(jì)算程序。

2. 注意事項(xiàng)

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

3. 算例分析

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

圖1平面應(yīng)力模型示意圖

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

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

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

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

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