基于損傷的滲流-應力耦合模型在工程地下水遷移研究中的應用

張 萊,陸桂華

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

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

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

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

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

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

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

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

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

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

式中:ν為泊松比。

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

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

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

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

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

式中:D1,D2,D3為相應于3個主應力方向的損傷分量。因D1≠D2≠D3,所以共有5個獨立變量。

2 基于損傷的滲流-應力耦合模型

2.1 巖體滲流有限元計算基本方程

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

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

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

2.2 有限元計算格式

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

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

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

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

式(14)對已知水頭邊界節(jié)點形成常數(shù)項。通過式(14)計算以后,將常數(shù)項移到等號右端,得N個未知節(jié)點的線性代數(shù)方程組。對時間項取隱式有限差分,寫為矩陣形式有

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

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

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

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

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

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

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

3 工程實例

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

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

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

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

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

圖5 三維滲流計算模型效果示意圖

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

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

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

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

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

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

e.總體滲流量偏大,右岸比左岸滲流量小,以壩肩繞滲為主。壩基滲流量以防滲帷幕底部繞滲為主,滲流量屬中等。相比而言左岸滲流量較大,主要是因為山梁部位有多個強滲透斷層帶貫通上下游,因此該部位斷層帶的防滲處理應該加強。

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

4 結(jié)論與討論

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

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

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

:

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

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

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

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

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

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

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

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

[9]陳祖安,伍向陽.巖石滲透率隨靜水壓力變化的關(guān)系研究[J].巖石力學與工程學報,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)模型及工程應用[J].巖土工程學報,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]安民,俞茂宏,吳熹,等.廣義雙剪屈服準則在巖石力學中的應用[J].巖土力學,1991,12(1):17-26.

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