考慮滲流應力耦合的XFEM水壓致裂模擬

李雙 夏曉舟 裴磊

摘要：推導了含裂隙飽和多孔介質的水壓致裂全耦合問題的擴展有限元增量迭代求解格式。通過對含中心裂紋無限大平板的裂尖位移場解析解的分析，按I型和Ⅱ型變形特征進行提取，構建了一種新型的內聚裂尖加強函數，使傳統擴展有限元中受裂尖加強影響的節點由每個節點8個裂尖加強自由度縮減到2個（平面問題），大大減小了計算規模。以含中心斜裂紋的多孔介質平板為例，驗證了內聚裂尖加強模式的正確性，并對兩種裂尖加強模式的雅克比矩陣條件數進行了計算。結果表明：內聚裂尖加強模式下擴展有限元的計算精度與傳統擴展有限元的計算精度同階，且網格較密時，本文處理方法的計算精度高于傳統模式的。最后將改進的內聚裂尖加強模式應用到混凝土重力壩的水壓致裂模擬中，取得了良好的計算效果。

關鍵詞：擴展有限元;飽和多孔介質;滲流一應力耦合;水壓致裂;裂尖加強函數

中圖分類號：TV139.14

文獻標志碼：A

doi：10.3969/j .issn. 1000- 1379.2019.05.023

多孔介質[1-2]是由不同大小和形狀的固體顆粒骨架構成的，并至少有一種流體介質填充在介質的孔隙中。全部或部分飽和多孔介質的力學行為一般由固體骨架與孔隙流體的相互作用決定，因此如果考慮在恒定荷載作用下的土壤固結，通常會遇到流體流動和固體骨架變形的耦合問題。水壓致裂是油藏最常用的增產措施[2-3]。在傳統有限元法中，裂紋路徑局限于單元間邊界，通常遇到網格依賴性問題。擴展有限元法（ Extended Finite Element Method，簡稱XFEM） [4-5]是一種強大而精確模擬不連續介質的方法，該方法使整個裂紋的幾何形狀獨立于網格，并且完全避免了裂紋擴展時網格的重新劃分。對于裂紋擴展、孔洞、復合材料等含不連續問題的介質，可以用擴展有限元法有效模擬。XFEM通過引入具有特殊性質的加強函數及與之相應的加強自由度，模擬含間斷問題。

本文將多孔介質耦合控制方程與XFEM結合，完成對含裂隙飽和多孔介質的流固耦合數值模擬。同時對裂尖位移場加強模式進行改進，依據線彈性斷裂力學裂尖位移場的解析解形式，按I型和Ⅱ型裂紋張開變形特征進行提取，構建了一種新型的內聚裂尖加強函數，模擬裂尖的奇異性[6]，同時將原來裂尖加強模式中每個裂尖影響節點的8個裂尖加強自由度減到2個。通過算例分析證明，內聚裂尖加強函數不僅減小了矩陣規模，而且能更好地模擬含裂隙多孔介質的流固耦合問題，并運用內聚裂尖加強模型分析了混凝土大壩的流固耦合問題。

1 含裂隙多孔介質擴展有限元建模

1.1 單元位移及壓力模式

XFEM基于單元分解的基本思想，在傳統有限元法的近似解中引入合適的加強函數項來反映研究對象的某些特征，如裂紋或材料界面等的不連續性，同時使研究對象獨立于網格劃分。

1.1.1 XFEM單元位移模式

為了描述裂隙多孔介質中裂隙處的位移跳躍，假定位移場在裂紋邊緣上是強不連續的，并用Heaviside函數和裂尖漸進加強函數，分別描述裂紋面的不連續性及裂尖的奇異性，得到單元內任意一點的位移：

1.1.2 XFEM單元壓力模式

為了表示流體流量在垂直于裂縫方向上的跳躍，假設流體壓力場是弱不連續的，即壓力場連續，而其垂直于裂縫的梯度在裂縫邊緣是不連續的，并用修正的水平集函數描述壓力梯度的不連續性，單元中壓力表達式為[1]

1.2 裂尖加強函數的改進

由式（2）中的基函數組成的裂尖漸進函數F1（x），每個受裂尖影響的節點都有8個裂尖附加自由度[7]，增大了整體剛度矩陣的組裝及方程求解的計算量。為了減少裂尖附加自由度個數，本文以斷裂力學裂尖位移場的解析解為基函數，按I型和Ⅱ型裂紋張開的變形特征進行提取，在局部坐標系下重新構造能反映裂尖奇異性的內聚裂尖位移加強模式。整體坐標系、裂尖局部坐標及極坐標系示意見圖1。

1.3 含裂隙多孔介質的擴展有限元離散

1.3.1 多孔介質的控制方程

變形多孔介質的控制方程包括流固混合相及各相線性動量平衡方程和連續性方程，運用Biot理論得到[9]。

忽略流體相的相對加速度，流固混合相的線性動量平衡方程可表示為

1.3.3控制方程的XFEM空間離散

用XFEM的位移場及壓力場模式對弱形式的控制方程（23）進行空間離散，根據Bubnov Galerkin方法，可以得到積分方程即式（23）的離散形式：

2 含裂隙多孔介質擴展有限元程序計算

筆者依據含裂隙多孔介質的流固全耦合方程，編寫了裂尖加強模式改進前后的擴展有限元程序，計算了兩種裂尖加強模式下雅克比矩陣的條件數（2范數）。

2.1 程序可靠性驗證

為驗證編寫的多孔介質流固全耦合擴展有限元程序的可靠性，將程序的計算結果與已有文獻計算結果進行對比。計算模型（見圖3）為10 mxlo m的正方形斷裂域，中心含一長度為2m的裂紋，為平面應變模型。模型的位移邊界條件是左右兩邊約束U=0，底部約束U_v=0，頂部不約束;壓力場邊界條件是左右兩邊不排水，頂部為排水邊界，壓力設置為0，底部為流量人口邊界，q= lxl0^-4m/s。彈性模量E=9 GPa，泊松比為0.4，滲透率k=lxl0^-5m/s，孔隙度為0.3，固相壓縮模量為lx10²⁷Pa，流體壓縮模量為2.lx10⁹Pa，固相密度為2 000 kg/m3，流體相密度為1 000 kg/m³，流體黏度μU=lxl0^-3Pa.s。時間步長為0.01 s，計算總時間T= 10 s。

T=8 s時含30°斜裂紋的多孔介質豎向位移云圖如圖4所示?？梢钥闯?，裂尖加強模式改進前后的擴展有限元程序所得云圖分布規律與文獻[1]結果吻合，證明了程序的正確性。兩種裂尖加強模式下的壓力云圖如圖5所示。由位移和壓力云圖可以看出，兩種裂尖加強模式都能很好地模擬含30°斜裂紋的多孔介質流固全耦合問題。

2.2 雅克比矩陣條件數對比

為了比較兩種裂尖加強模式的優越性，分別以水平中心裂紋、30°及45°中心斜裂紋為例，對比分析兩種加強模式下豎向位移和水平位移云圖，見圖6-圖8。可以看出：對于水平裂紋和30°斜裂紋，兩種裂尖加強模式都能很好地模擬多孔介質的耦合問題。但對于45°斜裂紋，原來的加強模式所得結果出現上下裂紋面靠近的趨勢，而改進的加強模式所得結果并沒有出現，所以改進后模擬效果更好。

還可通過程序在4組不同的網格密度下分別查看裂尖加強模式改進前后雅克比矩陣的條件數。水平中心裂紋、30°及45°中心斜裂紋的雅可比矩陣條件數分別見表1-表3。

從雅克比矩陣的條件數來看，網格密度為11 x11時，改進前的雅可比矩陣條件數大于改進后的條件數，網格密度為2lx21時結果相反。網格密度小于3lx31時，兩種裂尖加強模式的條件數并無規律，但當網格密度增大時，裂尖加強模式改進后的條件數比加強模式改進前的條件數小，并且網格越大，兩者的差距越明顯。

3 混凝土重力壩水壓致裂模擬

將內聚裂尖加強模式運用到混凝土重力壩水壓致裂的滲流一應力全耦合分析中，并針對不同流體黏度和滲透系數對水壓致裂的裂縫影響進行對比，重力壩模型及參數見圖9、圖10及表4。將混凝土視為多孔介質材料，庫水按牛頓流體處理。

3.1 流體黏度對孔隙水壓力的影響

不同流體黏度下混凝土大壩裂隙的孔隙水壓力分布如圖11所示。

可以看出，流體黏度從μ= lx 10^-1Pa.s到μ=1×10^-6Pa.s，裂隙注水口附近的孔隙水壓力在逐漸變小，即流體黏度對注水口的孔隙水壓力有影響，水壓力隨黏度的增大而變大。

3.2 滲透系數對孔隙水壓力及張開度的影響

混凝土大壩的滲透能力由裂隙的滲透系數表示，滲透系數是影響混凝土重力壩水壓致裂的主要因素之一，可以通過改變介質的滲透系數來控制裂隙的萌生及裂縫擴展過程。不同滲透系數下大壩的孔隙水壓力及裂隙張開位移隨滲透系數的變化見圖12和圖13。

可以看出，在注水開始一段時間后，滲透系數為30x 10^-9m/s時的孔隙水壓力最小，為lx 10^-9m/s時孔隙水壓力最大，可見在這一段時間里孔隙水壓力與滲透系數成反比關系[11]。同樣，在圖13中，裂隙張開位移隨滲透系數的增大而減少，與文獻[11]的結果相符。滲透系數越大，孔隙水壓力越小，對裂隙的張開作用越小，則相對裂隙張開位移越小。

4 結論

多孔介質的滲流應力耦合問題一直是研究的熱點，本文對含裂隙飽和多孔介質流固全耦合問題進行了數值建模，并結合擴展有限元法建立了含裂隙多孔介質的有限元模型。為了減少裂尖附加自由度個數，對裂尖位移場加強模式進行改進，以線彈性斷裂力學裂尖位移場的解析解為基函數，構建內聚裂尖加強模式，將裂尖加強自由度由8個減到2個。通過含中心裂隙的無限平板算例驗證了模型的可靠性，查看了改進前后雅克比矩陣的條件數，比較了兩種裂尖加強模式的優越性，同時將改進的裂尖加強模式運用到了含裂隙混凝土大壩中，對大壩的參數進行了簡單分析，證明了改進的裂尖加強在運用中的可行性。

（1）由位移云圖和壓力云圖可知，在裂紋傾角較小時，改進后的加強模式和改進前的加強模式都能很好地模擬含裂隙多孔介質的流固耦合問題。但當裂紋傾角較大時，改進后的內聚裂尖加強模式模擬效果更好。

（2）由條件數可知，當網格密度較稀疏時，兩種裂尖加強模式的條件數并無規律，但當網格密度增大時，2個附加自由度的加強模式條件數比8個附加自由度的加強模式條件數小，并且網格越大，兩者的差距越明顯。

（3）運用內聚裂尖加強對混凝土大壩的裂隙進行模擬，得出孔隙水壓力受流體黏度及滲透系數的影響，同時隨著水壓力一起變化的還有裂隙張開位移，孔隙水壓力越小，裂隙張開位移越小。

綜上，內聚裂尖加強模式下的擴展有限元計算精度與傳統擴展有限元計算精度同階，且網格較密時內聚裂尖加強模式的計算精度高于傳統模式的，這種方法不僅縮小了矩陣規模，而且能更好地模擬含裂隙多孔介質的流固全耦合問題，將其應用到混凝土重力壩的水壓致裂模擬能得到良好的計算效果。

參考文獻：

[1] KHOEI A R.Extended Finite Element Method：‘rheory andApplications[M].[S.l.]：[s.n.]， 2014： 409-460.

[2] 柳占立，莊茁，孟慶國，等，頁巖氣高效開采的力學問題與挑戰[J].力學學報，2017，49（3）：507-516.

[3] 楊克己，韓理安，樁基工程[M].北京：人民交通出版社，1992：2- 10.

[4]

MOHAMMADI S.Extended Finite Element Method： for FractureAnalysis of Structures[J].Xfem Fracture Analysis of Composites，2007，45（22）：5675-5687.

[5]李錄賢，王鐵軍，擴展有限元法（XFEM）及其應用[J].力學進展，2005（1）：5-20.

[6]夏曉舟，趙睿祺，章青，基于擴展有限元方法的含裂紋結構在動力荷載作用下的數值仿真[ C]//南方計算力學聯絡委員會，江蘇省力學學會，第十一屆南方計算力學學術會議（ SCCM-II）摘要集，南方計算力學聯絡委員會，江蘇省力學學會，2017：1.

[7] 江守燕，杜成斌，一種XFEM斷裂分析的裂尖單元新型改進函數[J].力學學報，2013，45（1）：134-138.

[8] 王榮華，章青，夏曉舟，含裂隙飽和多孔介質流一固耦合的擴展有限元分析[J].巖土力學，2017，38（5）：1489-1496.

[9] 王忠福，劉漢東，周敏，等，黃河大堤丁壩流固耦合分析研究[J].人民黃河，2008，30（8）：16-17.

[10] RETHORE J，DE BORST R，ABELLAN M A.ATwo-Scale Approach for Fluid Flow in Fractured Porous Media[J] .Intemational Joumal for Numerical Methods in Engi-neering， 2007， 71（7）：780- 800.

[11]

MOHAMMADNEJAD T，KHOEI A R.An Extended FiniteElement Method for Hydraulic Fracture Propagation in De-formable Porous Media with the Cohesive Crack Model[J].Finite Elements in Analysis and Design， 2013，73：77-95.