李卓徽，周宗青，高成路，張道生，白松松

（1. 浙江大學(xué)建筑工程學(xué)院，浙江 杭州 310058；2. 山東大學(xué)齊魯交通學(xué)院，山東 濟(jì)南 250002；3. 陸軍勤務(wù)學(xué)院軍事設(shè)施系，重慶 401311）

近幾十年來介質(zhì)中的流固耦合數(shù)值研究已成為力學(xué)、土木和環(huán)境工程等領(lǐng)域的主要研究課題，該研究涉及到固體變形和裂隙演化擴(kuò)展以及裂隙區(qū)流體流動。為了求解這類問題，出現(xiàn)了各種復(fù)雜的數(shù)值模型，可簡單分為3 類：連續(xù)、離散和混合方法。典型的連續(xù)方法主要是基于有限元法，但此類方法難以解決裂紋尖端應(yīng)力場奇異問題；離散方法通常采用離散元法，在宏細(xì)觀參數(shù)標(biāo)定方面存在復(fù)雜性和局限性；為了結(jié)合連續(xù)方法和離散方法的優(yōu)點(diǎn)，一些混合方法被提出［1］，但這些方法均在處理多裂紋、裂紋分叉問題時仍面臨挑戰(zhàn)，特別是在三維條件下。

近場動力學(xué)（PD）理論最早是由美國Sandia 國家實(shí)驗(yàn)室的Silling［2］在2000年提出。其基本思想是以非局部作用的積分模型代替?zhèn)鹘y(tǒng)理論的微分模型，可以解決有限元方法中尖端奇異性的問題，同時，兼具分子動力學(xué)方法和無網(wǎng)格方法的優(yōu)點(diǎn)，又突破了經(jīng)典分子動力學(xué)方法在計算尺度上的局限［2］。這使得近場動力學(xué)理論可以描述從連續(xù)到不連續(xù)、從微觀到宏觀的一系列力學(xué)行為。Zhou等［3］提出了一種考慮壓縮荷載作用下巖體峰后階段的微觀彈塑性本構(gòu)模型。Gao 等［4-6］將近場動力學(xué)應(yīng)用于隧道開挖，分析了在不同強(qiáng)度折減系數(shù)下不同節(jié)理傾角的節(jié)理巖體的破壞模式，以及模擬了在隧道開挖過程中巖體漸進(jìn)破壞導(dǎo)致突水通道形成的演化過程。與此同時，近場動力學(xué)理論還能模擬流體滲流問題，在統(tǒng)一框架下實(shí)現(xiàn)流-固耦合模擬是近場動力學(xué)理論的又一大突出優(yōu)勢。因此，越來越多的學(xué)者將近場動力學(xué)理論用于巖石材料流固耦合問題的模擬研究中。壽云東［7］針對巖土工程中的溫度場-滲流場-應(yīng)力場耦合問題，引入最大拉應(yīng)力強(qiáng)度準(zhǔn)則、莫爾-庫侖強(qiáng)度準(zhǔn)則和雙剪強(qiáng)度準(zhǔn)則，建立了基于近場動力學(xué)理論的裂隙巖體溫度場-滲流場-應(yīng)力場耦合數(shù)值模型，實(shí)現(xiàn)了巷道圍巖失穩(wěn)破壞過程的模擬。王允騰［8］專注鍵型近場動力學(xué)理論研究，針對巖體熱-水-力-化（THMC）耦合作用機(jī)理，提出了THMC 耦合作用下的鍵型近場動力學(xué)數(shù)值模型，從而實(shí)現(xiàn)了THMC 耦合作用下巖體破裂過程模擬。然而，由于近場動力學(xué)涉及到非局部的相互作用，與經(jīng)典模型相比，近場動力學(xué)理論在離散化時需要更多的計算資源，這一點(diǎn)在模擬流固耦合問題時尤為明顯。為了在合理的計算時間內(nèi)模擬各種多尺度物理過程，國內(nèi)外學(xué)者提出了近場動力學(xué)與經(jīng)典模型耦合的不同方案。Macek和Silling［9］利用嵌入節(jié)點(diǎn)和單元將PD 網(wǎng)格與有限元（FEM）網(wǎng)格耦合。Kilic 和Madenci［10］在PD 域和有限體積（FV）域的界面處引入了一個重疊區(qū)域，在 該 區(qū) 域 內(nèi) 求 解2 個 方 程。 Agwai 等［11］和Oterkus［12］采用了FEM 進(jìn)行整體建模，又采用PD進(jìn)行預(yù)測材料失效的子建模。Lubineau 等［13］采用基于能量等效的變形函數(shù)進(jìn)行耦合。這些函數(shù)只影響本構(gòu)參數(shù)，因此允許模型作為純粹非局部的、純 粹 局 部 的 或 混 合 的。Liu 和Hong［14］將FEM 和PD 域與界面元耦合，設(shè)計了2 種不同的耦合方案。Ni 等［15］提出了一種基于PD-FEM 耦合的模擬飽和多孔介質(zhì)中水力裂縫擴(kuò)展的混合建模方法?？梢钥闯隼媒鼒鰟恿W(xué)和連續(xù)介質(zhì)理論耦合的方法在不降低求解精度的情況下有效地減少計算量已經(jīng)成為了目前研究的一大熱點(diǎn)問題。

受Ni 等［15］啟發(fā)，提出一種有限體積法與近場動力學(xué)耦合的方法來模擬飽和孔隙裂隙介質(zhì)中水力裂縫擴(kuò)展問題，以此提高計算效率，并開展室內(nèi)試驗(yàn)尺度下多工況巖體水力壓裂模擬，驗(yàn)證該方法在流體驅(qū)動下模擬飽和裂隙孔隙介質(zhì)中裂紋擴(kuò)展的能力。

1 基本原理

1.1 普通態(tài)型近場動力學(xué)

近場動力學(xué)根據(jù)物質(zhì)點(diǎn)間相互作用的不同分為：鍵型（BB-PD）、普通態(tài)型（OSB-PD）和非普通態(tài)型（NOSB-PD）3 種不同類型［16-19］。鍵型近場動力學(xué)模型中對點(diǎn)力只考慮物質(zhì)點(diǎn)間相對位置，力的大小和方向均相同，而在描述物體力學(xué)性質(zhì)的過程中，鍵型近場動力學(xué)最棘手的問題在于材料的泊松比固定，即二維情況下所描述的材料泊松比恒等于1/3［19］；三維情況下所描述的材料泊松比恒等于1/4［19］，這極大地限制了近場動力學(xué)在工程領(lǐng)域的應(yīng)用范圍。

Silling等［17］在2007年提出了態(tài)型近場動力學(xué)理論。在該理論中，Silling 指出計算區(qū)域內(nèi)離散物質(zhì)點(diǎn)的彈性應(yīng)變能可分解為各向同性膨脹應(yīng)變能和形狀改變應(yīng)變能［17-18，20］，提出了“狀態(tài)”的概念，類似于傳統(tǒng)連續(xù)介質(zhì)力學(xué)張量概念，并通過引入變形矢量狀態(tài)和力密度矢量狀態(tài)使近場動力學(xué)能夠與傳統(tǒng)連續(xù)介質(zhì)力學(xué)理論有效結(jié)合，從而有效解決泊松比限制問題。

普通態(tài)型近場動力學(xué)的運(yùn)動方程為［17］

根據(jù)普通態(tài)型近場動力學(xué)理論模型［12］，可將力矢量狀態(tài)寫為

式中：θ為近場動力學(xué)理論中的體積膨脹；-ed為伸長狀態(tài)的偏張量；k、α為近場動力學(xué)材料參數(shù)，二者對應(yīng)于傳統(tǒng)彈性理論中的材料體積模量和剪切模量；m為加權(quán)體積，m5。

二維平面應(yīng)變問題下的普通態(tài)型近場動力學(xué)力密度標(biāo)量狀態(tài)表達(dá)式為［20］

其中

式中：K和G分別為體積模量和剪切模量；ν為泊松比。

1.2 有限體積法概述

有限體積法（Finite Volume Method，F(xiàn)VM）借鑒了有限元的部分優(yōu)點(diǎn)，以有限差分法（FDM）為基礎(chǔ)發(fā)展而來。同F(xiàn)EM 和FDM 一樣，將計算區(qū)域劃分成有限體積大小的離散網(wǎng)格，即進(jìn)行離散化處理，每個網(wǎng)格節(jié)點(diǎn)按照一定的方式形成表征該節(jié)點(diǎn)所對應(yīng)的控制體積，如圖2 所示。其關(guān)鍵思想在于針對每一個控制體積將待解的微分方程進(jìn)行積分形式處理，從而得出一組離散方程。其中的未知數(shù)是網(wǎng)格點(diǎn)上的因變量的數(shù)值。

圖2 有限體積法示意Fig.2 Schematic diagram of finite volume method

有限體積法的基本思路易于理解，具有具體的物理含義。離散方程實(shí)際上表示的是控制體積的通量平衡，也就是通量在有限大小的控制體積中的守恒原理，如同微分方程表示因變量在無限小的控制體積中的守恒原理一樣。有限體積法得出的離散方程使得因變量的積分守恒對任意控制體積都滿足，也就是說有限體積法的局部表述可以實(shí)現(xiàn)局部守恒，那么整個計算區(qū)域自然也得到滿足，進(jìn)而能夠以自然、直接的方法來解決對流占主導(dǎo)的流動問題的離散化，則特征變量φ在控制體積內(nèi)的守恒關(guān)系為

式中：φt為φ隨時間的變化量；φc為邊界對流進(jìn)入控制體積的量；φd為邊界擴(kuò)散進(jìn)入控制體積的量；φs為源項(xiàng)產(chǎn)生的量。

對比有限元法，可以發(fā)現(xiàn)有限元法的一個缺點(diǎn)是，對于連續(xù)試函數(shù)和基函數(shù)，沒有局部守恒，只能保證全局守恒。換句話說，只能保證域邊界上的凈通量是平衡的。另一個缺點(diǎn)是無法控制局部通量，這意味著穩(wěn)定對流占主導(dǎo)地位的流動的離散化并不簡單。此外，有限體積法推導(dǎo)其離散方程是通過控制體積的積分方程作為出發(fā)點(diǎn)，這與有限差分法直接從微分方程推導(dǎo)完全不同。所以，有限差分法僅當(dāng)網(wǎng)格極其細(xì)密時，離散方程才滿足積分守恒；而有限體積法即使在粗網(wǎng)格情況下，也能顯示出準(zhǔn)確的積分守恒。

綜上，有限體積法的特點(diǎn)可以總結(jié)為以下四點(diǎn)：①控制方程以積分形式表示，這有別于有限差分法。②積分方程體現(xiàn)了局部守恒性，即變量在任意控制體積內(nèi)守恒，這與有限元法不同。③離散化后，各個離散方程的每一項(xiàng)均有相應(yīng)的物理解釋，這是有限體積法最具優(yōu)勢的地方。④計算區(qū)域內(nèi)離散的各個網(wǎng)格節(jié)點(diǎn)之間控制體積不重疊，從而離散的局部變量守恒保證了整個計算區(qū)域內(nèi)場變量的守恒，實(shí)現(xiàn)整體守恒。

2 近場動力學(xué)與有限體積法耦合理論

2.1 固體模塊計算

為了模擬固體變形與孔隙壓力變化之間的相互作用，需要考慮流體壓力的影響和變形的影響，對近場動力學(xué)理論的運(yùn)動方程進(jìn)行修正。參考在單相流體飽和多孔介質(zhì)理論［21-22］中的有效應(yīng)力原理，其表達(dá)式為

式中：σeff為有效應(yīng)力張量；σtotal為總應(yīng)力張量；α為Biot系數(shù)；p為孔壓；I為單位張量。

在近場動力學(xué)理論中的標(biāo)量力密度狀態(tài)可看作是一個膨脹項(xiàng)和一個偏差項(xiàng)的加和，而且該膨脹項(xiàng)對應(yīng)于經(jīng)典連續(xù)介質(zhì)力學(xué)理論中的體積應(yīng)力，利用式（11）所示的有效應(yīng)力原理，將孔隙壓力項(xiàng)直接添加進(jìn)標(biāo)量力密度狀態(tài)的膨脹項(xiàng)中，引入標(biāo)量有效力密度狀態(tài)，則標(biāo)量力密度狀態(tài)與標(biāo)量有效力密度狀態(tài)之間的關(guān)系可寫成類似于式（11）所給出的總應(yīng)力張量與有效應(yīng)力張量之間的關(guān)系。其中，-ttotal在三維情況下對應(yīng)于式（3）中的-t，在二維平面應(yīng)變情況下對應(yīng)于式（4）中的-t，在二維平面應(yīng)力情況下對應(yīng)于（7）中的-t。

參考文獻(xiàn)［15，23-25］思路，將整個計算區(qū)域劃分為固體域（Rs）、過渡域（Rt）和裂隙域（Rf）三部分。同時，運(yùn)用近場動力學(xué)損傷場作為劃分依據(jù)，設(shè)置2個指標(biāo)用來標(biāo)識3個區(qū)域，如圖3所示。所得線性指標(biāo)函數(shù)為

圖3 計算區(qū)域劃分Fig.3 Division of calculation area

則利用線性指標(biāo)函數(shù)對固體域和裂隙域的滲透率進(jìn)行插值，可以得到過渡域的固體滲透率

另外，式（12）中的標(biāo)量有效力密度狀態(tài)取決于物質(zhì)點(diǎn)對所處的計算區(qū)域，可分為以下3種情況［26］：

（1）物質(zhì)點(diǎn)對處于固體域（Rs）。物質(zhì)點(diǎn)對之間的鍵未斷裂，且各物質(zhì)點(diǎn)鄰域內(nèi)的其余鍵未出現(xiàn)斷裂，即各物質(zhì)點(diǎn)損傷值為零，則三維和二維情況下的標(biāo)量有效力密度狀態(tài)如式（12）所示。

（2）物質(zhì)點(diǎn)對處于過渡域（Rt）。物質(zhì)點(diǎn)對之間的鍵出現(xiàn)斷裂，而各物質(zhì)點(diǎn)鄰域內(nèi)的其余鍵未出現(xiàn)斷裂。此時，式（12）中的物質(zhì)點(diǎn)與物質(zhì)點(diǎn)之間的相互作用力項(xiàng)就會消失，但是由于物質(zhì)點(diǎn)對斷裂之后孔隙空間是連續(xù)的，因此孔隙壓力項(xiàng)是存在的，則三維和二維情況下的標(biāo)量有效力密度狀態(tài)表示為

（3）物質(zhì)點(diǎn)對處于裂隙域（Rf）。裂隙空間的產(chǎn)生需要滿足2 個條件：一是物質(zhì)點(diǎn)之間的鍵的能量密度超過了臨界能量密度，換句話說就是物質(zhì)點(diǎn)對之間的鍵斷裂；二是各物質(zhì)點(diǎn)之間的損傷值超過了臨界損傷值。通過研究統(tǒng)計可知，當(dāng)物質(zhì)點(diǎn)損傷值超過0.4時，裂隙開始生成。故而，三維和二維情況下的標(biāo)量有效力密度狀態(tài)表示為

2.2 流體模塊計算

用達(dá)西定律來描述飽和孔隙裂隙介質(zhì)中的流場，則可表示為

式中：v為體積流速；μf為流體黏度；P為流體壓力。對于均勻各項(xiàng)同性體，K=k I。

在經(jīng)典理論中，流體計算節(jié)點(diǎn)只與相鄰節(jié)點(diǎn)相互作用，且有限體積單元穩(wěn)態(tài)達(dá)西流遵循流量守恒方程。因此，控制方程如式（19）所示：

式中：ρf為流體密度；v為體積流速；S為源項(xiàng)；?為散度算子。

對于裂隙滲流，采用立方定律來計算裂縫域的滲透率［27-28］

式中：b為裂隙開度。如圖1所示，根據(jù)有限體積法，網(wǎng)格節(jié)點(diǎn)的局部通量守恒方程可以離散為

圖1 影響域示意Fig.1 Schematic diagram of area influenced

式中：aP、aW、aS、aE、aN分別為網(wǎng)格節(jié)點(diǎn)與相鄰網(wǎng)格節(jié)點(diǎn)的離散方程系數(shù)；φP、φw、φS、φE、φN分別為網(wǎng)格節(jié)點(diǎn)與相鄰網(wǎng)格節(jié)點(diǎn)的特征通量；sP為網(wǎng)格節(jié)點(diǎn)源項(xiàng)數(shù)值。aP=aW+aS+aE+aN。

對所有網(wǎng)格節(jié)點(diǎn)均可列出對應(yīng)如式（21）的離散方程，對計算域邊界處的離散方程按邊界條件修正各個系數(shù)，再將其寫成矩陣形式，可得到

2.3 耦合計算方案

如圖4所示，在平面離散中，固體部分由近場動力學(xué)物質(zhì)點(diǎn)離散，流體部分由結(jié)構(gòu)性網(wǎng)格離散?？紤]流體壓力和固體變形的相互影響，需要在2 個離散網(wǎng)格之間建立過渡層，從而通過過渡層實(shí)現(xiàn)流體壓力與固體變形之間的相互傳遞。

圖4 耦合過程中的雙向影響示意Fig.4 Schematic diagram of bidirectional influence in coupling process

如圖4 所示，當(dāng)流體壓力向固體部分傳遞時，過渡層由近場動力學(xué)物質(zhì)點(diǎn)控制，通過物質(zhì)點(diǎn)與流體部分網(wǎng)格節(jié)點(diǎn)位置來確定物質(zhì)點(diǎn)所受到的流體壓力由哪一個網(wǎng)格節(jié)點(diǎn)控制。當(dāng)固體部分變形向流體網(wǎng)格傳遞時，此時過渡層由流體部分網(wǎng)格節(jié)點(diǎn)所控制，如流體壓力向固體部分傳遞一般，通過確定流體部分網(wǎng)格節(jié)點(diǎn)控制體積內(nèi)所包含的近場動力學(xué)物質(zhì)點(diǎn)，利用這些物質(zhì)點(diǎn)間的滲透系數(shù)向流體網(wǎng)格傳遞固體變形。算法流程圖如圖5所示。

圖5 算法流程Fig.5 Flowchart of algorithm

3 流固耦合模擬應(yīng)用與驗(yàn)證

通過3個不同的數(shù)值模擬算例說明本文所提出的近場動力學(xué)與有限體積法耦合方法適用于模擬巖石材料流固耦合及裂紋擴(kuò)展破裂過程，揭示相應(yīng)的裂紋擴(kuò)展演化力學(xué)機(jī)制。同時，將該耦合方法的數(shù)值模擬結(jié)果與其他數(shù)值或解析結(jié)果相對比，進(jìn)一步說明所提出的耦合方法能夠有效模擬巖石流固耦合的力學(xué)擴(kuò)展破裂過程。

3.1 多孔介質(zhì)滲流模擬

為了驗(yàn)證近場動力學(xué)與有限體積法耦合方法能對基質(zhì)巖石中的滲流問題進(jìn)行有效的數(shù)值模擬，將基質(zhì)巖石視作多孔介質(zhì)，對其滲流過程進(jìn)行了模擬對比。二維多孔巖石介質(zhì)模型的幾何尺寸如圖6所示，將試件簡化為一定厚度的矩形，其幾何尺寸L、W分別為1.0m、0.2m。在二維多孔巖石介質(zhì)模型的左右兩側(cè)分別施加恒定水壓，即PL0=9.5MPa，PR0=4.5MPa。另外多孔巖石介質(zhì)的滲透率為kr=1.0×10-12m2。

圖6 模型的幾何尺寸Fig.6 Geometric dimensions of model

在固體近場動力學(xué)數(shù)值模型中，采用8 000個近場動力學(xué)物質(zhì)點(diǎn)表征相應(yīng)的多孔巖石介質(zhì)的計算域。其中，相鄰物質(zhì)點(diǎn)的間距為Δx=0.000 5 m，鄰域半徑δ=3.015Δx，材料的彈性模量為E=22GPa，密度為ρ=2 462kg·m-3，在此次模擬中可不用設(shè)置虛擬物質(zhì)點(diǎn)作為固體層的邊界區(qū)域。在流體層設(shè)置8 000個網(wǎng)格表征流體的計算域，將固體層中的物質(zhì)點(diǎn)與流體層中的網(wǎng)格節(jié)點(diǎn)一一對應(yīng)，另流體密度為ρ=1 000kg·m-3，流體黏度系數(shù)為μ=1.0×10-3Pa·s。

多孔巖石介質(zhì)滲流條件下的孔隙水壓分布如圖7 所示?？梢悦黠@從圖中看出，多孔介質(zhì)的水壓力從左右2 個邊界向中心匯聚，且?guī)r石左側(cè)水壓高于巖石右側(cè)水壓。為了驗(yàn)證模擬結(jié)果的正確性，將多孔巖石介質(zhì)滲流模擬結(jié)果中的水壓分布數(shù)值模擬結(jié)果與解析解［29］對比，如圖8 所示?？梢杂^察到數(shù)值模擬結(jié)果與已有的解析解結(jié)果基本吻合，且通過與用近場動力學(xué)方法進(jìn)行滲流模擬的計算時間［8］對比，可以看出近場動力學(xué)與有限體積法耦合方法能夠高效進(jìn)行巖石材料相關(guān)滲流問題的模擬。

圖7 模擬結(jié)果云圖Fig.7 Cloud chart of simulation results

圖8 沿中心軸水壓分布PD-FVM數(shù)值計算結(jié)果與解析解對比Fig.8 Comparison of PD-FVM numerical calculation results and analytical solutions of water pressure distribution along central axis

3.2 流體驅(qū)動的裂縫擴(kuò)展模擬

3.2.1 預(yù)制水平裂隙試樣

為了驗(yàn)證近場動力學(xué)與有限體積法耦合方法能對基質(zhì)巖石中的滲流問題進(jìn)行有效的數(shù)值模擬，將對含有預(yù)制水平裂隙的巖石基質(zhì)中的滲流過程進(jìn)行了模擬驗(yàn)證。二維含水平預(yù)制裂隙的巖石介質(zhì)模型的幾何尺寸如圖9 所示，將試件簡化為一定厚度的矩形，其幾何尺寸L、W分別為1.0m、1.0m。在該巖石介質(zhì)模型的中間位置預(yù)制一個水平裂隙，該預(yù)制裂隙長度為0.2m。將巖石介質(zhì)模型的4 個邊界均設(shè)為排水邊界，即P0=0MPa，并且在預(yù)制水平裂隙內(nèi)部施加水壓p=1×104Pa·s-1。另外多孔巖石介質(zhì)的滲透率為kr=1.0×10-12m2，而裂隙中的等效滲透率為kf=1.333×10-6m2。

圖9 含有預(yù)制水平裂隙模型的幾何尺寸Fig.9 Geometric dimensions of model with a preexisting horizontal flaw

在固體近場動力學(xué)數(shù)值模型中，采用10 000個近場動力學(xué)物質(zhì)點(diǎn)表征相應(yīng)的多孔巖石介質(zhì)的計算域。其中，相鄰物質(zhì)點(diǎn)間距為Δx=0.000 5 m，鄰域半徑δ=3.015Δx，材料的彈性模量為E=22GPa，密度為ρ=2 462kg·m-3，泊松比ν=0.25。在此次模擬中可不用設(shè)置虛擬物質(zhì)點(diǎn)作為固體層的邊界區(qū)域，且各物質(zhì)點(diǎn)的滲透率各個方向均相同。在流體層設(shè)置10 000 個網(wǎng)格表征流體的計算域，將固體層中的物質(zhì)點(diǎn)與流體層中的網(wǎng)格節(jié)點(diǎn)一一對應(yīng)，另流體密度ρ=1 000kg·m-3，流體黏度系數(shù)μ=1.0×10-3Pa·s。因?yàn)榧僭O(shè)各物質(zhì)點(diǎn)的滲透率各向同性，故各網(wǎng)格節(jié)點(diǎn)的滲透率也各向相同。

用解析解對數(shù)值結(jié)果進(jìn)行驗(yàn)證。Sneddon 和Lowengrub［30］給出了在y=0 平面上長度為2lc的單個預(yù)制裂隙，在平面應(yīng)變假設(shè)下，y方向壓力驅(qū)動的位移為

圖10 數(shù)值模擬結(jié)果與解析解沿裂紋的垂直位移對比Fig.10 Comparison of vertical displacement along crack of numerical simulation results and analytical solution

3.2.2 預(yù)制傾斜裂隙試樣

圖11 含有預(yù)制水平裂隙試樣在不同時間步下裂紋擴(kuò)展與水壓分布Fig.11 Crack propagation and water pressure distribution of specimen containing a pre-existing horizontal flaw at different time steps

為了驗(yàn)證近場動力學(xué)與有限體積法耦合方法能對基質(zhì)巖石中的滲流問題進(jìn)行有效的數(shù)值模擬，對含有預(yù)制傾斜裂隙的巖石基質(zhì)中的滲流過程進(jìn)行了模擬驗(yàn)證。二維含預(yù)制傾斜裂隙的巖石介質(zhì)模型的幾何尺寸如圖12 所示，將試件簡化為一定厚度的矩形，其幾何尺寸L、W分別為0.15m、0.15m。在該巖石介質(zhì)模型的中間位置預(yù)制一個角度為θ=30°的傾斜裂隙，該預(yù)制裂隙長度為0.015m。將巖石介質(zhì)模型的4 個邊界均設(shè)為排水邊界，即P0=0MPa，并且在預(yù)制傾斜裂隙內(nèi)部施加水壓p=1×104Pa·s-1。另外多孔巖石介質(zhì)的滲透率為kr=1.0×10-12m2，而裂隙中的等效滲透率為kf=1.333×10-6m2。

圖12 含有預(yù)制傾斜裂隙模型的幾何尺寸Fig.12 Geometric dimensions of model with a preexisting inclined flaw

在固體近場動力學(xué)數(shù)值模型中，采用10 000 個近場動力學(xué)物質(zhì)點(diǎn)表征相應(yīng)的多孔巖石介質(zhì)的計算域。其中，相鄰物質(zhì)點(diǎn)的間距為Δx=0.001 5 m，鄰域半徑δ=3.015Δx，材料的彈性模量為E=15GPa，密度為ρ=2 685kg·m-3，泊松比ν=0.15。在此次模擬中可不用設(shè)置虛擬物質(zhì)點(diǎn)作為固體層的邊界區(qū)域。在流體層設(shè)置10 000個網(wǎng)格表征流體的計算域，將固體層中的物質(zhì)點(diǎn)與流體層中的網(wǎng)格節(jié)點(diǎn)一一對應(yīng)，另流體密度為ρ=1 000kg·m-3，流體黏度系數(shù)為μ=1.0×10-3Pa·s。為了能更好地體現(xiàn)傾斜裂隙的滲透特性，此時需要考慮各物質(zhì)點(diǎn)的滲透率為各向異性，故各網(wǎng)格節(jié)點(diǎn)的滲透率也各向異性。

含有預(yù)制傾斜裂隙的巖石介質(zhì)在注水條件下當(dāng)施加壓力達(dá)到9.006MPa 時預(yù)制裂紋起裂，則裂紋擴(kuò)展過程中的裂紋分布和壓力分布的變化如圖13所示。可以明顯從圖中看出，巖石介質(zhì)中的水壓力從注水孔中心向四周擴(kuò)散，且?guī)r石注水孔處的水壓高于四周水壓。綜合來看，與預(yù)制水平裂隙模擬結(jié)果相似，新生的水力裂紋從預(yù)制裂隙尖端起裂，隨著水壓增大，裂紋沿著預(yù)制裂隙方向擴(kuò)展，且裂紋寬度無明顯變化。同時，最大水壓力存在于裂紋內(nèi)部，水壓在沿裂紋擴(kuò)展方向的壓力梯度遠(yuǎn)小于垂直于裂紋擴(kuò)展方向的壓力梯度。可以將裂紋內(nèi)部的水壓力視作均勻分布。通過與已有的水力壓裂裂紋擴(kuò)展試驗(yàn)［30］對比，模擬計算中裂紋的起裂壓力9.006MPa與試驗(yàn)中試件起裂壓力8.908MPa 的結(jié)果基本接近，且模擬的裂紋擴(kuò)展形態(tài)與文獻(xiàn)中的試驗(yàn)結(jié)果基本吻合，如圖14 所示，說明了近場動力學(xué)與有限體積法耦合方法在考慮滲透率各向異性的情況下，能夠更加貼合實(shí)際地對巖土工程中的水力壓裂問題進(jìn)行高效模擬。

圖13 含有預(yù)制傾斜裂隙試樣在不同時間步下裂紋擴(kuò)展與水壓分布Fig.13 Crack propagation and water pressure distribution of specimen containing a pre-existing inclined flaw at different time steps

圖14 實(shí)驗(yàn)室預(yù)制30°傾斜裂隙水力壓裂試驗(yàn)結(jié)果［31］Fig.14 Hydraulic fracturing test results of laboratory prefabricated 30°inclined fracture[31]

依據(jù)現(xiàn)有參數(shù)開展了不同圍壓對水力壓裂裂紋擴(kuò)展的影響研究。一共設(shè)置了5 種工況，具體工況情況如表1所示。

表1 不同工況的圍壓設(shè)置情況Tab.1 Confining pressure setting under different conditions

模擬結(jié)果如圖15～19 所示，對比5 個工況可以知道，圍壓的存在會改變裂隙起裂壓力，而圍壓差的存在對水力壓裂裂縫擴(kuò)展具有誘導(dǎo)作用，裂紋總是沿著最大主應(yīng)力方向延伸，這一現(xiàn)象與試驗(yàn)結(jié)果幾乎吻合。對比起裂壓力，可以看出，當(dāng)圍壓差相同的情況下，圍壓越大，巖石的起裂壓力就越大；圍壓越小，則巖石越容易被裂隙內(nèi)水壓力壓裂。

圖15 工況1裂紋擴(kuò)展和水壓分布Fig.15 Crack propagation and water pressure distribution under Condition 1

圖16 工況2裂紋擴(kuò)展和水壓分布Fig.16 Crack propagation and water pressure distribution under Condition 2

圖17 工況3裂紋擴(kuò)展和水壓分布Fig.17 Crack propagation and water pressure distribution under Condition 3

圖18 工況4裂紋擴(kuò)展和水壓分布Fig.18 Crack propagation and water pressure distribution under Condition 4

圖19 工況5裂紋擴(kuò)展和水壓分布Fig.19 Crack propagation and water pressure distribution under Condition 5

4 結(jié)語

針對在近場動力學(xué)統(tǒng)一框架下進(jìn)行流固耦合模擬存在的計算效率問題，開展了普通態(tài)型近場動力學(xué)與有限體積法耦合模擬分析方法研究，建立了一種用于模擬飽和裂隙孔隙介質(zhì)中水力裂縫擴(kuò)展的混合有限體積法和普通態(tài)型近場動力學(xué)建模分析方法。普通態(tài)型近場動力學(xué)用來描述固體骨架的變形和捕捉裂紋的發(fā)展，而有限體積法則用來以一種計算高效的方式描述基于達(dá)西定律的流體滲流問題。

首先，通過多孔介質(zhì)滲流模擬算例驗(yàn)證了所提方法的有效性，其數(shù)值結(jié)果與解析解吻合較好，且計算效率大幅提高，能夠很好改善近場動力學(xué)進(jìn)行流固耦合模擬存在的計算效率問題。

其次，通過幾個數(shù)值算例進(jìn)一步驗(yàn)證了該方法在流體驅(qū)動下模擬飽和裂隙孔隙介質(zhì)中裂紋擴(kuò)展的能力和主要特點(diǎn)，在流體驅(qū)動的水力壓裂模擬中觀察到的現(xiàn)象與實(shí)驗(yàn)結(jié)果一致。

最后，通過對5個工況的模擬，展現(xiàn)了圍壓差對水力壓裂裂縫擴(kuò)展具有誘導(dǎo)作用的現(xiàn)象，而且模擬結(jié)果表明裂紋總是沿著最大主應(yīng)力方向延伸，這一現(xiàn)象與試驗(yàn)結(jié)果吻合良好。

作者貢獻(xiàn)聲明：

李卓徽：程序研發(fā)，模型構(gòu)建、數(shù)據(jù)分析及論文撰寫。

周宗青：研究思路指導(dǎo)，論文審閱與修改。

高成路：研究內(nèi)容及方案指導(dǎo)。

張道生：提供編程幫助。

白松松：提供理論幫助。