王 慶，姚 俊，譚文祿

(中電建生態(tài)環(huán)境集團(tuán)有限公司，廣東 深圳 518102)

在邊坡[1-2]、壩基[3]、隧洞[4]等巖石工程中，經(jīng)常會(huì)遇到地下水與巖石(體)之間的耦合問(wèn)題，其作用機(jī)理是：一方面，當(dāng)巖石(體)變形時(shí)會(huì)引起孔隙率、節(jié)理厚度等變化，從而引起滲流場(chǎng)的變化；另一方面滲流場(chǎng)變化會(huì)引起孔隙水壓力和有效應(yīng)力的變化，從而引起巖石(體)應(yīng)力應(yīng)變場(chǎng)的變化。對(duì)應(yīng)力滲流耦合最早的研究是Biot的三維固結(jié)理論[5]，但該理論中滲透系數(shù)在變形過(guò)程中為常數(shù)。有關(guān)學(xué)者的研究表明[6-9]，滲透系數(shù)在巖石的變形破壞過(guò)程中，會(huì)發(fā)生較大的變化，可能會(huì)影響工程的安全性，法國(guó)Malpasset拱壩失事正是沒(méi)有考慮到滲透系數(shù)變化而導(dǎo)致的[10]，因此有必要研究巖石變形過(guò)程中滲透系數(shù)的變化規(guī)律。

目前應(yīng)力滲流耦合的研究方法主要分為宏觀唯象和細(xì)觀2類(lèi)方法。宏觀唯象方法又可以分為直接法和間接法。直接法是指通過(guò)室內(nèi)室外實(shí)驗(yàn)對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合得到經(jīng)驗(yàn)公式的方法[11-12]，這些公式只能反映巖石變形過(guò)程中滲透系數(shù)的變化的某一階段，不能體現(xiàn)應(yīng)力峰值后的變化規(guī)律。為了克服經(jīng)驗(yàn)公式的這一缺點(diǎn)，學(xué)者們開(kāi)始研究全應(yīng)力-應(yīng)變過(guò)程中的滲透系數(shù)變化規(guī)律[13-16]，但由于公式的參數(shù)太多，形式受巖石性質(zhì)、實(shí)驗(yàn)條件等影響太大，而且不能解釋滲透系數(shù)變化的物理力學(xué)機(jī)理，因此不利于在工程中應(yīng)用。間接法則是通過(guò)引入中間變量，間接建立應(yīng)力(應(yīng)變)與滲透系數(shù)之間的關(guān)系，最常用的中間變量是孔隙率[17]。間接法雖然能減少模型參數(shù)數(shù)量，但巖石應(yīng)力應(yīng)變?nèi)^(guò)程不同階段的滲透系數(shù)變化的機(jī)理不同，如峰前滲透系數(shù)急劇增大的主要原因是微裂紋擴(kuò)展，而不是孔隙率變化，因此用一個(gè)中間變量很難描述全過(guò)程的滲透系數(shù)變化規(guī)律。細(xì)觀方法是通過(guò)描述巖石在變形過(guò)程中細(xì)觀結(jié)構(gòu)(孔隙、微裂紋等)的變化，從而引起滲透系數(shù)變化的方法[18-20]。細(xì)觀方法可以解釋巖石變形過(guò)程中的滲透系數(shù)變化機(jī)理，但計(jì)算過(guò)程中需要記錄每一個(gè)微裂紋的傾向、傾角、開(kāi)度等參數(shù)，因此需要借助計(jì)算機(jī)程序?qū)崿F(xiàn)。

為了提出既能解釋滲透系數(shù)的變化機(jī)理，又方便應(yīng)用的模型，許多學(xué)者開(kāi)始從等效連續(xù)介質(zhì)的角度將巖石細(xì)觀結(jié)構(gòu)的變化等效為損傷演化[21-23]。將損傷作為連接應(yīng)力(應(yīng)變)和滲透系數(shù)之間的橋梁，建立巖石、混凝土、水泥砂漿等類(lèi)巖石脆性材料的彈性損傷-滲流耦合模型，這個(gè)耦合過(guò)程中有2個(gè)重要的準(zhǔn)則，一個(gè)是損傷的演化準(zhǔn)則，另一個(gè)則是損傷過(guò)程中的滲透系數(shù)變化準(zhǔn)則，本文在傳統(tǒng)的彈性損傷模型和滲透性變化模型基礎(chǔ)上，提出改進(jìn)的彈性損傷-滲流耦合模型。

1 改進(jìn)的彈性損傷模型

傳統(tǒng)的巖石彈塑性損傷模型，雖然有可以考慮塑性和可以區(qū)分拉壓破壞的優(yōu)點(diǎn)，但是該本構(gòu)關(guān)系也有一些缺點(diǎn)。首先是輸入?yún)?shù)過(guò)多，彈塑性損傷本構(gòu)關(guān)系中拉壓曲線(xiàn)是分開(kāi)的，而且要輸入塑性強(qiáng)化的應(yīng)力～應(yīng)變曲線(xiàn)以及損傷因子隨著開(kāi)裂應(yīng)變的變化曲線(xiàn)，這些曲線(xiàn)的參數(shù)過(guò)多；另外一點(diǎn)是使用彈塑性損傷本構(gòu)關(guān)系計(jì)算，結(jié)果很難收斂，而且收斂性也不能控制。鑒于這兩個(gè)缺點(diǎn)，這里在Mazars模型的基礎(chǔ)上，改進(jìn)其收斂性和網(wǎng)格依賴(lài)性，提出一種改進(jìn)的彈性損傷模型。Mazars模型是一種描述混凝土損傷的本構(gòu)關(guān)系，在該模型中，應(yīng)力應(yīng)變關(guān)系如下：

(1)

(2)

(3)

(4)

其中，ABS()是取絕對(duì)值函數(shù)，如果f為正，結(jié)果為f；如果為負(fù)，結(jié)果為0。

在Mazars模型中，等效應(yīng)變作為內(nèi)變量來(lái)判斷損傷的初始化和演化，但是這里是將有效應(yīng)力分解為拉壓部分，按照這種思想，將拉壓有效應(yīng)力的等效應(yīng)力作為損傷初始化準(zhǔn)則，形式如下：

(5)

(6)

(7)

(8)

(9)

式中：

(10)

(11)

(12)

(13)

式(6)、(8)是受拉損傷的初始化準(zhǔn)則，其形式和最大拉應(yīng)力準(zhǔn)則類(lèi)似，式(7)、(9)是受壓損傷的初始化準(zhǔn)則，由于受壓破壞中含有壓縮和剪切破壞的組合，其形式類(lèi)似于Drucker-Prager屈服準(zhǔn)則，K是Drucker-Prager屈服準(zhǔn)則的錐角，R0是雙軸壓縮強(qiáng)度和單軸壓縮強(qiáng)度比，f±是拉伸和壓縮強(qiáng)度。

Mazars模型中的損傷演化準(zhǔn)則形式如下：

(14)

(15)

材料的軟化階段會(huì)導(dǎo)致邊值問(wèn)題的不適定，式(14)、(15)中的損傷演化準(zhǔn)則對(duì)網(wǎng)格的依賴(lài)性很大，為了減小網(wǎng)格依賴(lài)性，這里采用斷裂能的手段[25]來(lái)修正損傷演化準(zhǔn)則。常數(shù)A+和B-可以寫(xiě)成如下斷裂能的形式：

(16)

(17)

(18)

式中G±——拉壓斷裂能；lc——微裂紋的特征長(zhǎng)度；Ve——單元體積；E0——初始彈性模量。

考慮拉損傷和壓損傷的總的損傷因子定義如下：

d=1-(1-d+)(1-d-)

(19)

以上的改進(jìn)策略可以減小網(wǎng)格依賴(lài)性，但是材料軟化在隱式分析中，如ABAQUS Standard, 還存在收斂性問(wèn)題，這里采用正則化準(zhǔn)則來(lái)處理收斂性問(wèn)題。正則化準(zhǔn)則的主要目的是使軟化材料的切線(xiàn)模量在增量步足夠小的時(shí)候保持正定。根據(jù)式(1)中的應(yīng)力-應(yīng)變關(guān)系，可以得到Jacobian矩陣的形式如下:

(20)

為了提高收斂性，ABAQUS子程序UMAT中使用一種損傷因子的黏性正則化準(zhǔn)則，此時(shí)，不直接使用式(14)、(15)來(lái)計(jì)算損傷因子，而是使用如下形式的正則化損傷因子：

(21)

(22)

(23)

(24)

從式(23)和式(24)中可以看出：

(25)

因此，使用正則化損傷因子計(jì)算得到的Jacobian矩陣形式如下：

(26)

需要注意黏性系數(shù)η的選擇，過(guò)小的黏性系數(shù)導(dǎo)致收斂所需要的增量步很??；過(guò)大的黏性系數(shù)可能會(huì)過(guò)度延遲剛度的弱化，導(dǎo)致軟化階段的應(yīng)力-應(yīng)變曲線(xiàn)偏離實(shí)際過(guò)多。

2 損傷過(guò)程中滲透系數(shù)的變化

巖石的壓縮過(guò)程中滲透性試驗(yàn)中反應(yīng)出滲透系數(shù)的變化通常具有幾個(gè)階段。根據(jù)Zhang[26]的試驗(yàn)數(shù)據(jù)描繪的典型滲透系數(shù)變化的幾個(gè)階段見(jiàn)圖1。A是應(yīng)力-應(yīng)變曲線(xiàn)中損傷開(kāi)始的點(diǎn)，對(duì)應(yīng)于滲透系數(shù)開(kāi)始快速增大的點(diǎn)，B點(diǎn)是應(yīng)力-應(yīng)變曲線(xiàn)中的峰值應(yīng)力點(diǎn)，對(duì)應(yīng)于滲透系數(shù)曲線(xiàn)的轉(zhuǎn)折點(diǎn)，B點(diǎn)之后滲透系數(shù)增長(zhǎng)速度有所降低，然后達(dá)到滲透系數(shù)最大的點(diǎn)C，對(duì)應(yīng)于應(yīng)力-應(yīng)變曲線(xiàn)的轉(zhuǎn)折點(diǎn)，C點(diǎn)之后應(yīng)力下降速度減緩。OA階段滲透系數(shù)變化不大，但是有時(shí)候在開(kāi)始的時(shí)候滲透系數(shù)有所減小，原因是在壓應(yīng)力作用下孔隙變化和初始微裂紋的閉合。AC階段滲透系數(shù)的增大主要是因?yàn)槲⒘鸭y數(shù)目的增多以及微裂紋的張開(kāi)，C點(diǎn)之后滲透系數(shù)的減小主要是因?yàn)樵嚰茐闹笤趪鷫旱淖饔孟潞暧^裂紋閉合造成的。

圖1 脆性巖石的典型應(yīng)力-應(yīng)變曲線(xiàn)和滲透系數(shù)-應(yīng)變曲線(xiàn)

目前描述脆性巖石壓縮過(guò)程中滲透系數(shù)的變化的經(jīng)驗(yàn)?zāi)Ｐ秃屠碚撃Ｐ椭荒苊枋鰸B透系數(shù)變化的某個(gè)階段，通常是增大的階段，并不能同時(shí)描述滲透系數(shù)變化的整個(gè)過(guò)程。假設(shè)巖石在壓縮過(guò)程中的滲透系數(shù)的變化主要由孔隙閉合引起的滲透系數(shù)減小和局部損傷引起的滲透系數(shù)增大這2個(gè)因素造成的。描述孔隙閉合引起的滲透系數(shù)減小現(xiàn)象的理論有Kozeny-Carman模型，形式如下：

(27)

式中Kp——孔隙閉合引起的滲透系數(shù)變化；c——一個(gè)與細(xì)觀結(jié)構(gòu)有關(guān)的常數(shù)；n——孔隙度。

假設(shè)巖石顆粒不可壓縮，Vp是孔隙的體積；Vs是巖石顆粒的體積；e0是初始孔隙比；Δe是孔隙比的變化，則體積應(yīng)變?yōu)椋?/p>

(28)

考慮到孔隙比和孔隙率的關(guān)系e=n/(1-n)，代入式(28)可以得到孔隙率和體積應(yīng)變的關(guān)系：

(29)

式中n0——初始孔隙率。

將式(29)代入式(27)中，得到滲透系數(shù)隨著體積應(yīng)變變化關(guān)系為：

(30)

(31)

在描述損傷引起的滲透系數(shù)增大方面，本文將從局部損傷的角度推導(dǎo)出滲透系數(shù)的變化，可以看出在精度要求不高的時(shí)候，理論公式和經(jīng)驗(yàn)公式的形式類(lèi)似。在細(xì)觀結(jié)構(gòu)中存在著初始裂紋和后來(lái)發(fā)展的微裂紋，細(xì)觀結(jié)構(gòu)中的滲流主要受到這些微裂紋的影響，微裂紋中的滲流滿(mǎn)足平行板定理，其滲透系數(shù)為：

(32)

式中u——微裂紋的張開(kāi)度；ξ——微裂紋的粗糙度。

雖然式(32)比較簡(jiǎn)單，但是需要知道微裂紋的幾何產(chǎn)狀，而且本文的計(jì)算方法是基于連續(xù)介質(zhì)的，所以按照微裂紋的思想，用損傷帶代替微裂紋來(lái)描述滲流通道的改變，見(jiàn)圖2。一個(gè)面積為S的截面受到單位水頭梯度的作用，材料沒(méi)有損傷部分的滲透系數(shù)為K0，材料的局部損傷帶的寬度為λlc，在具體計(jì)算過(guò)程中可以看作是單元的特征長(zhǎng)度。等效的思想是讓損傷帶和微裂紋的流量相等，具體如下：

(33)

a) 微裂紋 b) 損傷帶

由式(33)可以得到損傷帶的滲透系數(shù)為：

(34)

ξ代表了微觀滲流通道的彎曲程度。最重要的是微裂紋開(kāi)度u的計(jì)算，本文通過(guò)改進(jìn)的彈性損傷模型來(lái)得到裂紋開(kāi)度。在擴(kuò)展有限元中，裂紋的張開(kāi)是通過(guò)位移跳躍來(lái)表示的，而在漸進(jìn)損傷理論中，裂紋的張開(kāi)是通過(guò)損傷帶的法向累積位移表示的，所以這里只考慮拉損傷造成的滲透系數(shù)變化，在不考慮卸載造成的裂紋閉合的情況下，裂紋的張開(kāi)可以表示為：

(35)

假設(shè)損傷帶內(nèi)的應(yīng)變是均勻的，式(35)的積分結(jié)果為:

(36)

要得到裂紋張開(kāi)度和損傷的關(guān)系，需要求得式(14)的反函數(shù)的解，這里反函數(shù)沒(méi)有顯式解，需要通過(guò)迭代等數(shù)值方法得到，假設(shè)反函數(shù)的解已經(jīng)得到，形式如下:

(37)

將式(37)代入到式(36)中得到微裂紋的開(kāi)度：

(38)

將式(38)代入到式(34)中就可以得到損傷帶的滲透系數(shù)，但因?yàn)槭?37)沒(méi)有顯式解，所以需要簡(jiǎn)化。之前提到過(guò)理論推導(dǎo)的損傷帶滲透系數(shù)的公式在計(jì)算精度要求不高的時(shí)候，和指數(shù)型經(jīng)驗(yàn)公式有類(lèi)似的形式，通過(guò)這一關(guān)系可以將損傷帶滲透系數(shù)簡(jiǎn)化為指數(shù)型公式。將式(14)的指數(shù)項(xiàng)展開(kāi)成泰勒級(jí)數(shù)，并取一次項(xiàng)，得到反函數(shù)的簡(jiǎn)化解為：

(39)

至此，由孔隙閉合引起的滲透系數(shù)減小和局部損傷引起的滲透系數(shù)增大都考慮到了。所以，能同時(shí)考慮滲透系數(shù)變化不同階段的模型可以通過(guò)結(jié)合Kozeny-Carman公式和指數(shù)型公式得到，形式如下：

K=Kp(1+αdβ)

(40)

式中，α、β是與微裂紋幾何產(chǎn)狀有關(guān)的常數(shù)。將式(30)代入式(40)中得到損傷過(guò)程中滲透系數(shù)變化公式為：

(1+αdβ)

(41)

3 數(shù)值算例

3.1 模型邊界和參數(shù)

選取水泥砂漿材料作為研究對(duì)象，具體成分是砂子、水泥和水，其水灰比為0.4，砂子重量占40%，水泥漿和砂子的密度分別為2 000、2 650 kg/m3，在不考慮氣孔的情況下，細(xì)觀結(jié)構(gòu)中基質(zhì)和砂子的體積比分別為67%、33%。砂子采用的粒徑小于0.5 mm的河砂，假設(shè)砂子的粒徑滿(mǎn)足在區(qū)間[0.3,0.5]之間的均勻分布。有限元模型大小為4×4×4(mm3)，位移邊界條件為在模型下邊界施加固定約束；由于需要研究峰值后的應(yīng)力應(yīng)變規(guī)律，故采用位移邊界來(lái)代替荷載邊界，在模型上邊界勻速施加0.04 mm的位移，對(duì)應(yīng)的最大豎向應(yīng)變?yōu)?.01；滲流邊界為在模型上邊界施加1 kPa恒定不變的孔隙水壓力，下邊界為自由透水邊界，其他邊界為不透水邊界，見(jiàn)圖3?；|(zhì)和顆粒的力學(xué)和滲流參數(shù)見(jiàn)表1、2。

圖3 水泥砂漿的細(xì)觀模型

表1 細(xì)觀成分的力學(xué)參數(shù)

表2 細(xì)觀成分的滲流參數(shù)

3.2 結(jié)果分析

由于水泥砂漿材料的細(xì)觀模型中含有基質(zhì)(水泥漿)和顆粒(砂子)2種組分，分析水泥砂漿的應(yīng)力應(yīng)變和滲流特性必須綜合基質(zhì)和顆粒的應(yīng)力場(chǎng)和滲流場(chǎng)，因此結(jié)果分析中將采用等效宏觀應(yīng)力和等效宏觀滲透系數(shù)來(lái)表示水泥砂漿的特性。等效宏觀應(yīng)力是細(xì)觀模型應(yīng)力場(chǎng)的體積平均值：

(42)

等效宏觀滲透系數(shù)的處理方式與等效宏觀應(yīng)力有所不同，首先根據(jù)細(xì)觀模型的流速場(chǎng)計(jì)算等效宏觀流速：

(43)

然后根據(jù)達(dá)西定律求解等效宏觀滲透系數(shù)：

(44)

單軸壓縮應(yīng)力-應(yīng)變曲線(xiàn)和滲透系數(shù)變化曲線(xiàn)見(jiàn)圖4。可以看出宏觀應(yīng)力-應(yīng)變曲線(xiàn)和室內(nèi)實(shí)驗(yàn)數(shù)據(jù)吻合得很好，室內(nèi)實(shí)驗(yàn)中的試件由于初始裂紋的存在，導(dǎo)致在彈性階段初有一個(gè)壓密階段，峰后階段改進(jìn)的損傷模型軟化曲線(xiàn)比較平滑。由于滲透實(shí)驗(yàn)數(shù)據(jù)的缺乏，所以這里不進(jìn)行滲透系數(shù)曲線(xiàn)的對(duì)比，但可以看出滲透系數(shù)的變化曲線(xiàn)與Zhang[26]的試驗(yàn)數(shù)據(jù)一致，可以分為3個(gè)階段，但應(yīng)力峰值后的滲透系數(shù)減小階段比Zhang的試驗(yàn)數(shù)據(jù)平緩，其原因可能是本文中的模型不考慮卸載后細(xì)觀裂紋的閉合效應(yīng)。相應(yīng)于應(yīng)力應(yīng)變曲線(xiàn)上A和B點(diǎn)的損傷和流速云，見(jiàn)圖5、6。在峰值應(yīng)力點(diǎn)A，損傷在基質(zhì)和顆粒接觸的附近區(qū)域開(kāi)始發(fā)展。在殘余應(yīng)力點(diǎn)B，損傷區(qū)域逐漸發(fā)展到充滿(mǎn)基質(zhì)。由于顆粒不透水，所以滲流只發(fā)生在基質(zhì)中，顆粒中沒(méi)有流速場(chǎng)。

a)應(yīng)力-應(yīng)變曲線(xiàn)

b)滲透系數(shù)-應(yīng)變曲線(xiàn)

a)A點(diǎn)

b)B點(diǎn)

a)A點(diǎn)

b)B點(diǎn)

4 結(jié)論

類(lèi)巖石材料的應(yīng)力-損傷-滲流耦合模型需要研究變形過(guò)程中的損傷演化和滲透系數(shù)隨損傷變化這兩個(gè)規(guī)律。本文提出一種改進(jìn)的彈性損傷模型；同時(shí)將孔隙壓密階段和裂紋擴(kuò)展階段的滲透系數(shù)變化過(guò)程統(tǒng)一起來(lái)，提出改進(jìn)的全應(yīng)力應(yīng)變過(guò)程中滲透系數(shù)-損傷演化模型。最后通過(guò)在水泥砂漿材料的細(xì)觀模型上實(shí)施單軸壓縮數(shù)值實(shí)驗(yàn)，得出以下結(jié)論。

a) 數(shù)值實(shí)驗(yàn)的應(yīng)力應(yīng)變曲線(xiàn)和室內(nèi)實(shí)驗(yàn)結(jié)果吻合，表明改進(jìn)的彈性損傷模型用來(lái)模擬類(lèi)巖石脆性材料的有效性。

b) 數(shù)值結(jié)果可以明顯識(shí)別滲透系數(shù)在應(yīng)力應(yīng)變?nèi)^(guò)程中的不同階段，在應(yīng)力達(dá)到峰值前由于孔隙不斷壓密，滲透系數(shù)不斷減??；在應(yīng)力峰值附近，應(yīng)變?yōu)?.004時(shí)，滲透系數(shù)突然增大，是因?yàn)榇藭r(shí)開(kāi)始出現(xiàn)損傷區(qū)域；在應(yīng)力峰值附近，應(yīng)變?yōu)?.007時(shí)，滲透系數(shù)達(dá)到最大值。

c) 由于缺少滲透系數(shù)的室內(nèi)實(shí)驗(yàn)數(shù)據(jù)，本文的滲透系數(shù)-損傷演化模型的參數(shù)還需進(jìn)一步通過(guò)實(shí)驗(yàn)確定，但滲透系數(shù)曲線(xiàn)趨勢(shì)基本能反映出Zhang[24]的試驗(yàn)數(shù)據(jù)規(guī)律。