深埋隧洞滲流與應力耦合分析

張繼勛 盛 亮 任旭華 于牧萍

（1.水文水資源與水利工程科學國家重點實驗室，南京 210098；2.河海大學 水利水電學院，南京 210098；3.河海大學 土木與交通學院，南京 210098）

深部巖體一般處于高地應力、高滲透水壓、高地溫等特殊地質環境中，因此巖體的力學行為和擾動響應明顯有別于淺層巖體［1]，巖體的基本性質和擾動響應不是單個因素造成的，而是多種因素相互耦合共同作用的結果，在開展工程數值分析的時候也只能就地質條件選擇主要的因素進行分析.應力場和滲流場是富水區深埋隧洞影響圍巖穩定性的兩個主要因素，20世紀70年代，Louis首次研究了自然狀態下非破壞巖體中的滲流－應力耦合作用，他提出的滲透系數和正應力呈負指數關系也為其他學者所證實［2]；Witherspoon［3]提出了可變性介質中滲流－應力的模型；王媛等［4]在裂隙巖體滲流－應力耦合的實驗、模擬和計算方法的研究也開展了大量卓有成效的工作，當前的研究成果雖然很豐富，但在工程實際應用上由于一些計算規模和巖體實際條件所限還受到限制.本文以高水位富水區隧道圍巖為研究對象，以有效應力原理為基礎，應用滲透系數和孔隙度隨體積應變改變而改變的動態演化模型，實現了滲流與應力完全耦合，分析了滲流場與應力場的相互影響程度，給工程實際的滲流場與應力場耦合分析以及圍巖穩定性研究提供了依據.

1 應力平衡方程和滲流連續方程

1.1 平衡方程

巖土介質應力平衡方程可采用虛功原理來表示，即在某一時刻巖土體的虛功作用在該巖土體上作用力（體力和面力）產生的虛功相等，即

式中，t為面力，f為體力；δε、δu分別為虛位移和虛應變.

由于滲流連續性方程含有時間項，為了將應力和滲流進行耦合，需對虛功方程進行時間的求導，經過推導具體表示式為

1.2 連續方程

考慮某一體積巖石，根據質量守恒原理，在dt時間流入該體積內的水量應等于其內部儲水量的增加，流體的滲流用達西定律來描述，經過推導可得出滲流的連續方程為

式中，k0為初始滲透系數張量與水密度的乘積；kr為比滲透系數，可為飽和度、應力、應變或損傷變量等的函數；g為重力加速度矢量；n為孔隙度；kw為水的體積模量.

1.3 有限元方程的離散

定義形函數

根據式（4）和式（2），經過簡化可以得到固相有限元列式

在滲流場分析中，有兩類邊界條件：一個是流量邊界條件，另一個是孔壓邊界條件.流量邊界條件可以表示為

式中，n為流量邊界的單位法向；qw是單位時間內流過邊界的水流量.

孔壓邊界條件可表示為

式中，Pwb是已知邊界處的孔隙壓力值.

使用Gakerkin方法，即

式中，a、b為任意函數為控制方程為通過邊界的連續方程.

經過簡化應力－滲流耦合方程可表示為

2 裂隙巖體應力滲流耦合本構模型

2.1 基于多孔介質的裂隙巖體有效水力原理

應用李培超［5]等在論文中的觀點，假定工程巖體為多孔介質，推導出基于多孔介質的有效應力原理.該原理包含多孔介質的結構參數孔隙度φ，代替了其他有效應力公式中用的較多的經驗參數（比如常用的biot常數），其有效應力原理為

該原理較好地體現了多孔介質中的流固耦合效應，對于飽和巖體而言，p就是孔隙水壓力.

2.2 裂隙巖體滲透系數動態演化模型

對于流固耦合工程巖體，孔隙度和滲透系數率K等參數將隨巖體的應力狀態不同而發生動態變化，因此有必要建立流固耦合作用下的動態模型.

根據李培超等［5]建立的飽和多孔介質流固耦合滲流的數學模型，可得到多孔介質孔隙度與體積應變、溫度、應力等有如下的關系

式中，φ0為初始孔隙度；εv為體積應變.其表達式為為多孔介質骨架固體顆粒的體積彈性壓縮模量；βs為熱膨脹系數.

若不考慮滲流工程中溫度和骨架顆粒的體積變化，無擴容現象時，其孔隙度可由式（12）得到

在擴容條件，可以得到在壓縮條件下孔隙度φ的動態演化模型為

由Kozeny－Carman方程導出滲透系數體積應變和關系式為

同樣，如果不考慮溫度為材料骨架顆粒的體積變化，則可以得等溫滲流過程中滲透系數的動態演化模型為

對于擴容后的情況，采用同樣的分析方法，可以得到壓縮條件下滲透系數的動態演化模型為

公式（14）～（17）反映了巖體在應力作用下使裂隙閉合時將孔隙度和滲透系數減小的特點，以及當微裂隙擴展使得巖體產生擴容現象時，其孔隙度和滲透系數將相應的增大的特性.

3 算 例

大型商業軟件ABAQUS有豐富的接口來支持二次開發，對于裂隙巖體流固耦合的巖體本構關系，可以通過ABAQUS的用戶子程序USDFLD來實現滲流與應力的完全耦合，在程序中可以定義巖體損傷、孔隙度或滲透系數變化等實現本構的改變［6].在耦合計算中，ABAQUS能將滲流場與應力場直接耦合，無需進行滲流場與應力場的反復迭代，只要按時間過程連續求解就可得到全部結果，即通過將節點位移和孔隙水壓力作為節點自由度進行空間離散，將應力平衡方程和滲流連續方程寫成矩陣形式，并對滲流連續方程引入時間積分，得到耦合控制方程，然后在每個時間步求解方程，并同時滿足位移邊界條件和滲流邊界條件.

3.1 計算模型與邊界條件

本算例的研究是一個埋深為1 180m的深埋隧洞的開挖過程，隧洞的洞徑為12m，計算范圍上下左右都取了洞徑的15倍，所以模型的大小為360m×360 m的二維模型，如圖1所示.本算例分析的目的是為了說明滲流與應力的耦合分析在深埋地下工程中的必要性，所以在分析的過程中用到了裂隙巖體滲透系數動態演化模型，即式（13）和（14）.巖體的彈模E為20 000MPa，泊松比μ為0.25，巖體的密度ρ為2 800 kg／m3，內聚力c為1.5MPa，內摩擦角φ為42°，抗拉強度σt為2MPa，孔隙度φ0為0.003，滲透系數K為4.63×10－9m／s，因為要考慮滲透數的動態演化，所以計算的過程中彈模和滲透系數是隨著體積應變的改變而改變的，巖體材料采用D－P屈服準則.

圖1 計算模型網格圖

邊界條件為：在模式的兩側施加x方向的位移約束，在底面上施加y方向上的位移約束，因為模型的上表面距地面有1 000m，所以在上面加了27.44MPa的均勻面壓力，上表面初始水頭為300m，下表面的初始水頭為660m，左右兩側施加沿重力方向上梯度變化的水頭壓力.

3.2 結果分析

分別分析了考慮耦合的情況和不考慮耦合的情況，如圖2～4所示.可以看出，一方面開挖應力場對滲透系數的影響距離洞周越遠影響越小，不同位置影響的規律略有差別，滲流場對應力場分布規律的影不大；另一方面隧洞開挖引起二次應力重分布對滲流場有明顯的影響，但是滲流場的變化對圍巖的應力分布的影響不大，因此在隧洞開挖過程中，應力場對滲流場的影響處于主要地位，而滲流場對應力場的影響處于次要地位.為了簡化分析，在滿足工程要求的前提下可以按照如下做法：滲流場分析宜考慮應力場變化的影響；應力場分析可以忽略滲流場影響.

4 結 論

介紹了滲流與應力耦合在深埋地下工程中應用的基本理論，對于裂隙巖體，介紹了裂隙巖體滲透系數的動態演化模型，通過算例說明了深埋地下工程中進行耦合分析的必要性.在滿足工程要求的前提下根據分析目的可取其重要的方面.

［1]王金昌，陳頁開.ABAQUS在土木工程中的應用［M].杭州：浙江大學出版社，2006：318.

［2]White W B.Geomorphology and Hydrology of Karst Terrains［M].Oxford：Oxford University Press，1988.

［3]Witherspoon.Validity of Cubic Law for Fluid Flow in a Deformable Rock Fracture［J].Water Resources Research，1980，16（6）：1016－1024.

［4]王 媛，劉 杰.基于敏感性分析的裂隙巖體滲流與應力靜態全耦合參數反演［J].巖土力學，2009，30（2）：311－317.

［5]李培超，李賢桂，盧德唐.飽和土體一維固結理論的修正－飽和多孔介質流固耦合滲流模型之應用［J].中國科學技術大學學報，2010，40（12）：1273－1276.

［6]陳衛忠，伍國軍，戴永浩，等.錦屏二級水電站深埋引水隧洞穩定性研究［J].巖土工程學報，2008（8）：1184－1190.