基于連續(xù)介質(zhì)的裂隙巖體流固耦合數(shù)值分析

唐佳，彭振斌，何忠明

?

基于連續(xù)介質(zhì)的裂隙巖體流固耦合數(shù)值分析

唐佳1, 2，彭振斌1，何忠明3

(1. 中南大學地球科學與信息物理學院，湖南長沙，410083；2. 中國電建集團中南勘測設計研究院有限公司地質(zhì)工程處，湖南長沙，410007；3. 長沙理工大學交通運輸工程學院，湖南長沙，410076)

基于庫區(qū)水位漲落或暴雨入滲，經(jīng)過演變易造成滑坡，研究滲透壓作用下裂隙損傷效應對滲透張量的影響，通過理論推導，提出裂隙巖體滲流場與損傷場耦合的分析模型；以湖北潘口水電站進水口邊坡為研究對象，建立基于連續(xù)介質(zhì)的裂隙巖體流固耦合滲流數(shù)值計算分析方法。研究結果表明：邊坡在開挖過程中變形以卸荷回彈變形為主，回彈變形范圍隨開挖逐級增大，各級邊坡開挖對鄰級邊坡變形的影響最顯著，對較遠坡段的影響逐步減弱；邊坡開挖面基本處于小拉應力狀態(tài)，邊坡淺表層拉應力區(qū)分布范圍及量值均很??；邊坡全部開挖完成后，邊坡的塑性區(qū)主要沿近壩的開挖馬道外緣分布，分布范圍小，在表層巖體內(nèi)；裂隙巖體流固耦合模型的三維數(shù)值計算能夠較真實地模擬邊坡開挖過程位移場、應力場、塑性區(qū)、拉應力區(qū)的分布規(guī)律。

滑坡；滲透壓力；連續(xù)介質(zhì)；流固耦合；數(shù)值計算

水電站庫區(qū)水位漲落或暴雨入滲易造成滑坡。通過對水庫滑坡進行實測分析，中村浩之等[1]認為造成水庫滑坡的主要因素有水庫浸水、水庫水位急劇下降與降雨3個方面；王思敬等[2]根據(jù)水庫滑坡產(chǎn)生的原因?qū)⑵鋭澐譃樘烊换潞驮谒畮r作用影響下的滑坡2類；王士天等[3]則認為，當庫水位已達到敏感水位或者庫水位消落且處在快速消落期時很容易發(fā)生滑坡。唐曉松等[4]發(fā)現(xiàn)水庫滑坡不同于一般山體滑坡, 其穩(wěn)定性受水位波動的影響十分明顯。嚴福章等[5]通過對清江隔河巖水庫坪滑坡進行研究，提出在水庫蓄水后，該滑坡發(fā)生了較長的緩慢變形，其根本原因是水庫蓄水產(chǎn)生的材料和水力學效應共同作用以及滑帶抗剪強度在滑體變形后降低，產(chǎn)生了應變軟化作用。張有天等[6]通過對60座水庫滑坡進行統(tǒng)計后指出，85%的滑坡在建設期、蓄水期或者在水庫存工程完工之后2 a內(nèi)發(fā)生。裂隙巖質(zhì)岸坡位移速率隨著庫水位上漲不斷加大，在庫水位降低時不斷減小甚至會停止。李佳偉等[7]通過巖石力學試驗，揭示了砂板巖巖體的高空隙水壓力效應，建立了主要力學參數(shù)預測模型?？梢?，為了探討地質(zhì)災害的真正原因，必須從巖石水力學的角度對巖體的滲透規(guī)律和損傷力學響應進行研究，并采取流固耦合分析法對裂隙岸坡失穩(wěn)機制進行分析。湖北潘口水電站庫岸大部分為白云母石英片巖，巖性較軟弱，巖石風化較強烈，加之區(qū)域構造活動頻繁，節(jié)理、構造發(fā)育，易產(chǎn)生順層滑坡和沿結構面與層面組合的滑坡，為此，本文作者以潘口水電站庫岸為研究對象，同時考慮靜水壓力和動水壓力的影響，提出基于等效連續(xù)介質(zhì)的裂隙巖體流固耦合滲流分析數(shù)值模型，同時對考慮滲流與應力相互作用下進水口邊坡進行數(shù)值分析，研究相應部位邊坡的位移場、應力場、塑性區(qū)、拉應力區(qū)分布規(guī)律，系統(tǒng)評價邊坡的穩(wěn)定性。

1 裂隙巖體連續(xù)介質(zhì)滲流模型

1.1 等效滲透系數(shù)的計算

同理，用表示裂縫法線測線的長度，表示裂縫的數(shù)目，則巖體線密度d為

(2)

基于巖塊的滲透系數(shù)m、寬度m及裂隙面的滲透系數(shù)為f、寬度c，由達西定律得到巖體的滲流量：

(4)

式中：為巖體寬度；f為單條裂隙寬度。

則巖塊寬度為

可得巖體等效滲透系數(shù)K：

(6)

垂直面其等效滲透系數(shù)為

則

(8)

采用線連續(xù)性系數(shù)將裂隙的實際滲透系數(shù)t修正為

K和K又可以表示為：

(10)

K和K反映了巖體滲透特性的各向異性。假定裂紋的角度為，則等效二階張量的滲透率為：

(12)

(14)

由此獲得了裂隙巖體二維等效滲透系數(shù)。

1.2 滲透壓作用下巖體的滲透系數(shù)

考慮巖體受到滲透水壓力和應力共同耦合作 用[8]，且其滲透系數(shù)一般隨應力環(huán)境的改變而變化，滲透系數(shù)張量可表示為

式中：為裂隙半徑；為裂隙隙寬；為裂隙法向矢量。相應滲透量為

(16)

(17)

(＞1) (18)

1.3 裂隙巖體流固耦合機理

據(jù)文獻[9?10]，假定巖體水頭值為，則靜水壓為

(20)

當單裂隙水平分布時，裂隙中的動水壓為

(22)

1.4 巖體滲流場和應力場耦合數(shù)學模型的建立

式中：巖體的等效滲透系數(shù)為應力場的函數(shù)，即；s為貯水率；為巖體的滲流區(qū)間；和分別為初始時刻及1上的水頭分布；1和2分別為已知水頭邊界及流量邊界；為2邊界上法向量；為2上的流量分布。巖體的平衡方程為

式中：f和f為分別為和的巖體質(zhì)量力。巖體的幾何方程為

(26)

巖體的本構方程為

巖體的應力邊界條件為

1.5 耦合方程的建立

同時在裂隙巖體滲流中考慮應力對滲流的影響，引進耦合方程：

2 流固耦合數(shù)值計算

在模擬巖體的流固耦合機理時，視流體滿足Biot方程，依據(jù)Darcy定律在巖體孔隙介質(zhì)中的流動[11?12]，得以下幾個主要方程。

2.1 平衡方程

巖體滲流平衡狀態(tài)時，

(31)

液體質(zhì)量平衡方程為

動量平衡方程為

(33)

其中：g為重力加速度在方向的分量；，為體積密度；s和w分別為固體和液體密度；(1?)s為巖石干密度即d。

2.2 運動方程

用達西定律描述流體的運動[13]：

式中：為介質(zhì)的滲透系數(shù)；f為流體密度：g為重力加速度在方向的分量；x為方向的坐標分量。

2.3 本構方程

孔隙介質(zhì)本構方程的增量形式為

無形的勁氣開始在天問大師和蕭飛羽身前流轉(zhuǎn)，漸漸地，勁氣將血醒的塵土帶著飛了起來。天問大師雙手如挽重物緩緩提到胸前，緊接著向外一翻，一股強大的勁氣應手而出。蕭飛羽沒有作勢，右掌隨意一翻?！稗Z！”一聲悶響，塵土被激起、飛射、旋轉(zhuǎn)、墜落。

2.4 液體響應方程

完全飽和時，=1，響應方程為

同時有

(37)

當＜1時，響應方程為

相對滲透性與飽和度的關系為

(39)

對于小變形響應方程為

(41)

2.5 相容方程

應變與速度的關系為

式中：u為介質(zhì)的速度。

2.6 邊界條件

耦合計算可采用以下透水邊界條件[14]：

式中：n為邊界外法線方向流速分量；為滲漏系數(shù)；e為滲流出口處的孔隙水壓力。

3 邊坡應力變形有限元分析

3.1 邊坡巖體工程特性

潘口水電站岸塔式進水口布置在左岸壩前。邊坡地形較陡，大部分地段基巖裸露，巖性為白云母石英片巖夾云母片巖，局部含炭質(zhì)條帶，富含石英脈；緩坡部位分布有殘、坡積物。巖層片理走向為320°~330°/ NE∠40°~55°，發(fā)育的節(jié)理主要有280°~290°/ SW∠40°~60°，345°~355°/SW∠40°~60°和60°~70°/ SE∠40°~70°共3組。邊坡平面布置及其相互關系如圖1所示(比例尺為1:500，圖中數(shù)字1~6代表剖面編號)，邊坡巖體的物理力學參數(shù)取值見表1。

3.2 有限元計算模型及邊界條件

采用的自主二次開發(fā)的有限元軟件EVP3D進行有限元建模及計算。模型的建造遵循“點→線→面→體”自下而上的建模技術。模型中采用六面體單元和部分四面體退化單元模擬。模型共剖分單元18 194個，結點14 932 個，其中開挖單元2 581 個，如圖2所示。邊坡開挖后的模型如圖3所示?？紤]模型所處的地形地貌條件及邊坡荷載方向，模型邊界條件，和向分別取法向支座約束。

圖1 進水口邊坡平面圖

表1 邊坡巖體力學參數(shù)

圖2 開挖前的進水口邊坡三維有限元模型

3.3 計算荷載及施工過程模擬

邊坡計算荷載包括巖體自重應力場、開挖荷載等。施工期開挖荷載是最主要的計算荷載，其計算式為[15]

式中：{}為開挖荷載向量；[]為單元幾何矩陣；[]為形函數(shù)矩陣；{}為初始應力；{}為體積力向量；E為開挖單元集合。本次計算在初始地應力場(主要考慮自重應力場)模擬的基礎上，按307，327，347，362，382，402和422 m共7個開挖量模擬進水口邊坡的開挖過程。

圖3 開挖后的進水口邊坡三維有限元模型

Fig. 3 Three-dimensional finite element model of intake-slope after excavation

3.4 位移成果分析

進水口邊坡在開挖過程中的整體位移云圖如圖4所示。為了分析進水口邊坡開挖過程中的變形規(guī)律，在邊坡開挖坡面上選取13 個關鍵點，如圖5所示。部分關鍵點(1，3，5，7，9，11和13)隨邊坡開挖的位移過程曲線如圖6所示。

從圖4~6可看出進水口邊坡在開挖過程中具有如下變形規(guī)律：

1) 邊坡在開挖過程中，各開挖坡面上的變形以卸荷回彈變形為主，回彈變形的方向垂直向上，并略指向山體內(nèi)部，邊坡完全開挖后產(chǎn)生的最大變形出現(xiàn)在307 m高程的基坑中央約14.06 mm處。

2) 邊坡各坡段在開挖過程中產(chǎn)生的變形增量與開挖量有關，開挖量越大，邊坡的變形增量也越大，且高高程坡面相比低高程坡面變形小。

3) 從開挖產(chǎn)生變形的變化過程可知開挖引起的回彈變形范圍，隨著開挖的深入而逐級增大；各級邊坡變形受鄰級邊坡開挖的影響最嚴重，受較遠邊坡開挖的影響則逐層減弱；而開挖坡段受后斷邊坡開挖影響邊坡變形遞增，最終則趨向平穩(wěn)。

(a) 開挖步1; (b) 開挖步3; (c) 開挖步5; (d) 開挖步7

圖5 進水口邊坡關鍵點位置示意圖

關鍵點：1—K1; 2—K3; 3—K5; 4—K7; 5—K9; 6—K11; 7—K13。

3.5 應力計算成果分析

邊坡在開挖過程中3?3 剖面第一主應力云圖如圖7所示，第三主應力云圖如圖8所示，全部開挖完成后的拉應力區(qū)分布如圖9所示。計算分析結果表明：邊坡面在開挖過程中基本上處于小拉應力狀態(tài)，且處于表層巖體之中深度約5 m以內(nèi)，這種應力狀態(tài)整體上對邊坡穩(wěn)定性影響很?。幻恳徊介_挖產(chǎn)生的拉應力基本上都處于表層巖體中，隨開挖量增大逐漸減小但變化不大，后2步(即327 m和307 m)開挖在開挖馬道坡面上產(chǎn)生很小范圍的壓應力區(qū)；每一步開挖產(chǎn)生的壓應力也基本上都處于表層巖體中，全部開挖完成后小主應力分布在307~402 m 之間的坡段上，最大值約為2.5 MPa，出現(xiàn)在307 m高程基坑內(nèi)側(cè)靠近坡腳處。

圖7 3?3剖面第7步開挖第一主應力云圖

圖8 3?3剖面第7步開挖第三主應力云圖

圖9 進水口邊坡完全開挖后拉應力區(qū)分布

3.6 塑性區(qū)分布

邊坡開挖完成后的整體塑性區(qū)如圖10所示。

分析圖10可知：邊坡的塑性區(qū)主要沿近壩馬道外緣分布，在高程402，382和362 m馬道外緣分布較集中，高程327 m馬道也有小范圍分布?？傮w上，開挖引起的塑性區(qū)范圍較小，深度在表層巖體8 m之內(nèi)，對開挖邊坡整體穩(wěn)定性影響不大。針對邊坡淺表的塑性變形，應及時進行淺表層的錨桿支護，以確保邊坡局部穩(wěn)定。

圖10 進水口邊坡完全開挖后塑性區(qū)分布

4 結論

1) 考慮滲透壓作用下裂隙損傷效應對滲透張量的影響，通過理論推導，提出了裂隙巖體滲流場與損傷場耦合的分析模型。

2) 建立了基于等效連續(xù)介質(zhì)的裂隙巖體流固耦合滲流數(shù)值計算分析方法。對潘口水電站進水口邊坡進行模擬驗證，結果表明邊坡在開挖過程中變形以卸荷回彈變形為主，范圍隨開挖逐級增大；各級邊坡開挖對鄰級邊坡變形的影響最顯著，對較遠坡段的影響逐步減弱；邊坡開挖面基本處于小拉應力狀態(tài)，淺表層拉應力區(qū)分布范圍及量值均很小。邊坡全部開挖完成后，邊坡的塑性區(qū)主要沿近壩的開挖馬道外緣分布，分布在表層巖體范圍內(nèi)。

3) 采用裂隙巖體流固耦合模型的三維數(shù)值計算能較真實地模擬邊坡開挖過程位移場、應力場、塑性區(qū)、拉應力區(qū)的分布規(guī)律，相對于傳統(tǒng)的二維極限平衡分析，其對工程設計與施工具有重要的指導意義。

參考文獻：

[1] 中村浩之, 王恭先. 論水庫滑坡[J]. 水土保持通報, 1990, 10(1): 53?64. Nakamura Hiroyuki, WANG Gongxian. Creep on reservoir landslide[J]. Bulletin of Soil and Water Conservation, 1990, 10(1): 53?64.

[2] 王思敬, 馬鳳山, 杜永康. 水庫地區(qū)的水巖作用及其地質(zhì)環(huán)境影響[J]. 工程地質(zhì)學報, 1996, 4(3): 1?9. WANG Sijing, MA Fengshan, DU Yongkang. On the rock-water interaction in reservoir areas and its geoenvironmental effect[J]. Journal of Engineering Geology, 1996, 4(3): 1?9.

[3] 王士天, 劉漢超, 張悼元, 等. 大型水域水巖相互作用及其環(huán)境效應研究[J]. 地質(zhì)災害與環(huán)境保護, 1997, 8(1): 69?89. WANG Shitian, LIU Hanchao, ZHANG Daoyuan, et al. Study on water rock interaction and its environmental effect in large water area[J]. Journal of Geological Hazards and Enviroment Preservation, 1997, 8(1): 69?89.

[4] 唐曉松, 鄭穎人, 唐輝明, 等. 水庫滑坡變形特征和預測預報的數(shù)值研究[J]. 巖土工程學報, 2013, 35(5): 940?947. TANG Xiaosong, ZHEN Yingren, TANG Huiming, et al. Numerical researches on deformation characteristics and prediction of reservoir landslides[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(5): 940?947.

[5] 嚴福章, 王思敬, 徐瑞春. 清江隔河巖水庫蓄水后茅坪滑坡的變形機理及其發(fā)展趨勢研究[J]. 工程地質(zhì)學報, 2003, 11c(1): 15?24. YAN Fuzhang, WANG Sijing, XU Ruichun. Creeping deformation of the Maoping landslidean ditssbabm in the Geheyan reservoir on the Qingjiang river[J]. Journal of Engineering Geology, 2003, 11c(1): 15?24.

[6] 張有天. 巖石水力學與工程[M]. 北京: 中國水利水電出版社, 2005: 11?47. ZHANG Youtian. Rock hydraulics and engineering[M]. Beijing: China Water and Power Press, 2005: 11?47.

[7] 李佳偉, 徐進, 王璐, 等. 砂板巖巖體力學特性的水巖耦合試驗研究[J]. 巖土工程學報, 2013, 35(3): 599?604. LI Jiawei, XU Jin, WANG Lu, et al. Water-rock coupling tests on mechanical properties of sandy slate rock mass[J]. Chinese Journal of Geotechnical Engineering, 2013, 35(3): 599?604.

[8] 趙延林, 曹平, 趙陽升, 等. 雙重介質(zhì)滲流場?應力場?溫度場耦合模型及三維數(shù)值研究[J]. 巖石力學與工程學報, 2007, 26(增2): 4024?4031. ZHAO Yanlin, CAO Ping, ZHAO Shengyang, et al. Dual media model for thermo-hydro-mechanical coupling and 3D numerical simulation[J]. Chinese Journal of Rock Mechanics and Engineering, 2007, 26(Suppl 2): 4024?4031.

[9] 易順民, 朱珍德. 裂隙巖體損傷力學導論[M]. 北京: 科學出版社, 2005: 9?37. YI Shunmin, ZHU Zhende. Introduction to damage mechanics of fractured rock mass[M]. Beijing: Science Press, 2005: 9?37.

[10] 朱維申, 張強勇. 節(jié)理巖體脆彈性斷裂損傷模型及其工程應用[J]. 巖石力學與工程學報, 1999, 18(3): 245?249. ZHU Weishen, ZHANG Qiangyong. Fracture damage model of jointed rock mass and its engineering application[J]. Chinese Journal of Rock Mechanics and Engineering, 1999, 18(3): 245?249.

[11] 楊延毅. 節(jié)理裂隙巖體損傷斷裂力學模型及其在巖體工程中的應用[D]. 北京: 清華大學土木工程學院, 1990: 13?41. YANG Yanyi. Damage fracture mechanics model of Joint fracture and application[D]. Beijing: Tsinghua University. School of Civil Engineering, 1990: 13?41.

[12] KEMENY J M, COOK N G W. Crack models for the failure of rocks in compression,constitutive laws for engineering materials[M]. New York: Elsevier, 1987: 879?887.

[13] 趙延林. 裂隙巖體滲流?損傷?斷裂耦合理論及應用研究[D]. 長沙: 中南大學資源與安全工程學院, 2009: 16?48. ZHAO Yanlin. Coupling theory of seepage?damage?fracture in fractured rock masses and application[D]. Changsha: Central South University. School of Resources and Safety Engineering, 2009: 16?48.

[14] 謝東海, 馮濤, 趙延林, 等. 裂隙煤巖體的流固耦合精細模型[J]. 中南大學學報(自然科學版), 2014, 45(6): 2014?2020. XIE Donghai, FENG Tao, ZHAO Yanlin, et al. Fluid-solid coupling detailed model of rock-coal mass[J]. Journal of Central South University (Science and Technology), 2013, 44(5): 2014?2020.

[15] 陳濤, 吳旭敏, 陳曉江. 錦屏二級上游調(diào)壓室洞室群開挖支護動態(tài)設計與研究[J]. 巖石力學與工程學報, 2014, 33(Sl): 2593?2599. CHEN Tao, WU Xumin, CHEN Xiaojiang. Dynamic design and research on excavationg and support design of upper surge chamber caverns of Jinping II hydropower station[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(Sl): 2593?2599.

(編輯 陳燦華)

Numerical analysis of fluid-solid coupling in fractured rock mass based on continuous medium

TANG Jia1, 2, PENG Zhenbin1, HE Zhongming3

(1. School of Geosciences and Info?Physics, Central South University, Changsha 410083, China;2. Geological Engineering Department, Powerchina Zhongnan Engineering Co. Ltd., Changsha 410007, China;3. School of Traffic and Transportation Engineering, Changsha University of Science and Technology, Changsha 410076, China)

Considering that reservoir water level fluctuation or rainstorm infiltration will cause landslide, the effect of crack damage effect on the permeability tensor was studied. Through theoretical deduction, the analysis model of coupling of seepage field and damage field in fractured rock mass was put forward. Then taking the Hubei Pankou hydropower station intake-slope as object, a numerical analysis method of fluid solid coupling seepage flow in fractured rock mass was established based on the equivalent continuous medium. The results show that deformation of slope is mainly based on unloading springback during excavation; deformation range increases with the increase of excavation; the influence of slope excavation at all levels on the deformation of adjacent level slope is the most significant; distance slope section is gradually weakened. After the completion of the excavation, slope of the plastic zone is mainly layed along the near dam outer edge of the distribution, and the distribution range is small and in the surface of rock mass. Three-dimensional numerical calculation of fluid solid coupling model of fractured rock mass can simulate the distribution law of the displacement field, stress field, plastic zone and tensile stress zone of the slope excavation process.

landslide; osmotic pressure; continuous medium; fluid?solid coupling; numerical calculation

10.11817/j.issn.1672-7207.2016.11.025

TU457

A

1672?7207(2016)11?3800?08

2016?01?10；

2016?03?22

國家自然科學基金資助項目(51508042) (Project(51508042) supported by the National Natural Science Foundation of China)

唐佳，博士研究生，工程師，從事地質(zhì)工程、巖土工程研究；E-mail: csutangjia@126.com