李雙 夏曉舟 裴磊

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

關(guān)鍵詞：擴展有限元;飽和多孔介質(zhì);滲流一應(yīng)力耦合;水壓致裂;裂尖加強函數(shù)

中圖分類號：TV139.14

文獻標志碼：A

doi：10.3969/j .issn. 1000- 1379.2019.05.023

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

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

1 含裂隙多孔介質(zhì)擴展有限元建模

1.1 單元位移及壓力模式

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

1.1.1 XFEM單元位移模式

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

1.1.2 XFEM單元壓力模式

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

1.2 裂尖加強函數(shù)的改進

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

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

1.3.1 多孔介質(zhì)的控制方程

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

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

1.3.3控制方程的XFEM空間離散

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

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

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

2.1 程序可靠性驗證

為驗證編寫的多孔介質(zhì)流固全耦合擴展有限元程序的可靠性，將程序的計算結(jié)果與已有文獻計算結(jié)果進行對比。計算模型（見圖3）為10 mxlo m的正方形斷裂域，中心含一長度為2m的裂紋，為平面應(yīng)變模型。模型的位移邊界條件是左右兩邊約束U=0，底部約束U_v=0，頂部不約束;壓力場邊界條件是左右兩邊不排水，頂部為排水邊界，壓力設(shè)置為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°斜裂紋的多孔介質(zhì)豎向位移云圖如圖4所示。可以看出，裂尖加強模式改進前后的擴展有限元程序所得云圖分布規(guī)律與文獻[1]結(jié)果吻合，證明了程序的正確性。兩種裂尖加強模式下的壓力云圖如圖5所示。由位移和壓力云圖可以看出，兩種裂尖加強模式都能很好地模擬含30°斜裂紋的多孔介質(zhì)流固全耦合問題。

2.2 雅克比矩陣條件數(shù)對比

為了比較兩種裂尖加強模式的優(yōu)越性，分別以水平中心裂紋、30°及45°中心斜裂紋為例，對比分析兩種加強模式下豎向位移和水平位移云圖，見圖6-圖8?？梢钥闯觯簩τ谒搅鸭y和30°斜裂紋，兩種裂尖加強模式都能很好地模擬多孔介質(zhì)的耦合問題。但對于45°斜裂紋，原來的加強模式所得結(jié)果出現(xiàn)上下裂紋面靠近的趨勢，而改進的加強模式所得結(jié)果并沒有出現(xiàn)，所以改進后模擬效果更好。

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

參考文獻：

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

[2] 柳占立，莊茁，孟慶國，等，頁巖氣高效開采的力學問題與挑戰(zhàn)[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）及其應(yīng)用[J].力學進展，2005（1）：5-20.

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

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

[8] 王榮華，章青，夏曉舟，含裂隙飽和多孔介質(zhì)流一固耦合的擴展有限元分析[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.