含缺陷空間變異性巖體裂紋擴展的近場動力學模擬*

仉文崗 孟凡勝 盧志堂 何昌杰 李建新 文海家

(①重慶大學土木工程學院， 重慶 400045， 中國)(②重慶大學山地城鎮建設與新技術教育部重點實驗室， 重慶 400045， 中國)(③重慶大學庫區環境地質災害防治國家地方聯合工程研究中心， 重慶 400045， 中國)(④合肥工業大學資源與環境工程學院， 合肥 230009， 中國)(⑤中國建筑第五工程局有限公司， 長沙 410004， 中國)(⑥中建五局第三建設有限公司， 長沙 410004， 中國)

0 引 言

三峽庫區分布著大量的巖質邊坡，邊坡巖體結構復雜，內部含有裂紋孔隙等原始缺陷(王成等， 2004； 湯明高等， 2019； 周家文等， 2019)。巖體結構面作為巖體的重要組成部分，其存在不僅破壞了巖體完整性，而且對巖體的強度和變形破壞具有很大影響(蘭恒星等， 2019； 趙海軍等， 2019； Chen et al.，2020)。巖體破壞的外因是外部荷載擾動導致巖體屈服，內因則是內部裂紋的萌生、擴展、貫通導致巖體承載力下降。因此對于含缺陷巖體的裂紋發展過程，從細觀力學角度對其受力變形和損傷演化過程進行研究，將其與巖石宏觀斷裂力學機制聯系起來，是當今巖石力學與工程領域的重點問題，對于庫岸巖體邊坡穩定分析也具有十分重要的現實意義(趙程等， 2015)。

長期以來，許多學者對含缺陷巖體的損傷破壞進行了研究，其中包括理論分析，模型試驗，數值模擬等方法。目前理論發展尚不完全成熟，模型試驗成本較高且含缺陷試樣制作復雜，因此數值模擬方法得到了更廣泛的應用(蔣明鏡等， 2014； 王洪建等， 2015； 陳鵬宇， 2018)。而在數值模擬方法中，有限元法會遇到裂紋尖端奇異性的問題，擴展有限元法在模擬不連續變形時需要引入外部準則，離散元法在模擬連續區域時不能反映連續變化的特點。由于含缺陷巖體材料兼有連續區域與非連續區域，使得大多數數值模擬方法難以滿足要求。Silling教授提出的近場動力學方法(Peridynamics， PD)較好地解決了這個問題(Silling， 2000)。近場動力學理論基于非局部作用思想，采用空間積分進行物質內部作用的描述。它兼有無網格方法和分子動力學的優點，相比其他方法，該方法對于微觀與宏觀、連續與非連續的力學行為具有統一的表述，在數值方面具有無網格屬性和不連續求解功能。因此對于兼有連續性與非連續性的含缺陷巖體材料，該方法具有很大的適用性。Yang et al. (2020)利用基于狀態的近場動力學模型研究了含非直線預制裂隙的紅砂巖試樣在不同裂隙角下的斷裂行為。馬鵬飛等(2019)基于改進的近場動力學模型，對含單裂縫試樣在拉伸荷載條件下的裂紋發展過程進行了模擬。李錚等(2019)通過5個不同的數值算例說明了近場動力學算法在處理脆性巖體材料斷裂問題方面的有效性和準確性。盧志堂等(2016)提出了近場動力學和有限差分的聯合算法，并通過對比證明了該方法在模擬巖體層裂破壞的優越性。然而在利用近場動力學方法模擬巖體裂紋發展的文獻中，考慮巖體材料空間變異性的還幾乎沒有。

巖土體在形成的歷史過程中由于地質環境及物理化學作用的不同，在不同的空間位置處會表現出空間變異性。大量研究表明巖土體空間變異性對邊坡，基坑，隧道等各類建筑物的穩定性有重大影響。如Liu et al. (2019)采用蒙特卡洛模擬結合極限平衡法和材料點法研究了巖土體空間變異性邊坡變形過程中多種破壞模式的發生和演化過程，并且強調了真實巖土體空間變異性數據的價值。黃廣龍等(2010)研究發現采用可靠度指標評價基坑整體穩定性比安全系數更為合理，并且可靠度指標受黏聚力、內摩擦角及支護樁嵌固深度變異性的影響較大。王長虹等(2018)將隨機場理論引入盾構隧道地表沉降的可靠度指標分析中，發現與隨機變量模型相比，基于隨機場理論的模型可以得到更精確的可靠指標以優化隧道設計。最近也有學者開始考慮空間變異性對巖體自身變形損傷方面的影響，提出了一種結合隨機場理論和離散元方法的概率隨機離散元分析(RDEA)方法，研究單軸壓縮下巖石碎裂的機理(Zhao et al.， 2020)。但總體來說，考慮空間變異性對巖體失效及裂紋發展影響的研究還很少。

本文采用了一種隨機近場動力學方法(RPD)來模擬單軸壓縮下含缺陷巖體裂紋的發展。用近場動力學方法模擬巖體裂紋發展過程，隨機場方法表征巖體材料的空間變異性。以二維含預制單裂縫矩形巖體試樣為例，對計算程序進行了準確性驗證，并比較了不同預制裂縫角度后續裂紋的發展情況。最后通過隨機場統計參數的變化，重點研究了空間變異性對后續裂紋發展的影響，所得結論可為庫區裂隙巖體穩定性分析提供一定參考價值。

1 理論方法及程序驗證

1.1 近場動力學理論

1.1.1 基本理論概述

2000年，由美國Sandia實驗室的Silling教授提出了近場動力學方法的基本思想。近場動力學經過多年的發展已形成以鍵為基礎、普通狀態為基礎與非普通狀態為基礎等多種分支(朱其志等， 2016)，本文僅考慮以鍵為基礎的近場動力學。如圖1所示，某物體占據空間域R，假設某一時刻t，空間內任一物質點x與其周圍空間一定半徑范圍內的其他物質點x′之間存在著相互作用，這種相互作用稱為鍵，即質點之間的力通過鍵傳遞，若相互作用力為f，則

圖1 物質點間相互作用

f=f(x，x′，u(x，t)，u(x′，t)t)

(1)

式中：u為物質點位移。

一個物質點的運動狀態是外力和其近場范圍內所有物質點共同作用的結果，則根據牛頓第二定律可得到質點x的運動方程為

(2)

式中：ρ為物質密度；b為單位體積物質所受的外載荷；Hx為質點x的近場域范圍，可用下式表示：

Hx=H(x，δ)={x′-x≤δ|x′∈δ}

(3)

式中：δ為近場域半徑。

質點的位置矢量為：

ξ=x′-x

(4)

t時刻質點的位移矢量為：

η=u′(x′，t)-u(x，t)

(5)

則質點之間力的表達式可簡化為：

f=f(x，x′，u(x，t)，u(x′，t)，t)=f(η，ξ)

(6)

為了區分力的方向并且將其推廣到物體內的任意質點，可定義質點間的力為：

(7)

對于線彈性材料，f(η，ξ)為鍵的微觀彈性應變能密度ω對該鍵相對位移矢量η的導數：

(8)

鍵的微觀彈性應變能密度為：

(9)

鍵的伸長率為：

(10)

所以可得基于鍵的近場動力學模型的本構函數：

(11)

式中：c為微觀模量；μ為表征鍵的連接狀態的特征函數； |η+ξ|與|ξ|分別為鍵變形后的長度與初始長度。

1.1.2 本構關系

與傳統理論類似，PD 本構關系也是物質力的狀態和變形狀態之間的關系，并且PD本構關系中材料參數可以由傳統本構關系中材料參數推導得出。在二維平面應力問題中，PD本構函數中微觀模量c與傳統本構函數中楊氏模量E存在關系式：

(12)

需要注意的是在PD中，材料的泊松比被限定為ν=1/3(Huang et al.，2015)。

在傳統本構模型中，對于彈脆性材料，當應力到達峰值時材料突然破壞，因此可將PD本構函數中的特征函數μ定義為：

(13)

式中：s0為極限伸長率，即鍵斷裂的臨界變形值，當鍵的變形超過臨界變形值，鍵斷裂且不再傳遞荷載。

對于二維平面應力問題，鍵的極限伸長率為：

(14)

式中：G0為物質點的能量釋放率。

因此通過鍵的斷裂情形，可以定義物質點的非局部損傷參數，即表征物質點的損傷情況：

(15)

可見近場動力學方法本身就存在破壞準則，不需要再引入其他準則來判斷材料是否發生破壞。

1.2 隨機場理論

Vanmarcke(1997)率先將隨機場方法引入巖土工程可靠度領域。隨機場是一個描述空間變異性的概率模型，可以用均值、標準差、自相關函數和波動范圍等予以描述。大量研究表明，空間任意兩點巖土體間存在著自相關性和變異性。由于巖土工程中現場試驗數據的有限性，基于這些數據難以建立起表征巖土體參數自相關性的自相關函數，因此一般選擇采用理論自相關函數來近似代替真實自相關函數。另一方面，巖土體參數大多數遵循非高斯分布，并且基于喬列斯基分解的中點法在模擬非高斯隨機場中易于理解，便于編程實現。

因此本文采用基于喬列斯基分解的中點法模擬相關非高斯隨機場(蔣水華等， 2014)，對數正態變量m的隨機場，表達式如下所示：

Rm(x，y)=exp[μlnm+σlnm·Gm(x，y)]

(16)

式中：μlnm和σlnm分別為變量lnm的均值和標準差；Gm(x，y)為正態分布隨機場，可通過下式實現：

(17)

式中：n為參數隨機場數目；Zk為獨立標準正態隨機樣本矩陣；L為喬列斯基分解計算的三角矩陣，如下所示：

L×LT=C

(18)

(19)

式中：C為相關矩陣，τxij=|xi-yj| 和τyij=|yi-yj|分別為空間任意兩點間的水平和垂直方向相對距離。

采用的二維指數型自相關函數，表達式如下所示：

(20)

式中：δh和δv分別為水平和豎直方向的波動范圍，波動范圍越大，表示參數的空間自相關程度越強。由式(20)可知自相關函數ρ(τx，τy)只與空間兩點間的相對距離有關，與兩點絕對位置無關。

1.3 程序驗證

本文利用MATLAB進行編程，對Wang et al. (2017)中含單裂縫模型在單軸壓縮條件下的裂紋擴展過程進行了模擬。如圖2所示，矩形試樣高150mm，寬75mm，試樣中心包含一長為12.7mm，寬為1.27mm的缺陷，缺陷是通過刪除該處材料點來實現的。將模型簡化為二維平面應力問題，并進行離散化處理。離散近場動力學材料點間距為dx=0.83mm，影響域半徑為3.015dx，在試樣上下端部設置限制側向位移的端部約束條件，然后施加速度邊界條件，通過改變最外層物質點的速度，將邊界條件通過鍵的作用逐漸傳遞給材料區域內的所有物質點，采用自適應動態松弛方法進行時間積分，從而實現單軸壓縮的模擬過程。假定巖體試樣為彈脆性材料，各參數取表1材料均值，當預制裂縫傾角為45°時，模擬結果如圖3所示。首先在預制裂縫尖端出現損傷，然后裂紋朝向兩端發展，而由于在上下端部存在約束，使得裂紋最終沒有朝向最大壓應力方向擴展，但總體來說模擬結果與原文獻基本一致，包括起裂角度，初始擴展路徑等。因此認為本文采用的近場動力學計算程序在模擬巖體裂紋發展方面是適用的。

圖2 單軸壓縮下含單一缺陷巖體試樣

表1 巖體試樣力學參數

圖3 近場動力學模擬結果

為了模擬巖體試件的空間變異性，基于隨機場理論，利用MATLAB編程，將彈性模量E和能量釋放率G0定義為隨機場變量，統計特性參數如表1所示。當水平波動范圍δh分別為0.10m和1.00m時，彈性模量E的隨機場分布如圖4所示，可見本文所用隨機場程序可以較好實現巖體試樣的空間變異性。

圖4 不同水平波動范圍下E的隨機場分布

1.4 計算流程

考慮空間變異性存在時巖體損傷情況計算流程如下，做出程序計算流程圖如圖5所示：

圖5 程序計算流程圖

(1)確定物質點間距dx與近場范圍δ，模型離散化處理，確定巖體參數統計特性；

(2)根據離散化模型，創建隨機場單元網格；

(3)生成不同統計參數下的巖體參數的對數正態分布隨機場；

(4)將生成的隨機場作為部分輸入參數導入PD計算程序，確定邊界條件；

(5)計算粒子間作用力，判斷鍵是否斷裂；

(6)計算每時步的所有粒子間作用力，計算所有時步，得出結果。

2 結果討論

2.1 預制裂縫傾角的影響

為了比較不同預制裂縫傾角下后續裂紋的發展情況，選擇預制裂縫傾角為30°， 45°， 60°時，分別對考慮隨機場存在與不存在的情況進行了模擬，隨機場參數如表1所示。不同角度下均質巖體裂紋發展圖像如圖6所示，考慮空間變異性時裂紋發展圖像如圖7所示。選擇預制裂縫尖端一側相鄰3行3列共9個物質點，監測各物質點在各時步下的位移，做出9個物質點的位移均值與時步關系曲線如圖8所示，通過位移時程曲線的比較說明不同情況下的裂紋發展速度。結合圖6與圖7對比發現考慮空間變異性時兩側裂紋發展不再完全對稱且模型整體破壞更為嚴重。圖8則表明不管是否考慮巖體的空間變異性，都是預制裂縫傾角為45°時后續裂紋發展略大于預制裂縫傾角為60°時，當預制裂縫為30°時后續裂紋發展則相對較慢。當然上述裂紋發展模擬結果除了與巖體參數有關，也會受端部約束等設置因素的影響。并且通過對比可以發現空間變異性的存在會減弱預制裂縫傾角的影響，空間變異性存在時不同預制裂縫角度下監測點最終位移差比不考慮空間變異性時更小。圖8結果也表明，不管預制裂縫角度如何變化，考慮空間變異性的存在時后續裂紋的發展都更快，說明若把巖體當做均質材料時將低估巖體的受損傷程度。

圖6 均質下不同角度預制裂縫模型裂紋發展

圖7 空間變異性下不同角度預制裂縫模型裂紋發展

圖8 預制裂縫傾角對裂紋發展的影響

2.2 變異系數的影響

為了考慮E和G0的變異程度對后續裂紋發展的影響，在隨機場模型中將E和G0的變異系數COV分別設置為0.1， 0.2， 0.3， 0.4， 0.5，水平波動范圍δh和豎向波動范圍δv分別取0.20m和0.10m，預制裂縫傾角選擇45°，預制裂縫尖端附近物質點的位移時程曲線如圖9所示。結果表明E和G0的變異系數對后續裂紋的發展速率有顯著影響，變異系數越高，預制裂縫尖端附近物質點位移越大，即后續裂紋發展更快。因為隨著變異系數的增大，參數變異性增強，同時巖體更多部分會表現出參數E和G0數值較小的情況，因此隨著變異系數的增大，裂紋發展速度明顯提高。由圖9可知1000時步時COV為0.5時比COV為0.1時物質點位移增長了92.796%。

圖9 變異系數對裂紋發展的影響

2.3 水平波動范圍的影響

圖10所示為隨機場水平波動范圍δh對后續裂紋發展的影響。Ching et al.(2013)研究表明，對于指數型自相關函數，隨機場單元尺寸小于0.018～0.054倍的波動范圍時模擬效果較好。考慮到隨機場單元及模型整體尺寸，水平波動范圍δh分別選擇0.10m， 0.15m， 0.20m， 0.60m， 1.00m進行模擬，E和G0的變異系數取0.3，豎向波動范圍δv取0.10m，預制裂縫傾角為45°，物質點位移時程曲線結果表明：整體來看，隨著水平波動范圍δh的增大，后續裂紋發展更快，原因歸結為當δv一定時，δh越大，越接近為軟硬交互的水平層狀巖體，層理弱面相對而言剛度小、變形大，加載方向與層理弱面垂直，此時層理弱面對軸向壓縮變形影響較大，因此隨著δh的增大裂紋發展更快，這與鄧華鋒等(2020)的研究結果是一致的。但是δh從0.10m增長到0.20m的影響要大于從0.20m增長到1.00m的影響，原因則可能是囿于模型尺寸的限制，δh較大時對隨機場模型影響已不明顯。

圖10 水平波動范圍對裂紋發展的影響

2.4 豎向波動范圍的影響

圖11所示為隨機場豎向波動范圍δv對后續裂紋發展的影響。同樣基于隨機場單元及模型整體尺寸的考慮，豎向波動范圍δv分別選取0.05m， 0.08m， 0.10m， 0.3m， 0.5m，E和G0的變異系數為0.3，水平波動范圍δh取0.2m，在預制裂縫傾角為45°模型中，位移時程曲線結果表明：整體來看隨著δv的增大后續裂紋發展稍快，但是在不同δv下物質點位移區別并不大，δv對后續裂紋發展影響有限。原因一方面是模型尺寸的影響，另一方面則是δh一定時，δv越大越接近為垂直層理巖體，加載方向與層理弱面平行，層理弱面分布對軸向變形影響較小，因此表現為在不同δv下，裂紋發展速度區別不大。

圖11 豎向波動范圍對裂紋發展的影響

3 結 論

本文采用了一種新的隨機近場動力學方法(RPD)模擬研究了單軸壓縮下含缺陷巖體的損傷情況。用近場動力學方法(PD)模擬巖體材料受壓損傷過程中裂紋發展過程，用隨機場方法(RFM)表征巖體材料彈性模量E和能量釋放率G0的空間變異性，并利用已有模型驗證了計算程序的準確性，得出主要結論如下：

(1)近場動力學方法可以有效模擬巖體材料損傷及裂紋發展問題，當考慮巖體材料彈性模量E和能量釋放率G0的空間變異性時，裂紋發展速度明顯加快。

(2)預制裂縫傾角對巖體損傷有較大影響，在任何預制裂縫角度下考慮空間變異性存在時，后續裂紋發展都更迅速，并且當預制裂縫為45°時，巖體最容易受壓破壞。

(3)隨機場參數中，E和G0的變異系數和水平波動范圍δh都對裂紋發展速度有明顯影響。隨著變異系數的增大，水平波動范圍的增大，后續裂紋發展速度都有明顯提高，而在不同豎向波動范圍δv下，裂紋發展速度沒有太大區別。因此確定合理的隨機場參數值特別是變異系數及水平波動范圍，對確定單軸壓縮下巖體強度及損傷至關重要。