張 萊,陸桂華

(1.河海大學(xué)水文水資源學(xué)院,江蘇南京 210098;2.中國(guó)煤炭地質(zhì)總局水文地質(zhì)工程地質(zhì)環(huán)境地質(zhì)勘查院,河北邯鄲 056004;3.江蘇省水利廳,江蘇南京 210029)

巖石損傷是工程巖體在水的滲流與應(yīng)力相互作用下的變形和破壞,是許多工程共同關(guān)心的問(wèn)題。在巖土邊坡、礦井突水、油汽開發(fā)、水電站大壩運(yùn)行、核廢料地下處理的污染泄漏等工程中,地下巖石與地下溶洞受高應(yīng)力、高水壓的作用,尤其是裂隙,在高壓、高水頭和復(fù)雜的化學(xué)侵蝕作用下,裂隙的發(fā)育發(fā)展影響到整個(gè)工程的安全。當(dāng)前國(guó)內(nèi)外有許多學(xué)者致力于滲流-應(yīng)力耦合問(wèn)題的研究[1-3],并且提出了多種理論分析和數(shù)值模擬模型[4-7]。多數(shù)學(xué)者的研究偏向于流場(chǎng)(流體力學(xué))作用于固體(巖體)產(chǎn)生反應(yīng)(包括巖體的變形、破壞、穩(wěn)定性等)[8-10],只有少數(shù)學(xué)者研究固體骨架(巖體)的純力學(xué)作用。

為探求一種清晰明了的耦合關(guān)系模型,本文基于巖石受載過(guò)程強(qiáng)度弱化源于內(nèi)部缺陷生長(zhǎng)的概念,推導(dǎo)了基于應(yīng)變(變形)的損傷演化方程,描述了巖石強(qiáng)度變化的非線性特征,并采用應(yīng)變-損傷-滲透性關(guān)系方程描述耦合關(guān)系,開發(fā)了工程計(jì)算效果較好的滲流-應(yīng)力耦合的非線性計(jì)算程序。

1 巖石損傷本構(gòu)關(guān)系及損傷演化方程

1.1 三維損傷變量及演化方程

引入Lemaitre損傷變量εc描述三維損傷并構(gòu)建演化方程:

式中:D為三維損傷變量;ε為軸向應(yīng)變;εc為常量;n為材料脆性參數(shù),與材料脆性有關(guān)。

圖1 損傷材料應(yīng)力應(yīng)變關(guān)系

如圖1[11]所示,εc是當(dāng)n→∞時(shí)的彈性應(yīng)變上界值,包括破壞前應(yīng)變和破壞后的殘余應(yīng)變。即損傷的演化取決于應(yīng)變狀態(tài)的變化,應(yīng)變狀態(tài)的變化又與損傷狀態(tài)有關(guān)。一維狀態(tài)的損傷演化方程可以在主應(yīng)力空間擴(kuò)展到三維。對(duì)脆性巖石材料而言,在裂紋尖端的局部區(qū)域內(nèi)可能不會(huì)出現(xiàn)塑性屈服,而是脆性破壞后的微小碎末區(qū)域包含在未破壞的完整巖石之中。因此,在描述損傷演化的時(shí)候就不能認(rèn)為巖石破壞后仍然具有承載能力,殘余應(yīng)變也就無(wú)從談起。εc為包含殘余應(yīng)變的容許值,應(yīng)當(dāng)去除殘余應(yīng)變,其值等于峰值應(yīng)力對(duì)應(yīng)的應(yīng)變。采用最大拉應(yīng)變準(zhǔn)則來(lái)判定材料是否破壞是合適的。臨界狀態(tài)為

試驗(yàn)表明混凝土材料在單軸壓縮時(shí)的損傷是張應(yīng)力引起的,損傷方向與載荷方向正交。這表明混凝土的損傷是各向異性或正交各向異性的[12]。堅(jiān)硬巖石材料與混凝土材料性質(zhì)相似,可認(rèn)為巖石材料的損傷是正交各向異性的。在主應(yīng)力空間,應(yīng)力主軸、應(yīng)變主軸和損傷主軸互相重合。材料初始狀態(tài)是各向同性的線彈性體,損傷后表現(xiàn)出正交各向異性性質(zhì)。在此假定下,損傷張量D和應(yīng)力張量σ均為二階張量,設(shè)3個(gè)主損傷分量為D1,D2,D3,則根據(jù)熱力學(xué)第二定律、一維損傷的定義式及損傷變量的物理意義,并令

式中:ν為泊松比。

在單軸壓縮狀態(tài)(σ1＞0,σ2=σ3=0),有 ε1＞0,ε2=ε3=ε′,ε′＜0,則式(3)變?yōu)?/p>

從巖石試件的單軸壓縮試驗(yàn)觀察到試件破壞的方向大體與壓應(yīng)力方向一致,與式(4)中D1＜D2=D3相符。另外,按式(4)σ3方向也有損傷。這是因?yàn)楫?dāng)微裂紋出現(xiàn)后,微裂紋的方向并不完全是破壞時(shí)的宏觀裂紋方向,初始微裂紋的方向存在一定的隨機(jī)性,這也正是D1的含義。在數(shù)量上D1遠(yuǎn)遠(yuǎn)小于D2和D3。

1.2 彈脆性材料損傷本構(gòu)方程

無(wú)損狀態(tài)下柔度矩陣為

則由廣義虎克定律及式(6)推出彈脆性損傷材料的本_構(gòu)關(guān)系:

式中:D1,D2,D3為相應(yīng)于3個(gè)主應(yīng)力方向的損傷分量。因D1≠D2≠D3,所以共有5個(gè)獨(dú)立變量。

2 基于損傷的滲流-應(yīng)力耦合模型

2.1 巖體滲流有限元計(jì)算基本方程

符合Darcy定律的非均質(zhì)非各向同性可壓縮巖土體的三維空間非穩(wěn)定滲流不考慮巖土體的壓縮性(Ss=0),得到飽和穩(wěn)定滲流的控制方程:

初始條件為計(jì)算開始水頭場(chǎng)分布狀態(tài),邊界條件包括已知水頭邊界和已知流量邊界:

式中:kl為法向上的滲透系數(shù);q為過(guò)潛流面的已知單位面積流量,q=0為不透水邊界,q≠0為潛流邊界;μ為給水度;Γ1為已知水頭邊界;Γ2為已知流量邊界;Γ3為自由面邊界。

2.2 有限元計(jì)算格式

根據(jù)變分原理,上述定解問(wèn)題等價(jià)于下列泛函極小值問(wèn)題:

將計(jì)算域離散化,以Ie(h)表示單元體 Ωe的泛函,即

分別對(duì)任意單元e進(jìn)行Ie1,Ie2,Ie3求導(dǎo)數(shù)和極小值變換,將所有單元的泛函求微分后疊加,并利用I(h)極小值的條件,有

式中:N′i為以i為公共節(jié)點(diǎn)的單元數(shù)。

式(14)對(duì)已知水頭邊界節(jié)點(diǎn)形成常數(shù)項(xiàng)。通過(guò)式(14)計(jì)算以后,將常數(shù)項(xiàng)移到等號(hào)右端,得N個(gè)未知節(jié)點(diǎn)的線性代數(shù)方程組。對(duì)時(shí)間項(xiàng)取隱式有限差分,寫為矩陣形式有

這就是有限單元法最后求解的線性代數(shù)方程。

2.3 滲流-應(yīng)力耦合模型有限元網(wǎng)格劃分與計(jì)算

有限元計(jì)算采用等參數(shù)六面體單元作為基本單元,特殊部位使用六面體單元的退化單元。標(biāo)準(zhǔn)六面體等參數(shù)單元及其局部坐標(biāo)如圖2所示,各種退化單元形狀如圖3所示。

圖2 六面體八節(jié)點(diǎn)等參數(shù)單元

圖3 六面體等參數(shù)單元的各種退化形式

對(duì)六面體八節(jié)點(diǎn)等參數(shù)單元的基函數(shù)進(jìn)行局部坐標(biāo)與整體坐標(biāo)之間的變換[14],利用高斯積分法進(jìn)行數(shù)值積分得

式中:Al,Am,Ak為加權(quán)系數(shù);ξl,ηm,ζk為積分點(diǎn)局部坐標(biāo)。

3 工程實(shí)例

某水電站壩址區(qū)河谷呈開闊不對(duì)稱的V形,右岸山體雄厚,左岸為一較單薄的山脊梁,河谷深切,兩岸卸荷強(qiáng)烈,傾倒變形嚴(yán)重。由于壩址區(qū)分布的較軟弱巖層傾角近直立,河谷不斷被侵蝕下切,河谷兩岸巖體原有的平衡被破壞,兩岸巖層均向河谷產(chǎn)生傾倒變形。為了評(píng)價(jià)壩區(qū)和山梁巖層及斷層帶的滲透性質(zhì)、檢驗(yàn)防滲措施的有效性,提出改良意見,制定了如下三維有限元滲流及滲透穩(wěn)定計(jì)算方案。

根據(jù)前面推導(dǎo)的公式和有限元網(wǎng)格劃分方法,采用南京水利科學(xué)研究院于1974年開發(fā)并不斷完善的UNSS3程序進(jìn)行三維滲流計(jì)算,應(yīng)用于工程實(shí)踐。計(jì)算采用容重替代法處理靜水壓力,等效節(jié)點(diǎn)力方法處理動(dòng)水壓力。三維有限元計(jì)算單元以八節(jié)點(diǎn)六面體(圖2)為主,并包括了5種形式的退化單元(圖3),以便適應(yīng)復(fù)雜的地下結(jié)構(gòu)與巖(斷)層產(chǎn)狀。滲流計(jì)算模型共劃分27390個(gè)單元、30 080個(gè)節(jié)點(diǎn)。剖分后的三維網(wǎng)格見圖4。

圖4 三維滲流計(jì)算模型網(wǎng)格剖分

邊界條件:三維滲流計(jì)算邊界條件包括邊界地下水位、地表水位、地表出滲及特殊工況的內(nèi)部邊界。設(shè)計(jì)擬定的特征水位如下:正常蓄水位為1408.00m,相應(yīng)的下游水位為1309.80m;設(shè)計(jì)洪水位為1408.26m,相應(yīng)的下游水位為1324.74m;校核洪水位為1411.92m,相應(yīng)的下游水位為1327.15m。

三維滲流計(jì)算模型范圍:上游邊界至壩腳垂直距離約為300m;下游河床至壩腳垂直距離約為350m。左右兩岸由于山高坡陡,分水嶺距壩址較遠(yuǎn),且較遠(yuǎn)處缺乏地質(zhì)勘察資料,因此右岸取到右壩肩以外約500m的位置、左岸取到溢洪道以外約700m的位置作為計(jì)算域邊界。左右邊界以地質(zhì)勘察地下水位作為已知水頭邊界條件進(jìn)行三維滲流計(jì)算。考慮到河床部位灌漿帷幕的深度,計(jì)算模型底部邊界取到1160.00m高程。三維計(jì)算模型的效果圖見圖5。

圖5 三維滲流計(jì)算模型效果示意圖

通過(guò)三維天然滲流場(chǎng)參數(shù)反演和極限工況下三維滲流數(shù)值計(jì)算(表1),得出以下結(jié)論:

a.反演計(jì)算結(jié)果表明巖層及斷層帶滲透性具有各向異性的性質(zhì)。

b.斷層帶具有良好的導(dǎo)水性,承擔(dān)的坡降小于允許比降,但流速較大,應(yīng)注意防止在巖體裂隙附近出現(xiàn)集中滲漏造成的沖刷破壞。左壩肩附近斷層(F109,F121,F122,F123)具備導(dǎo)水特征,因此應(yīng)在上游采取有效的擋水措施,以保證滲透穩(wěn)定;溢洪道部位斷層(F124,F125,F126,F127)具有排水作用,不影響溢洪道的滲透穩(wěn)定。

表1 反演計(jì)算得到的滲透參數(shù) 10-4cm/s

c.左岸山梁部位排水孔幕的布置能夠有效地降低山梁內(nèi)浸潤(rùn)面高度,山梁部位出逸比降不大,能夠滿足滲透穩(wěn)定要求。

d.壩體浸潤(rùn)面低,由于面板及防滲帷幕承擔(dān)了絕大部分水頭差,壩體堆石部分的浸潤(rùn)面基本在建基面以上5.0m的位置,出逸點(diǎn)位置也不高(高程1310.62m),不存在滲透失穩(wěn)問(wèn)題。

e.總體滲流量偏大,右岸比左岸滲流量小,以壩肩繞滲為主。壩基滲流量以防滲帷幕底部繞滲為主,滲流量屬中等。相比而言左岸滲流量較大,主要是因?yàn)樯搅翰课挥卸鄠€(gè)強(qiáng)滲透斷層帶貫通上下游,因此該部位斷層帶的防滲處理應(yīng)該加強(qiáng)。

從三維滲流有限元計(jì)算成果來(lái)看,工程滲控措施基本上是合適的,能夠有效降低壩頭及左岸山體的浸潤(rùn)面高度。但是,右岸帷幕的防滲效果不顯著,建議在下階段進(jìn)行右岸帷幕的延展長(zhǎng)度及深度的敏感性分析?？傮w上流量比較大,主要是左岸山梁的強(qiáng)導(dǎo)水?dāng)鄬訋Я髁看?局部大流量會(huì)帶來(lái)沖刷破壞問(wèn)題,因此下階段設(shè)計(jì)工作應(yīng)該加強(qiáng)對(duì)斷層帶的防滲處理措施。建議取消左岸山梁引水管下游的灌漿帷幕,或在該位置建設(shè)排水孔幕,這樣效果可能會(huì)更好。

4 結(jié)論與討論

通過(guò)研究巖石破壞過(guò)程中的損傷演化規(guī)律,以及與之相伴而生的滲透性變化規(guī)律,引入隨機(jī)概念描述巖石細(xì)觀強(qiáng)度分布,得到了初始強(qiáng)度非均勻的巖石試樣,采用應(yīng)變主軸、損傷主軸和附加滲透主軸一致的原則,建立了基于應(yīng)變 -損傷- 附加滲透性、反映巖石強(qiáng)度變化非線性特征的損傷演化模型,使得耦合計(jì)算過(guò)程物理意義清晰明了,更準(zhǔn)確地表達(dá)了耦合過(guò)程的非線性關(guān)系。通過(guò)工程實(shí)例表明筆者建立的基于損傷的滲流- 應(yīng)力耦合模型是正確可行的。

巖石的破壞是從局部開始的,但是當(dāng)巖石被作為均勻介質(zhì)時(shí),基于統(tǒng)計(jì)方法的損傷力學(xué)很難真實(shí)地描述巖石損傷的局部化現(xiàn)象。基于拉應(yīng)變的損傷變量既可以描述壓縮狀態(tài)的巖石試件,又適用于拉伸狀態(tài)的巖石試件。在材料內(nèi)部應(yīng)力集中的區(qū)域,應(yīng)變異常,損傷也相應(yīng)地發(fā)生異常,這就是巖石破壞局部化的原因。但是有限元計(jì)算微元體損傷后的力學(xué)形態(tài)還需要做進(jìn)一步研究,使之能夠準(zhǔn)確模擬損傷的局部化擴(kuò)展。

巖石損傷的最后結(jié)果是巖石斷裂破壞,基于連續(xù)介質(zhì)力學(xué)體系的損傷力學(xué)還不能準(zhǔn)確地描述巖石的斷裂過(guò)程以及斷裂狀態(tài)。將損傷和斷裂理論相結(jié)合來(lái)研究巖石的破壞過(guò)程將是值得研究的技術(shù)方法。在這一點(diǎn)上,損傷因子和斷裂指標(biāo)的結(jié)合以及由損傷到斷裂的物理機(jī)制將是斷裂損傷聯(lián)合研究的重點(diǎn),也是難點(diǎn),還有待于深入研究。

:

[1]黃秋楓,胡海浪.滲流-應(yīng)力-損傷耦合研究現(xiàn)狀[J].災(zāi)害與防治工程,2009,66(5):27-34

[2]楊天鴻,唐春安,徐濤,等.巖石破裂過(guò)程的滲流特性:理論、模型與應(yīng)用[M].北京:科學(xué)出版社,2004:32.

[3]張文捷,程榮蘭,詹美禮,等.巖體滲流的一種改進(jìn)數(shù)學(xué)模型[J].河海大學(xué)學(xué)報(bào):自然科學(xué)版,2010,38(1):52-57.

[4]王環(huán)玲,徐衛(wèi)亞,童富國(guó).泄洪霧雨區(qū)裂隙巖質(zhì)邊坡飽和-非飽和滲流場(chǎng)與應(yīng)力場(chǎng)耦合分析[J].巖土力學(xué),2008,29(8):2397-2403.

[5]宋曉晨,徐衛(wèi)亞.裂隙巖體滲流模擬的三維離散裂隙網(wǎng)格數(shù)值模型(Ⅰ):裂隙網(wǎng)格的隨機(jī)生成[J].巖石力學(xué)與工程學(xué)報(bào),2004,23(12):2015-2020.

[6]盧剛,周志芳.降雨入滲下互層狀裂隙巖體非飽和滲流分析[J].巖土工程學(xué)報(bào),2008,30(9):1399-1403.

[7]黃勇,陳靜.巖體裂隙介質(zhì)各向異性耦合模型研究[J]河海大學(xué)學(xué)報(bào):自然科學(xué)版,2009,37(2):171-174.

[8]王俊光,梁冰.滲透動(dòng)水壓力作用下裂隙巖體滲流與應(yīng)力耦合分析[J].遼寧工程技術(shù)大學(xué)學(xué)報(bào),2009,28(4):178-180.

[9]陳祖安,伍向陽(yáng).巖石滲透率隨靜水壓力變化的關(guān)系研究[J].巖石力學(xué)與工程學(xué)報(bào),1995,14(2):155-159.

[10]ODA M T,TAKEMURA A,AOKI T.Damage growth and permeability change in triaxial compression tests of india granite[J].Mechanics of Materials,2002,34:313-331.

[11]周維垣,剡公瑞,楊若瓊.巖體彈脆性損傷本構(gòu)模型及工程應(yīng)用[J].巖土工程學(xué)報(bào),1998,20(5):54-57.

[12]PLESHA M E.Constitutive modeling of rock joints with dilation[C]//26th US Symposium on Rock Mechanics.Rapid City,SD(USA):American Rock Mechanics Association,1985:387-392.

[13]安民,俞茂宏,吳熹,等.廣義雙剪屈服準(zhǔn)則在巖石力學(xué)中的應(yīng)用[J].巖土力學(xué),1991,12(1):17-26.

[14]魏澤光.三維有限元法中等參數(shù)六面體單元的單元分析[J].高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào),1979,1(1):3-9.