高成路，李術(shù)才，周宗青，2，李利平，張道生，厲明浩，劉冠男

（1.山東大學(xué)齊魯交通學(xué)院，山東濟(jì)南 250002；2.陸軍勤務(wù)學(xué)院軍事設(shè)施系，重慶 401311；3.深部巖土力學(xué)與地下工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇徐州 221116；4.中國(guó)礦業(yè)大學(xué)力學(xué)與土木工程學(xué)院，江蘇徐州 221116）

地下水是導(dǎo)致工程巖體失穩(wěn)破壞的主要原因之一[1］，地下水滲流與地應(yīng)力耦合作用下工程巖體的變形、損傷、破壞和失穩(wěn)機(jī)制是眾多研究人員關(guān)心的課題，尤其是巖體水力壓裂的干預(yù)和控制成為當(dāng)前研究熱點(diǎn)。巖體水力壓裂是指由于水壓升高引起巖體中裂隙發(fā)生與擴(kuò)展的物理現(xiàn)象，在提高油氣產(chǎn)量方面具有重要的應(yīng)用價(jià)值[2］，同時(shí)，水力壓裂法也是國(guó)際巖石力學(xué)學(xué)會(huì)建議的地應(yīng)力測(cè)量方法之一。[3］。在巖土工程領(lǐng)域，巖體水力壓裂又作為一種地質(zhì)災(zāi)害廣泛存在，如大壩漏水、隧道突涌水等，往往造成嚴(yán)重的傷亡事故。因此，科學(xué)描述巖體水力壓裂過(guò)程、揭示巖體應(yīng)力-滲流耦合破壞機(jī)制，對(duì)于油氣開(kāi)采和災(zāi)害防控等均具有重要意義。

針對(duì)巖體水力壓裂機(jī)理，國(guó)內(nèi)外學(xué)者開(kāi)展了大量力學(xué)分析，基于彈性理論、損傷理論、斷裂理論和分形理論等先后提出了多種描述方法[4-6］，但是大量的現(xiàn)場(chǎng)實(shí)踐表明，現(xiàn)有理論模型很難對(duì)大型水力壓裂中觀測(cè)的多縫起裂和擴(kuò)展路徑進(jìn)行合理的解釋。在室內(nèi)試驗(yàn)方面，國(guó)內(nèi)外學(xué)者針對(duì)水力壓裂以及裂縫擴(kuò)展規(guī)律開(kāi)展了大量研究，其中真三軸模擬壓裂試驗(yàn)是最為常見(jiàn)的一種形式。陳勉等[7］采用大尺寸真三軸模擬試驗(yàn)系統(tǒng)結(jié)合聲發(fā)射監(jiān)測(cè)裝備對(duì)天然巖樣和人造巖樣進(jìn)行水力壓裂裂縫擴(kuò)展機(jī)理模擬實(shí)驗(yàn)，揭示了地應(yīng)力、斷裂韌性、節(jié)理和天然裂縫等因素對(duì)水力裂縫擴(kuò)展的影響規(guī)律。趙益忠等[8］利用真三軸模擬壓裂實(shí)驗(yàn)系統(tǒng)對(duì)多種巖樣進(jìn)行了水力壓裂裂縫起裂及裂縫擴(kuò)展模擬實(shí)驗(yàn)，得到了壓后裂縫幾何形態(tài)和壓裂過(guò)程中壓力隨時(shí)間的變化規(guī)律。郭印同等[9］、衡帥等[10］、李芷等[11］、侯振坤等[12］采用真三軸巖土工程模型試驗(yàn)機(jī)結(jié)合工業(yè)CT掃描技術(shù)建立了一套頁(yè)巖水力壓裂物理模擬與壓裂縫表征方法，并開(kāi)展了頁(yè)巖水力壓裂物理模擬試驗(yàn)研究，實(shí)現(xiàn)了水力壓裂裂縫的空間形態(tài)描述，并揭示了水力裂縫與天然裂縫的相互作用規(guī)律。雖然真三軸壓裂測(cè)試系統(tǒng)結(jié)合工業(yè)CT掃描、聲發(fā)射裝置、X衍射等先進(jìn)技術(shù)成為研究復(fù)雜裂縫網(wǎng)絡(luò)形成機(jī)理的主要試驗(yàn)手段，但是由于巖石材料內(nèi)部結(jié)構(gòu)和組成成分的復(fù)雜性，想要直觀展示水力裂縫的擴(kuò)展過(guò)程還是比較困難的。

理論、試驗(yàn)和模擬是研究巖石力學(xué)問(wèn)題的三大途徑。隨著數(shù)值模擬理論與方法的豐富和完善，研究人員采用有限元、邊界元和顆粒離散元等方法針對(duì)水力壓裂過(guò)程開(kāi)展了大量研究。謝興華[13］開(kāi)發(fā)了三維水力劈裂有限元計(jì)算程序，研究了開(kāi)裂對(duì)滲透性的影響規(guī)律。薛炳等[14］采用滲流-應(yīng)力耦合模型對(duì)水力壓裂過(guò)程進(jìn)行了三維有限元研究。袁志剛等[15］采用ANSYS軟件模擬了穿層鉆孔水力壓裂裂縫的擴(kuò)展。彪仿俊等[16］使用ABAQUS軟件針對(duì)水力壓裂水平裂縫影響參數(shù)開(kāi)展了數(shù)值模擬研究。Francisco等[17］利用擴(kuò)展有限元模擬了多孔巖石中水力裂縫和天然裂縫的交互作用，揭示了地應(yīng)力差和接近角等參數(shù)在裂縫網(wǎng)絡(luò)發(fā)展中的作用規(guī)律。王濤等[18］利用擴(kuò)展有限元法進(jìn)行了頁(yè)巖水力壓裂裂縫擴(kuò)展模擬。Olson[19］利用邊界元理論建立了水平井多段壓裂的裂縫擴(kuò)展模型。王理想等[20］基于連續(xù)-非連續(xù)單元法和中心型有限體積法提出了解決水力壓裂流固耦合問(wèn)題的二維混合數(shù)值計(jì)算模型。Yang等[21］利用顆粒離散元軟件PFC開(kāi)展了大量水力裂縫擴(kuò)展規(guī)律研究，揭示了射孔角度和射孔布置對(duì)同步壓裂裂縫擴(kuò)展機(jī)理和巖石破壞模式的影響規(guī)律。

雖然傳統(tǒng)數(shù)值方法在模擬巖體水力壓裂過(guò)程中取得了較好的研究結(jié)果，但是巖體水力壓裂過(guò)程是典型的連續(xù)-非連續(xù)應(yīng)力-滲流耦合過(guò)程，傳統(tǒng)連續(xù)變形分析方法和非連續(xù)變形分析方法將巖體材料視為完全連續(xù)體或者完全非連續(xù)體，均與實(shí)際情況有所出入。基于非局部作用思想的近場(chǎng)動(dòng)力學(xué)不再采用連續(xù)變形假設(shè)，而是通過(guò)積分方程求解物理力學(xué)問(wèn)題，突破了傳統(tǒng)方法裂紋擴(kuò)展路徑的限制，成為求解連續(xù)-非連續(xù)問(wèn)題的一種有效方法[22-23］。該方法不僅可以模擬固體材料的損傷演化機(jī)制，還能夠描述地下水在孔隙介質(zhì)、裂隙介質(zhì)中的滲流規(guī)律。本文將近場(chǎng)動(dòng)力學(xué)在模擬固體材料變形與地下水滲流兩方面的優(yōu)勢(shì)相結(jié)合構(gòu)建流體壓力驅(qū)動(dòng)條件下裂隙巖體應(yīng)力-滲流耦合的常規(guī)態(tài)型近場(chǎng)動(dòng)力學(xué)模擬方法，通過(guò)與室內(nèi)試驗(yàn)及傳統(tǒng)數(shù)值方法計(jì)算結(jié)果對(duì)比驗(yàn)證模擬方法的有效性，揭示天然巖體水力壓裂過(guò)程中水力裂隙的擴(kuò)展規(guī)律。

1 近場(chǎng)動(dòng)力學(xué)基本原理

近場(chǎng)動(dòng)力學(xué)通過(guò)將計(jì)算區(qū)域離散為有限數(shù)量帶有物性信息的、具有一定體積的物質(zhì)點(diǎn)，并且利用所有物質(zhì)點(diǎn)空間位置關(guān)系的變化和對(duì)外界擾動(dòng)信息的響應(yīng)來(lái)表征材料的變形與破壞。根據(jù)近場(chǎng)動(dòng)力學(xué)的非局部作用思想（示意圖見(jiàn)圖1），計(jì)算區(qū)域Rx內(nèi)的任意物質(zhì)點(diǎn)i的位置矢量xi與其一定范圍內(nèi)所有物質(zhì)點(diǎn)的位置矢量xj（xj∈Rx，‖ ‖xi-xj≤δ）存在相互作用，其中δ為鄰域（horizon），為一個(gè)正常數(shù)。xi只與其鄰域范圍內(nèi)的物質(zhì)點(diǎn)存在相互作用，而與超出這個(gè)局部區(qū)域的物質(zhì)點(diǎn)之間的相互作用為零。也就是說(shuō)，近場(chǎng)動(dòng)力學(xué)研究的是物質(zhì)點(diǎn)與其鄰域范圍內(nèi)的所有物質(zhì)點(diǎn)相互作用的物理現(xiàn)象，非局部作用思想使近場(chǎng)動(dòng)力學(xué)在模擬裂隙巖體破壞問(wèn)題時(shí)能夠?qū)崿F(xiàn)跨裂紋計(jì)算[24-25］。

圖1 近場(chǎng)動(dòng)力學(xué)非局部作用示意Fig.1 Schematic diagram of nonlocal action of peridynamics

在鍵型近場(chǎng)動(dòng)力學(xué)模型中，物質(zhì)點(diǎn)對(duì)之間存在大小相等、方向相反的力密度矢量。假設(shè)xi與其鄰域內(nèi)任意其他物質(zhì)點(diǎn)的xj之間的相互作用力分別為

式中：T為力密度矢量，N·m-6；u為物質(zhì)點(diǎn)位移矢量，m；t為時(shí)間，s；f為本構(gòu)力函數(shù)（constitutive function），包含材料的所有本構(gòu)信息，不再以傳統(tǒng)的應(yīng)力-應(yīng)變關(guān)系形式出現(xiàn)，不再需要位移場(chǎng)的連續(xù)性假設(shè)及對(duì)位移場(chǎng)進(jìn)行空間求導(dǎo)。因此，無(wú)論位移場(chǎng)的連續(xù)與否，都不再影響基本方程的求解，近場(chǎng)動(dòng)力學(xué)模型在求解不連續(xù)問(wèn)題時(shí)不再出現(xiàn)數(shù)值奇異性等病態(tài)特征。

近場(chǎng)動(dòng)力學(xué)基本方程[26］為

式中：ρ為物質(zhì)點(diǎn)密度，kg·m-3；u?為物質(zhì)點(diǎn)加速度，m·s-2；V為物質(zhì)點(diǎn)體積；m3；H為物質(zhì)點(diǎn)的非局部作用范圍；b為外載荷密度，N·m-3。據(jù)此可得由本構(gòu)力函數(shù)表示的xi的運(yùn)動(dòng)方程，如式（4）：

在近場(chǎng)動(dòng)力學(xué)中，物質(zhì)點(diǎn)鍵的伸長(zhǎng)率s可以由鍵的相對(duì)位置ξ=xj-xi和相對(duì)位移η=uj-ui表示，見(jiàn)式（5）：

在近場(chǎng)動(dòng)力學(xué)本構(gòu)關(guān)系中通過(guò)引入標(biāo)量函數(shù)μ(t，ξ)表征材料的破壞特性，如圖2所示，即當(dāng)鍵的伸長(zhǎng)率s超過(guò)極限值s0時(shí)，則認(rèn)為對(duì)應(yīng)的鍵斷裂，此時(shí)物質(zhì)點(diǎn)間相互作用力消失。

圖2 物質(zhì)點(diǎn)間相互作用力與伸長(zhǎng)率的關(guān)系[27］Fig.2 Relationship between interaction force and stretch of material points[27]

因此，對(duì)于微觀彈脆性（PMB）材料[27］，其本構(gòu)力函數(shù)為

式中：c為微彈性模量，Pa·m-4；μ(t，ξ)為具有歷史依賴性的標(biāo)量值函數(shù)，表示物質(zhì)點(diǎn)鍵的完整性，即當(dāng)μ(t，ξ)=1時(shí)鍵是完整的，而當(dāng)μ(t，ξ)=0時(shí)鍵是斷裂的，其表述為

式中：s0為臨界伸長(zhǎng)率，表征物質(zhì)點(diǎn)鍵所具有的最大伸長(zhǎng)量。需要注意的是，盡量PMB材料在初始狀態(tài)下是各向同性的，但是變形過(guò)程中某些特定方向鍵的斷裂將導(dǎo)致后續(xù)響應(yīng)的各向異性。

據(jù)此，可以通過(guò)統(tǒng)計(jì)物質(zhì)點(diǎn)斷裂鍵的數(shù)量與初始狀態(tài)的比值定義物質(zhì)點(diǎn)的局部損傷，如式（8）：

由此可見(jiàn)，物質(zhì)點(diǎn)的局部損傷是一個(gè)從0到1之間變化的值。當(dāng)φ=1時(shí)，所有與物質(zhì)點(diǎn)相互作用的鍵全部斷裂；當(dāng)φ=0時(shí)，所有與物質(zhì)點(diǎn)相互作用的鍵全部完整。局部損傷代表了材料內(nèi)部裂隙的形成過(guò)程，近場(chǎng)動(dòng)力學(xué)本構(gòu)力函數(shù)中包含了材料的損傷與斷裂描述，在模擬材料的損傷破壞時(shí)，不再需要附加的斷裂準(zhǔn)則。

為了描述材料的損傷程度，將計(jì)算區(qū)域圍巖損傷狀態(tài)與完整狀態(tài)之比定義為損傷指數(shù)[28］，如式（9）：

因此，當(dāng)損傷程度越大時(shí)，χ值越高。

鍵型近場(chǎng)動(dòng)力學(xué)模型為解決材料損傷破壞等不連續(xù)問(wèn)題提供了有效的模擬手段，但是依然存在3個(gè)主要問(wèn)題：一是過(guò)度簡(jiǎn)化，泊松比受到了限制；二是鍵型近場(chǎng)動(dòng)力學(xué)以對(duì)點(diǎn)力描述材料的本構(gòu)特征，無(wú)法與傳統(tǒng)應(yīng)力-應(yīng)變本構(gòu)關(guān)系直接建立聯(lián)系；三是無(wú)法區(qū)分材料的體積變形和形狀變形，因此無(wú)法描述材料的塑性不可壓縮性。因此，Silling等[29］通過(guò)引入“狀態(tài)”的數(shù)學(xué)概念建立態(tài)型近場(chǎng)動(dòng)力學(xué)模型。依據(jù)對(duì)點(diǎn)力是否沿鍵長(zhǎng)方向，態(tài)型近場(chǎng)動(dòng)力學(xué)模型又可以分為常規(guī)態(tài)型和非常規(guī)態(tài)型。

在常規(guī)態(tài)型近場(chǎng)動(dòng)力學(xué)模型中，物質(zhì)點(diǎn)間的相互作用力與鍵長(zhǎng)方向一致，其力密度矢量也滿足角動(dòng)量守恒定律，因此，力密度矢量可定義為

式中：y為變形后物質(zhì)點(diǎn)的位置矢量；A、B為近場(chǎng)動(dòng)力學(xué)參數(shù)，可以由應(yīng)變能密度變分得到，如式（12）、式（13）：

式中：ω為量綱一化影響函數(shù)，控制xi周圍質(zhì)點(diǎn)對(duì)其作用力的大??；θ為體積應(yīng)變；a、b、d為近場(chǎng)動(dòng)力學(xué)參數(shù)，可以由近場(chǎng)動(dòng)力學(xué)應(yīng)力能密度與經(jīng)典連續(xù)介質(zhì)力學(xué)應(yīng)變能密度等效求出。依據(jù)Madenci和Oterkus[30］的研究，不同維度條件下近場(chǎng)動(dòng)力學(xué)參數(shù)值為

式中：h為二維結(jié)構(gòu)厚度，m；A為物質(zhì)點(diǎn)橫截面積，m2；K和G分別為體積模量和剪切模量，Pa。

式中：E為彈性模量，Pa；μ為泊松比；ωij利用量綱分析可設(shè)定其大小為

同樣，利用近場(chǎng)動(dòng)力學(xué)力密度表示柯西應(yīng)力，并得到物質(zhì)點(diǎn)xi體積應(yīng)變的一般形式為

近場(chǎng)動(dòng)力學(xué)不僅可以描述固體材料的損傷破壞，也可以用來(lái)描述地下水滲流。在近場(chǎng)動(dòng)力學(xué)模型中，計(jì)算區(qū)域中任意一個(gè)物質(zhì)點(diǎn)xi與其鄰域內(nèi)其他所有物質(zhì)點(diǎn)xj通過(guò)鍵的形式相互作用，通過(guò)將物質(zhì)點(diǎn)看作是儲(chǔ)存水的容器，而將鍵看作是傳遞水的通道，可以模擬地下水在多孔介質(zhì)中的滲流過(guò)程[31］，其中鍵兩端物質(zhì)點(diǎn)存在的水力勢(shì)能差被認(rèn)為是驅(qū)動(dòng)水運(yùn)輸?shù)膭?dòng)力。在鍵型近場(chǎng)動(dòng)力學(xué)理論中，認(rèn)為每個(gè)鍵上的流體運(yùn)輸是相互獨(dú)立的，與其他鍵不產(chǎn)生相互影響。對(duì)于二維均勻各向同性多孔介質(zhì)的單向流動(dòng)過(guò)程，其近場(chǎng)動(dòng)力學(xué)基本方程為[32］

式中：Φ為物質(zhì)點(diǎn)的有效水力壓強(qiáng)，Pa；μw為流體黏度，Pa·s；k為材料滲透率，m2；ρw為流體密度，kg·m-3；φ為孔隙率；R為源匯項(xiàng)，kg·s-1；-ω為量綱一化影響函數(shù)，代表某個(gè)物質(zhì)點(diǎn)的xi的滲流狀態(tài)受其鄰域內(nèi)其他物質(zhì)點(diǎn)的xj影響程度的大小，取ω-(xi，xj)=1/‖ξxixj‖。在進(jìn)行二維均勻各向同性孔隙介質(zhì)穩(wěn)態(tài)滲流模擬時(shí)，地下水滲流場(chǎng)水力壓強(qiáng)可以表示為

式中：P為流體壓強(qiáng)，Pa；g為重力加速度，m·s-2；z為垂直方向水頭高度，m。對(duì)于恒溫條件下微小恒定壓縮流，流體密度可以表示為

式中：ρ0為流體密度在壓力為P0時(shí)的密度，kg·m-3；γ為流體的壓縮系數(shù)，是描述流體壓縮性大小的物理量，被定義為單位壓力變化時(shí)引起的液體單位體積的變化量，水體壓縮系數(shù)與壓力和溫度有關(guān)，Pa-1。

當(dāng)巖體滲流主要受裂隙控制時(shí)，一般采用離散裂隙網(wǎng)格模型來(lái)模擬巖體的滲流特征。該模型假設(shè)完整巖石的滲透性遠(yuǎn)遠(yuǎn)低于節(jié)理裂隙，因此巖石的滲透性可以忽略，認(rèn)為流體只在網(wǎng)絡(luò)狀的裂隙中運(yùn)動(dòng)，因此該模型比等效連續(xù)介質(zhì)模型更加符合實(shí)際情況。

由于裂隙是巖體滲流的主要通道，而巖體在地質(zhì)構(gòu)造的長(zhǎng)期影響下，裂隙的形狀差異十分嚴(yán)重，并且裂隙具有明顯的方向性，平行于裂隙方向的滲透系數(shù)顯然比垂直于裂隙方向的滲透系數(shù)大得多，因此裂隙巖體滲流具有明顯的各向異性特征。影響裂隙巖體滲流特性的主要因素有：裂隙的分組、空間分布、密度和規(guī)模等；裂隙組及單條裂隙的產(chǎn)狀、大小、粗糙度、連通性和充填特性等；裂隙之間的空間交互關(guān)系；巖塊孔隙率和滲透性；水的流動(dòng)狀態(tài)等[33］。由于這些因素的復(fù)雜性、隨機(jī)性和不確定性，在數(shù)值計(jì)算和工程分析中通常難以建立準(zhǔn)確的表征裂隙滲流特性的關(guān)系式，因此，通常將裂隙簡(jiǎn)單概化為由具有一定長(zhǎng)度、傾向、傾角和寬度的空間離散裂隙網(wǎng)格。立方定律是巖體水力學(xué)的理論基礎(chǔ)，將其寫(xiě)成經(jīng)典達(dá)西定律的形式為

式中：wd為裂隙的寬度，m；ν為流體運(yùn)動(dòng)黏滯系數(shù)，m2·s-1；J為水力梯度；kf為縫隙的滲透率，m2，有

由質(zhì)量守恒定律，根據(jù)式（18）同理可推導(dǎo)出當(dāng)ω-(xi，xj)=1/‖ξxixj‖時(shí)，二維裂隙滲流的近場(chǎng)動(dòng)力學(xué)基本方程為

式中：ρf為裂隙流體密度，即穿越裂隙的流體鍵兩端物質(zhì)點(diǎn)的流體密度，kg·m-3，一般可取為平均值；φf(shuō)為裂隙空間體積與物質(zhì)點(diǎn)體積的比值；Φf為裂隙表面物質(zhì)點(diǎn)的有效水力壓強(qiáng)，Pa；Rf為裂隙源匯項(xiàng)，kg·s-1；μf為裂隙流體黏度，Pa·s。

2 物質(zhì)點(diǎn)雙重覆蓋理論模型

天然巖體中賦存有大量的孔隙和裂隙，不僅導(dǎo)致巖體的力學(xué)性質(zhì)發(fā)生改變，也影響著其滲透特征。受地質(zhì)構(gòu)造的影響，地下工程圍巖不僅處于一定的地應(yīng)力環(huán)境中，同時(shí)還處于一定的滲流場(chǎng)中，應(yīng)力場(chǎng)與滲流場(chǎng)之間存在相互影響，這種相互影響被稱為應(yīng)力-滲流耦合作用（或流-固耦合作用）。裂隙巖體流固耦合作用主要表現(xiàn)為2個(gè)方面：一是地下水在裂隙巖體中滲流時(shí)導(dǎo)致巖體中產(chǎn)生滲流作用力，主要包括靜水壓力和動(dòng)力水壓力，從而改變巖體原始應(yīng)力狀態(tài)；二是裂隙巖體應(yīng)力狀態(tài)的改變導(dǎo)致巖體結(jié)構(gòu)發(fā)生變形甚至破壞，進(jìn)而改變巖體的滲透性能及地下水滲流路徑[34］。

在經(jīng)典土力學(xué)中，Terzaghi（太沙基）提出了著名的飽和土有效應(yīng)力原理，該原理將飽和土體看作由土顆粒和孔隙水2種介質(zhì)組成，土體的變形由土顆粒骨架所承擔(dān)的有效荷載控制[35］。Terzaghi有效應(yīng)力原理在分析飽和土體的力學(xué)特性時(shí)得到了符合客觀規(guī)律的認(rèn)識(shí)，對(duì)促進(jìn)土力學(xué)成為一門獨(dú)立的學(xué)科起到了重要的推動(dòng)作用。然而由于巖體材料的固相是連通的而非散體材料，因此Terzaghi有效應(yīng)力原理不再適用。為此，通過(guò)引入考慮飽和材料體積壓縮變形的修正系數(shù)α可以建立適用于巖體材料的有效應(yīng)力表達(dá)式[36］，如式（24）：

式中：σ為應(yīng)力，Pa；α為Biot系數(shù)，代表了孔隙水壓力對(duì)有效應(yīng)力的貢獻(xiàn)程度。程遠(yuǎn)方等[37］采用Crossplotting法、排水實(shí)驗(yàn)法與聲波動(dòng)態(tài)法測(cè)試了不同滲透性天然巖芯的Biot系數(shù)值，結(jié)果表明：高滲巖芯為0.90～1.00，中滲巖芯為0.80～0.90，而低滲巖芯為0.60～0.75。

利用有效應(yīng)力原理，本文提出的模擬裂隙巖體應(yīng)力-滲流耦合過(guò)程的近場(chǎng)動(dòng)力學(xué)物質(zhì)點(diǎn)雙重覆蓋理論模型如圖3所示。即物質(zhì)點(diǎn)不僅具有表征材料變形的固體力學(xué)性質(zhì)，同時(shí)也是儲(chǔ)水“容器”，地下水通過(guò)物質(zhì)點(diǎn)間的鍵傳遞水壓力。由于孔隙水壓力的作用，物質(zhì)點(diǎn)固體鍵發(fā)生了弱化作用，同時(shí)還需要承受裂隙表面的滲透壓力。與固體鍵的狀態(tài)完全相反，巖體損傷導(dǎo)致固體鍵的斷裂正好激活了流體鍵，此時(shí)，斷裂的鍵具有更高的滲流能力。

圖3 應(yīng)力-滲流耦合模擬的物質(zhì)點(diǎn)雙重覆蓋理論模型[28］Fig.3 Theoretical model of material point double repeated cap for stress-seepage coupling[28]

考慮地下水滲流時(shí)，近場(chǎng)動(dòng)力學(xué)物質(zhì)點(diǎn)間的相互作用力（式（10）、式（11））可以改進(jìn)為

對(duì)于裂隙表面，物質(zhì)點(diǎn)間不再發(fā)生相互作用，但是由于流體滲透壓力的存在，其運(yùn)動(dòng)狀態(tài)依然受到滲透作用力的影響。因此，當(dāng)物質(zhì)點(diǎn)間的鍵斷裂后，其力狀態(tài)可表示為[38］

可以看到，在有效應(yīng)力原理中，裂隙表面在水壓力的作用下相當(dāng)于牽引邊界，從而降低了材料的總應(yīng)力狀態(tài)。

3 數(shù)值迭代求解策略

近場(chǎng)動(dòng)力學(xué)運(yùn)動(dòng)方程和滲流方程的求解均需要借助數(shù)值方法來(lái)實(shí)現(xiàn)，通過(guò)在有限數(shù)量離散的物質(zhì)點(diǎn)上進(jìn)行時(shí)間積分，利用當(dāng)前時(shí)間步的位移增量與滲流增量預(yù)測(cè)下一時(shí)間步物質(zhì)點(diǎn)的位置、變形和孔隙水壓力，據(jù)此判斷物質(zhì)點(diǎn)鍵的損傷與材料的破壞。時(shí)間積分是近場(chǎng)動(dòng)力學(xué)求解的核心，在計(jì)算過(guò)程中驅(qū)動(dòng)著模擬的進(jìn)行和信息的傳遞，與經(jīng)典局部方法相比，近場(chǎng)動(dòng)力學(xué)時(shí)間積分方案并無(wú)顯著差異。

時(shí)間積分分為顯式時(shí)間積分和隱式時(shí)間積分。其中，顯式時(shí)間積分是通過(guò)大量的微小時(shí)間步長(zhǎng)來(lái)確定系統(tǒng)對(duì)模擬過(guò)程的響應(yīng)，通常用于瞬態(tài)模擬，適合于求解材料變形和破壞過(guò)程中的動(dòng)態(tài)問(wèn)題。隱式時(shí)間積分通常用于模擬具有相對(duì)較大的載荷步并且動(dòng)態(tài)影響可以忽略不計(jì)的情況，通過(guò)將加速度項(xiàng)設(shè)為零來(lái)得到準(zhǔn)靜態(tài)模擬的控制方程。

當(dāng)忽略裂隙巖體破壞過(guò)程中的動(dòng)力現(xiàn)象時(shí)，可以將滲透壓力作為驅(qū)動(dòng)載荷，逐級(jí)施加到模型中，對(duì)平衡方程進(jìn)行迭代求解，從而得到系統(tǒng)在當(dāng)前荷載條件下的穩(wěn)定狀態(tài)。在隱式時(shí)間積分方案中，系統(tǒng)的平衡狀態(tài)可以通過(guò)殘差來(lái)確定，當(dāng)殘差值降低到指定閾值時(shí)，準(zhǔn)靜態(tài)加載步被認(rèn)為是收斂的。因此，本文將近場(chǎng)動(dòng)力學(xué)流-固耦合求解過(guò)程劃分為固體模塊和流體模塊分別進(jìn)行，采用“混合”時(shí)間積分方案，即固體模塊采用隱式時(shí)間積分、流體模塊采用顯式時(shí)間積分，其計(jì)算流程如圖4所示。

圖4 近場(chǎng)動(dòng)力學(xué)流-固耦合模擬計(jì)算流程Fig.4 Flowchart of peridynamics fluid-solid coupling simulation

4 水力裂隙萌生-擴(kuò)展過(guò)程模擬

巖石注水試驗(yàn)是水力壓裂研究最常見(jiàn)的手段之一[39］。利用水力壓裂研究中的經(jīng)典試驗(yàn)案例對(duì)近場(chǎng)動(dòng)力學(xué)流-固耦合方法進(jìn)行可靠性分析和有效性驗(yàn)證?？招膱A柱體注水試驗(yàn)是經(jīng)典水力壓裂試驗(yàn)之一，采用平面應(yīng)力模型開(kāi)展二維條件下空心圓環(huán)水力壓裂過(guò)程模擬，如圖5所示。試件半徑r=25mm，其中間預(yù)留一個(gè)半徑r0=2mm的注水孔。模型采用環(huán)形等距方式離散為7 636個(gè)物質(zhì)點(diǎn)，物質(zhì)點(diǎn)間距Δ=0.5mm，鄰域大小δ=1.5mm。

圖5 空心圓柱體水力壓裂數(shù)值模型Fig.5 Numerical model of hydraulic fracturing of hollow cylinder

模型采用均勻、各向同性假設(shè)，巖體密度ρ=2 090kg·m-3，彈性模量E=5GPa，泊松比μ=0.3，臨界伸長(zhǎng)率s0=0.001。由于巖體水力壓裂過(guò)程耗時(shí)較短，巖體破壞過(guò)程中的滲流多以裂隙為主，因此采用裂隙介質(zhì)模型對(duì)材料滲流場(chǎng)進(jìn)行模擬，從而忽略巖石基質(zhì)的滲流作用。巖體材料孔隙率φ=0.5，Biot系數(shù)α=1.0，裂隙流體密度ρf=1 000kg·m-3，黏度μ=0.899Pa·s，壓縮系數(shù)γ=0.048 5。巖體滲透率與裂隙的寬度的平方成正比，即符合立方定律。試驗(yàn)注水速率采用壓力增量的形式進(jìn)行設(shè)置，其大小為ΔP=10kPa·s-1。通過(guò)固體模塊與流體模塊的交替迭代可以模擬空心圓柱體注水過(guò)程中水力裂隙的演化過(guò)程。

空心圓柱體水力壓裂過(guò)程中巖體的損傷、破壞與滲流演化規(guī)律如圖6所示。在注水過(guò)程中，圓孔內(nèi)壁孔隙水壓力逐級(jí)遞增，圓柱體受到水壓力的作用，在周向承受拉應(yīng)力作用，并引起材料發(fā)生拉伸變形。當(dāng)內(nèi)孔水壓力較低時(shí)，試件所受應(yīng)力未超過(guò)其強(qiáng)度，只發(fā)生變形，未出現(xiàn)損傷和破裂。隨著內(nèi)孔水壓力的升高，大約在P=6.6MPa時(shí)，內(nèi)孔周向拉應(yīng)力超過(guò)試件抗拉強(qiáng)度，裂紋由內(nèi)壁產(chǎn)生并向外壁擴(kuò)展，最終形成貫穿裂隙，試件劈裂破壞。

圖6 空心圓柱體水力壓裂模擬結(jié)果Fig.6 Simulation results of hydraulic fracturing of hollow cylinder

在該模擬中巖體的變形和破壞完全是由水力驅(qū)動(dòng)的，水力裂隙的產(chǎn)生是隨機(jī)的，不需要指定裂隙擴(kuò)展路徑。對(duì)比巖體損傷和孔隙水壓力分布云圖可以看出，流體沿著裂隙擴(kuò)展路徑滲流，并且僅在裂隙表面局部范圍內(nèi)進(jìn)行。由于采用裂隙滲流模型，在裂隙尖端孔隙水壓力為零，因此裂隙內(nèi)孔隙水壓力呈梯度分布。裂隙內(nèi)孔隙水壓力的升高，進(jìn)一步驅(qū)動(dòng)著裂隙由試件內(nèi)壁向外壁延伸，并最終形成劈裂破壞。

為了更好地分析空心圓柱體水力壓裂過(guò)程中巖體的損傷破壞程度，建立巖體損傷指數(shù)隨注水壓力的變化關(guān)系，如圖7所示?？梢钥闯觯诘退畨毫l件下（P＜6.5MPa），巖體保持良好的完整性，此時(shí)試件所受應(yīng)力未超過(guò)巖體破壞強(qiáng)度，并沒(méi)有局部損傷的產(chǎn)生。但是當(dāng)注水壓力超過(guò)臨界值時(shí)，巖體迅速損傷破壞。這表明，空心圓柱體在水力壓裂過(guò)程中，致使試件破裂的能量存在積蓄-釋放過(guò)程。材料的損傷和破壞消耗掉了儲(chǔ)存在巖體中的變形能，并隨著注水壓力的增大，進(jìn)一步產(chǎn)生能量的累積。但是由于裂隙的存在，在裂隙尖端出現(xiàn)應(yīng)力集中，導(dǎo)致材料更容易發(fā)生破壞。

圖7 空心圓柱體水力壓裂模擬過(guò)程中損傷指數(shù)變化規(guī)律Fig.7 Variation of damage index in hydraulic fracturing simulation of hollow cylinder

圖8給出了采用超高壓大流量滲流-應(yīng)力耦合試驗(yàn)儀開(kāi)展的空心圓柱體水力壓裂試驗(yàn)結(jié)果[39］以及應(yīng)用巖石破裂失穩(wěn)的滲流應(yīng)力耦合分析系統(tǒng)（FRFPA）開(kāi)展的數(shù)值模擬結(jié)果[40］。試驗(yàn)和數(shù)值模擬結(jié)果表明，注水過(guò)程中水力裂隙首先出現(xiàn)在內(nèi)壁表面，并隨著水壓力的增大，裂隙由內(nèi)壁向外壁擴(kuò)展，并最終形成貫通裂隙，試件整體發(fā)生劈裂破壞。本文利用近場(chǎng)動(dòng)力學(xué)流-固耦合模擬方法得到了相似的研究結(jié)論。但是，由于天然巖體材料的非均勻性，試驗(yàn)過(guò)程中裂隙首先出現(xiàn)在內(nèi)壁面材料強(qiáng)度較弱的某點(diǎn)上，并且由于注水壓力引起的應(yīng)力集中導(dǎo)致材料局部產(chǎn)生高拉應(yīng)力區(qū)，引起試件開(kāi)裂。F-RFPA采用Weibull分布將模型的力學(xué)參數(shù)進(jìn)行隨機(jī)賦值，以此來(lái)模擬巖體材料的非均勻性，得到了單條裂隙貫穿試件表面，而本文建立的近場(chǎng)動(dòng)力學(xué)模型對(duì)稱性較好，且采用了均勻性、各向同性假設(shè)，因此裂隙分布也呈現(xiàn)較好的對(duì)稱性。

圖8 空心圓柱體水力壓裂試驗(yàn)和RFPA模擬結(jié)果Fig.8 Results of hydraulic fracturing test and RFPA simulation for hollow cylinder

5 結(jié)語(yǔ)

基于有效應(yīng)力原理建立近場(chǎng)動(dòng)力學(xué)流-固耦合模擬的物質(zhì)點(diǎn)雙重覆蓋理論模型，并開(kāi)發(fā)了基于固體和流體時(shí)間步循環(huán)迭代的求解器，通過(guò)將水力壓裂動(dòng)態(tài)問(wèn)題轉(zhuǎn)化為準(zhǔn)靜態(tài)問(wèn)題實(shí)現(xiàn)水壓力驅(qū)動(dòng)作用下裂隙巖體水力壓裂過(guò)程的模擬。通過(guò)開(kāi)展裂隙巖體水力壓裂試驗(yàn)?zāi)M研究，并與前人研究結(jié)果進(jìn)行對(duì)比，驗(yàn)證該理論模型的準(zhǔn)確性，并取得了如下結(jié)論：

（1）基于有效應(yīng)力原理，通過(guò)將近場(chǎng)動(dòng)力學(xué)在模擬固體變形與流體滲流兩方面的優(yōu)勢(shì)結(jié)合起來(lái)可以構(gòu)建表達(dá)裂隙巖體流-固耦合效應(yīng)的物質(zhì)點(diǎn)雙重覆蓋理論模型，即物質(zhì)點(diǎn)不僅儲(chǔ)存了關(guān)于固體介質(zhì)變形和損傷的力學(xué)信息，同時(shí)也是儲(chǔ)存流體介質(zhì)的容器，鍵不僅反映物質(zhì)間相互作用力，同時(shí)也是流體滲流的通道。

（2）通過(guò)開(kāi)展空心圓柱體水力壓裂試驗(yàn)?zāi)M，并與室內(nèi)試驗(yàn)及前人研究結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文所建立的裂隙巖體應(yīng)力-滲流耦合模擬方法的準(zhǔn)確性和可靠性。

（3）研究結(jié)果顯示，空心圓柱體水力壓裂過(guò)程中裂紋首先由內(nèi)壁產(chǎn)生并向外壁擴(kuò)展，最終形成貫穿裂隙，試件劈裂破壞，該過(guò)程完全是由水力驅(qū)動(dòng)的，水力裂隙的產(chǎn)生是隨機(jī)的，不需要指定裂隙擴(kuò)展路徑，并且水力壓裂過(guò)程中致使試件破裂的能量存在積蓄-釋放過(guò)程，近場(chǎng)動(dòng)力學(xué)流-固耦合模型可以有效模擬水力裂隙萌生擴(kuò)展的客觀規(guī)律，為油氣開(kāi)采、地下水災(zāi)害防控等提供了分析手段。

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

高成路：程序研發(fā)及論文撰寫(xiě)。

李術(shù)才：研究思路指導(dǎo)。

周宗青：研究?jī)?nèi)容制定。

李利平：研究方案指導(dǎo)。

張道生：數(shù)值建模。

厲明浩：算例驗(yàn)證。

劉冠男：論文修訂。