考慮滲流應(yīng)力耦合的XFEM水壓致裂模擬

李雙 夏曉舟 裴磊

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

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

中圖分類號(hào)：TV139.14

文獻(xiàn)標(biāo)志碼：A

doi：10.3969/j .issn. 1000- 1379.2019.05.023

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

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

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

1.1 單元位移及壓力模式

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

1.1.1 XFEM單元位移模式

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

1.1.2 XFEM單元壓力模式

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

1.2 裂尖加強(qiáng)函數(shù)的改進(jìn)

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

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

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

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

忽略流體相的相對(duì)加速度，流固混合相的線性動(dòng)量平衡方程可表示為

1.3.3控制方程的XFEM空間離散

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

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

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

2.1 程序可靠性驗(yàn)證

為驗(yàn)證編寫的多孔介質(zhì)流固全耦合擴(kuò)展有限元程序的可靠性，將程序的計(jì)算結(jié)果與已有文獻(xiàn)計(jì)算結(jié)果進(jìn)行對(duì)比。計(jì)算模型（見圖3）為10 mxlo m的正方形斷裂域，中心含一長(zhǎng)度為2m的裂紋，為平面應(yīng)變模型。模型的位移邊界條件是左右兩邊約束U=0，底部約束U_v=0，頂部不約束;壓力場(chǎng)邊界條件是左右兩邊不排水，頂部為排水邊界，壓力設(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。時(shí)間步長(zhǎng)為0.01 s，計(jì)算總時(shí)間T= 10 s。

T=8 s時(shí)含30°斜裂紋的多孔介質(zhì)豎向位移云圖如圖4所示。可以看出，裂尖加強(qiáng)模式改進(jìn)前后的擴(kuò)展有限元程序所得云圖分布規(guī)律與文獻(xiàn)[1]結(jié)果吻合，證明了程序的正確性。兩種裂尖加強(qiáng)模式下的壓力云圖如圖5所示。由位移和壓力云圖可以看出，兩種裂尖加強(qiáng)模式都能很好地模擬含30°斜裂紋的多孔介質(zhì)流固全耦合問題。

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

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

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

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

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

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

3.1 流體黏度對(duì)孔隙水壓力的影響

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

參考文獻(xiàn)：

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

[2] 柳占立，莊茁，孟慶國(guó)，等，頁(yè)巖氣高效開采的力學(xué)問題與挑戰(zhàn)[J].力學(xué)學(xué)報(bào)，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]李錄賢，王鐵軍，擴(kuò)展有限元法（XFEM）及其應(yīng)用[J].力學(xué)進(jìn)展，2005（1）：5-20.

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

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

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