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

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

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

0 引 言

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

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

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

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

1 理論方法及程序驗證

1.1 近場動力學理論

1.1.1 基本理論概述

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

圖1 物質(zhì)點間相互作用

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

(1)

式中：u為物質(zhì)點位移。

一個物質(zhì)點的運動狀態(tài)是外力和其近場范圍內(nèi)所有物質(zhì)點共同作用的結(jié)果，則根據(jù)牛頓第二定律可得到質(zhì)點x的運動方程為

(2)

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

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

(3)

式中：δ為近場域半徑。

質(zhì)點的位置矢量為：

ξ=x′-x

(4)

t時刻質(zhì)點的位移矢量為：

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

(5)

則質(zhì)點之間力的表達式可簡化為：

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

(6)

為了區(qū)分力的方向并且將其推廣到物體內(nèi)的任意質(zhì)點，可定義質(zhì)點間的力為：

(7)

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

(8)

鍵的微觀彈性應(yīng)變能密度為：

(9)

鍵的伸長率為：

(10)

所以可得基于鍵的近場動力學模型的本構(gòu)函數(shù)：

(11)

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

1.1.2 本構(gòu)關(guān)系

與傳統(tǒng)理論類似，PD 本構(gòu)關(guān)系也是物質(zhì)力的狀態(tài)和變形狀態(tài)之間的關(guān)系，并且PD本構(gòu)關(guān)系中材料參數(shù)可以由傳統(tǒng)本構(gòu)關(guān)系中材料參數(shù)推導(dǎo)得出。在二維平面應(yīng)力問題中，PD本構(gòu)函數(shù)中微觀模量c與傳統(tǒng)本構(gòu)函數(shù)中楊氏模量E存在關(guān)系式：

(12)

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

在傳統(tǒng)本構(gòu)模型中，對于彈脆性材料，當應(yīng)力到達峰值時材料突然破壞，因此可將PD本構(gòu)函數(shù)中的特征函數(shù)μ定義為：

(13)

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

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

(14)

式中：G0為物質(zhì)點的能量釋放率。

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

(15)

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

1.2 隨機場理論

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

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

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

(16)

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

(17)

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

L×LT=C

(18)

(19)

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

采用的二維指數(shù)型自相關(guān)函數(shù)，表達式如下所示：

(20)

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

1.3 程序驗證

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

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

表1 巖體試樣力學參數(shù)

圖3 近場動力學模擬結(jié)果

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

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

1.4 計算流程

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

圖5 程序計算流程圖

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

(2)根據(jù)離散化模型，創(chuàng)建隨機場單元網(wǎng)格；

(3)生成不同統(tǒng)計參數(shù)下的巖體參數(shù)的對數(shù)正態(tài)分布隨機場；

(4)將生成的隨機場作為部分輸入?yún)?shù)導(dǎo)入PD計算程序，確定邊界條件；

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

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

2 結(jié)果討論

2.1 預(yù)制裂縫傾角的影響

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

圖6 均質(zhì)下不同角度預(yù)制裂縫模型裂紋發(fā)展

圖7 空間變異性下不同角度預(yù)制裂縫模型裂紋發(fā)展

圖8 預(yù)制裂縫傾角對裂紋發(fā)展的影響

2.2 變異系數(shù)的影響

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

圖9 變異系數(shù)對裂紋發(fā)展的影響

2.3 水平波動范圍的影響

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

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

2.4 豎向波動范圍的影響

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

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

3 結(jié) 論

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

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

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

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