崔 昊，鄭 宏，李春光，韓 月

（1. 中國電建集團華東勘測設(shè)計研究院有限公司，浙江 杭州 311122；2. 北京工業(yè)大學城市建設(shè)學部，北京 100124；3. 中國科學院武漢巖土力學研究所，湖北 武漢 430071）

外部荷載作用下巖體內(nèi)裂紋的擴展與貫通往往導致巖石工程發(fā)生失穩(wěn)破壞，因此預測巖體中裂紋的擴展過程對于揭示巖體的變形和破壞規(guī)律以及評價巖土工程的安全性具有重要意義。目前，受困于研究手段的不足，對巖體裂紋擴展規(guī)律的研究仍集中于二維空間內(nèi)。但實際工程中的巖體裂隙，無論深埋于巖體內(nèi)部或只露于淺層表面，其本質(zhì)都是三維裂紋擴展問題，因此，對巖體三維裂隙擴展的研究更具有實際意義［1］。

預測三維裂紋的擴展是一項極具難度的問題。在理論研究方面，由于三維裂紋擴展問題的復雜性、多樣性與求解上的巨大難度，至今仍未有理論上的相關(guān)突破。在室內(nèi)試驗方面，一般采用類巖石材料或CT 掃描方式對三維裂隙巖體的破裂機理進行研究［2］，但由于巖體內(nèi)部裂紋擴展觀測手段的匱乏，室內(nèi)試驗尚無法滿足工程問題的需要。

由于理論分析與室內(nèi)試驗的缺乏，學術(shù)界對三維裂紋擴展的研究集中于數(shù)值方法。有限元法（FEM）是最早用于模擬裂紋擴展問題的數(shù)值方法，但FEM在模擬裂紋擴展時需要不斷地重新劃分網(wǎng)格，極大地增加了計算工作量。對于復雜的三維裂紋問題，網(wǎng)格的切割劃分更具挑戰(zhàn)性。鑒于FEM的這些缺陷，出現(xiàn)了基于點近似的無網(wǎng)格方法［3］。該方法不需要隨裂紋擴展進行網(wǎng)格重構(gòu)工作，降低了裂紋擴展問題的難度。但無網(wǎng)格方法中也引入了一些新的缺陷［4］，如：①本質(zhì)邊界條件的施加相對復雜；②計算量更大；③配點型無網(wǎng)格計算穩(wěn)定性差等。

鑒于傳統(tǒng)數(shù)值求解方法在處理裂紋擴展問題時的不足，Silling［5］提出了基于非局部作用思想的近場動力學（Peridynamic，PD）方法。該方法以非局部作用的積分模型代替?zhèn)鹘y(tǒng)理論的微分模型，避免了因在位移不連續(xù)處求導引起的奇異性問題。2007 年Silling 等［6］提出了更為廣泛的非常規(guī)態(tài)基PD（Nonordinary state-based peridynamic，NOSB-PD）理論，該理論借助變形梯度張量F在PD 框架下引入經(jīng)典力學中的應力與應變，搭建了微觀鍵與宏觀強度準則間的關(guān)系，因而在巖土工程領(lǐng)域得到了廣泛的應用［7-8］。目前，NOSB-PD 方法中仍存在因零能模式與邊界效應等問題導致的計算精度不足的缺陷，且至今仍未有較為一致的解決方案［9］。

為解決PD方法中的零能模式與邊界誤差問題，首先簡介NOSB-PD 模型的基本理論及應用無網(wǎng)格伽遼金格式對該模型的重構(gòu)；隨后，引入PD 微分算子（PDDO）近似，并對比其與RKPM近似間的異同；基于對比結(jié)論，建立RKPM-PD耦合算法，并提出適用于該耦合算法的基于增量格式Newmark 算法的隱式迭代流程；最后，采用RKPM-PD耦合算法對若干三維裂紋擴展問題進行模擬分析。

1 近場動力學理論

1.1 非常規(guī)態(tài)基近場動力學

在鍵基PD 以及常規(guī)態(tài)基PD 模型中并沒有傳統(tǒng)意義上應力、應變的概念，這為近場動力學的應用帶來諸多不便。為解決這一問題，一種新的非常規(guī)態(tài)基近場動力學（Non-ordinary state-based peridynamic，NOSB-PD）方法被提出。在引入應力、應變概念后，不同材料的本構(gòu)關(guān)系可以很方便地應用于近場動力學方法。在物質(zhì)點xi處，NOSB-PD方法的運動方程為

為將傳統(tǒng)本構(gòu)模型應用于近場動力學框架，NOSB-PD方法通過變形梯度張量F建立起PD理論與傳統(tǒng)力學間的聯(lián)系。在NOSB-PD 模型中，變形梯度張量F被定義為

為充分應用NOSB-PD 方法與傳統(tǒng)本構(gòu)模型相耦合的優(yōu)勢，可在NOSB-PD 方法中通過物質(zhì)點的應力來定義鍵的損傷。定義鍵ξij的最大主應力為鍵的2 個端點最大主應力的平均值，當鍵上的最大主應力超過材料容許的最大強度時，該鍵發(fā)生斷裂。對于巖石材料，可通過最大拉應力準則及摩爾-庫侖準則來確定巖石材料的最大容許強度［10］。

1.2 基于伽遼金格式的NOSB-PD方法重構(gòu)

對于傳統(tǒng)無網(wǎng)格伽遼金格式，在采用節(jié)點積分后，其最終方程可表示為

式中：M為質(zhì)量矩陣；Kˉ為剛度矩陣；P為外荷載向量；U為所有節(jié)點的自由度向量?；谧冃翁荻葟埩縁的求解公式（4），節(jié)點應變求解矩陣Bi可表示為

式中：g=[gx gy]T=K ξij；m為xi支撐域內(nèi)節(jié)點xj的數(shù)目。

將式（7）代入式（6），可得到式（6）與NOSB-PD方法的全局方程式（3）完全等價的結(jié)論［11］。該結(jié)果揭示了NOSB-PD 方法實質(zhì)上是一種伽遼金弱形式方法?；谠摻Y(jié)論對NOSB-PD 方法中存在的零能模式與邊界誤差進行分析可得：①與其他無網(wǎng)格方法類似，NOSB-PD中的零能模式問題也是由于其在伽遼金弱形式下采用了節(jié)點積分；②NOSB-PD方法中的節(jié)點應變ε是支撐域內(nèi)最小二乘意義下的計算結(jié)果，而最小二乘下的近似一般不會滿足δ屬性，因此，NOSB-PD方法會在其邊界處出現(xiàn)較大的計算誤差。

2 近場動力學微分算子近似與重構(gòu)核函數(shù)近似

在非常規(guī)態(tài)基近場動力學中變形梯度張量F求解方式的基礎(chǔ)上，Madenci等［12］提出了一種被稱為近場動力學微分算子（PDDO）的理論，該理論同樣繼承了PD中以積分代替微分的思想，可通過積分形式求得高階微分算子，可視為式（4）的擴展形式。

2.1 近場動力學微分算子理論

依據(jù)多元函數(shù)泰勒公式，在d維空間中，標量函數(shù)uj在點xj處展開為

其中

在保證矩陣A可逆的情況下，通過式（11）即可計算出物質(zhì)點xi上任意維度任意階導數(shù)的PDDO近似。

2.2 二維空間下的PDDO近似

為便于理解，給出二維空間下一階導數(shù)的PDDO近似形式。在二維線性基條件下，式（9）可簡化為

式中：N，x，N，y分別為位移一階導數(shù)在PDDO 方法下的形函數(shù)。

2.3 RKPM近似函數(shù)

基于Liu等［13］研究成果，RKPM近似位移ui可寫為

對比PDDO 方案得到的位移近似函數(shù)式（14），可以看到PDDO 與RKPM 兩方案得到的位移近似函數(shù)完全一致。

盡管PDDO 與RKPM 兩方案得到的位移近似函數(shù)完全一致，但兩者對位移導數(shù)近似的求解并不一致。RKPM方法中位移導數(shù)的近似函數(shù)通過對位移近似函數(shù)求導得到，即

2.4 PDDO近似與RKPM近似的性質(zhì)對比

Krongauz 與Belytschko［14］指出，為獲得最優(yōu)的數(shù)值精度，伽遼金格式中的近似函數(shù)應滿足一致性條件與積分約束條件。

2.4.1 一致性條件

一般而言，一致性條件總是比較容易滿足的。容易證得，PDDO 與RKPM 均滿足一致性條件。在一維算例中，11 個固定節(jié)點均勻布置在區(qū)間［0，10］上，插值點間距0.01，基向量為p=[1x]T，影響函數(shù)取為三次樣條函數(shù)。假設(shè)節(jié)點位移為y=x，通過插值點計算得到的位移及位移導數(shù)的圖像如圖1所示。從圖中可以看出，當基函數(shù)一致時，PDDO 與RKPM計算得到的位移近似與位移導數(shù)近似具備相同階數(shù)的一致性。

圖1 一階導數(shù)為常數(shù)時RKPM與PDDO近似的精度比較Fig.1 Comparison of approximate accuracy between RKPM and PDDO at y=x

2.4.2 積分約束條件

積分約束條件也被稱作相容性（compatibility）條件，可寫為

式中：Ψ為伽遼金方法中的位移形函數(shù)；Ω為求解域；Γ為求解域邊界。對于全局伽遼金格式，積分約束條件是強制性滿足的。在RKPM 近似中，位移導數(shù)的形函數(shù)是通過對位移形函數(shù)N進行求導運算而得，因此RKPM 近似天然滿足積分約束條件。而PDDO 近似中位移導數(shù)形函數(shù)N，x與位移形函數(shù)N并不一定具有嚴格的求導關(guān)系，前者可能只是形函數(shù)導數(shù)的一種近似逼近，兩者間并不具備嚴格的導數(shù)與原函數(shù)的關(guān)系。因此，一般情況下，積分約束條件式（20）未必滿足。

為充分說明上述條件，取一維桿件進行分析。桿長L=1，節(jié)點間距dx=0.1，彈性模量E=1，一維桿件內(nèi)有均勻分布的應力σ，為保持平衡，桿件左右兩端需施加t=σ的外力荷載。顯然，模型的平衡方程等價為式（20）。模型內(nèi)力記為fint，并加以上標PDDO 及RKPM 以示區(qū)別；模型所受外力記為fout。計算結(jié)果如圖2所示，可以看到，fout=fRKPMint≠fPDDOint，即PDDO近似不滿足積分約束條件，而RKPM近似滿足該條件。

圖2 一維桿件內(nèi)外力分布Fig.2 Distribution of internal and external forces of one-dimensional bar

2.5 不同積分格式下RKPM與PDDO方法計算精度對比

為探究RKPM 與PDDO 方法在不同積分方案下的表現(xiàn)，采用與2.4.2 節(jié)相同的體力作用下的一維桿模型，并分別采用節(jié)點積分方案及兩點高斯積分方案對兩近似方法的收斂性進行分析。在不同積分方案下，兩方法的位移相對誤差收斂性如圖3 所示。對比不同積分方案，可以看到：

圖3 不同積分方案下PDDO與RKPM方法誤差對比Fig.3 Error comparison of PDDO and RKPM methods in different integration schemes

（1）對于RKPM 方法而言，采用高斯積分的精度遠遠高于節(jié)點積分。而對于PDDO 方法而言，積分階次的提高對結(jié)果的影響很小。

（2）在不同積分方案下，RKPM方法的精度均比PDDO方法的精度高。但具體而言，在節(jié)點積分時，兩者的精度比較接近；在采用高斯積分的情況下，兩者精度間的差距最為明顯。

3 RKPM-PD耦合算法

3.1 概述

傳統(tǒng)的無網(wǎng)格方法大多基于斷裂力學準則處理裂紋擴展問題。在數(shù)值計算中，裂紋的追蹤是非常棘手的問題，特別是在三維裂紋擴展問題中，對于裂紋面尖端的動態(tài)追蹤幾乎不可能實現(xiàn)。在PD 方法中，裂紋不再需要顯式的追蹤，而是通過鍵上的損傷來表示，當鍵的臨界伸長率超過一定極限時，兩點之間的作用鍵發(fā)生斷裂。在NOSB-PD 方法中，也可認為當鍵上的應力超過允許應力時，兩點之間的作用鍵發(fā)生斷裂。

基于PDDO 近似與RKPM 近似間的異同以及PD鍵形式的裂紋擴展方案，提出一種RKPM-PD耦合算法，即非裂紋計算部分采用基于RKPM 近似的無網(wǎng)格伽遼金格式以獲取更精確的計算精度，裂紋處理部分則采用PD 中基于鍵形式的損傷模型避免復雜裂紋面的追蹤問題。

在RKPM-PD 中，鍵的斷裂準則與裂紋表示方法與NOSB-PD理論一致，其不同之處在于，RKPMPD中除了更新節(jié)點與節(jié)點間鍵的狀態(tài)外，還需要更新積分點與節(jié)點間鍵的狀態(tài)。具體而言，在節(jié)點積分中只需考慮節(jié)點與節(jié)點間鍵的狀態(tài)，但對于背景網(wǎng)格積分（高斯積分），模型中存在積分點與節(jié)點兩類不同的物質(zhì)點。按照PD鍵的概念，算法中高斯點與節(jié)點間的鍵、節(jié)點與節(jié)點間的鍵，這2類不同的鍵均需考慮其斷裂與否。

3.2 增量格式Newmark算法

在時間迭代中，采用增量格式的Newmark 算法進行計算。t時刻的Newmark方程可寫為

式中：b0、b1、b2均為Newmark 迭代參數(shù)。在計算出ΔU后，即可得到Ut+Δt的值，并通過Ut+Δt計算出新的σt+Δt。之后，根據(jù)不同的斷裂準則，利用t+Δt時刻狀態(tài)信息判斷是否有新的鍵斷掉。對于是否有新鍵斷裂，可分為：

（1）如果當前步中沒有新的斷鍵出現(xiàn)，說明t+Δt時刻方程達到平衡，可進入下一步計算。在通過ΔU得到t+Δt時刻的速度及加速度后，進行變量回代，開啟下一時步的計算。

（2）如果當前步中有新的斷裂出現(xiàn)，說明t+Δt狀態(tài)并未達到平衡，方程（21）仍需進行迭代計算。此時，對于支撐域內(nèi)出現(xiàn)新斷鍵的物質(zhì)點，更新其鍵的狀態(tài)，并重新計算該點上的形函數(shù)[N]及對應的包含一階導數(shù)的矩陣[B]。在新的斷裂出現(xiàn)后，體系內(nèi)的不平衡力R為

根據(jù)式（26）的Ut+Δt重新計算出σt+Δt，并繼續(xù)判斷是否有新的斷鍵出現(xiàn)，如果依然出現(xiàn)新的斷鍵，則返回步驟（2）中繼續(xù)迭代，直至系統(tǒng)內(nèi)不再出現(xiàn)新的斷鍵。此時，結(jié)束當前時步內(nèi)的循環(huán)，轉(zhuǎn)至步驟（1）中。

4 三維裂紋擴展模擬計算

4.1 含傾斜預制裂紋的三維受拉平板

為體現(xiàn)RKPM-PD方法在處理三維裂紋擴展問題中的優(yōu)勢，考慮如圖4 所示模型，其幾何尺寸為2m×1m×0.2m，模型下表面（y=0m）固定，上表面（y=2m）受垂直于該平面的拉伸載荷。模型左側(cè)被一裂紋面完全貫穿，裂紋深度為0.2m，并與Oxy平面呈45°夾角。模型節(jié)點間距為0.03m，支撐域半徑為2 倍節(jié)點間距。模型材料力學參數(shù)為：彈性模量E=30GPa，泊松比ν=0.2，密度ρ=2 500 kg·m-3。裂紋擴展采用臨界伸長率準則，鍵的臨界伸長率sc=0.001。計算采用Newmark 積分，時間步長Δt=0.1ms，載荷Δσ=1×104Pa。

圖4 三維含傾斜裂紋受拉平板示意Fig.4 Sketch map of 3D tensile plate with inclined crack

圖5中給出預制傾斜裂紋面在受拉平板內(nèi)的擴展過程?？梢姡S著拉應力的增大，裂紋首先在傾斜裂紋面的邊緣向前擴展，隨后，裂紋面從傾斜狀態(tài)逐漸扭轉(zhuǎn)為豎直狀態(tài)；之后，裂紋面依然保持垂直于Oxy面的豎直狀態(tài)沿x方向向前延伸，最終貫穿整個平板。從圖中可見，裂紋穿透平板位置約在y=1處，該值近似于傾斜裂紋面y坐標的中心值。另外，在圖6 中給出了對應狀態(tài)下平板內(nèi)y方向位移的分布，可見，位移分布與裂紋面的擴展狀態(tài)相吻合。

圖5 含傾斜裂紋平板在拉伸載荷作用下的裂紋擴展路徑Fig.5 Propagations of 3D plate with inclined crack at tensile load

圖6 含傾斜裂紋平板在拉伸載荷作用下的縱向位移分布Fig.6 Snapshots of vertical displacement of 3D plate with inclined crack

4.2 含預制裂紋三維巴西圓盤

考慮不同厚度的帶預制裂紋的三維巴西圓盤試驗。如圖7 所示，模型在Oxy平面內(nèi)的直徑為100mm，z方向的厚度為10mm（厚度方向為垂直于紙面的方向）。圓盤上下兩端受壓，采用Newmark隱式積分方案，時間步長取為Δt=5×10-7s，平均加載速度為每步10-8m。巖石試樣力學參數(shù)為：彈性模量E=10GPa，泊松比v=0.25，密度ρ=2 500 kg·m-3，抗拉強度T=0.5MPa，黏 聚力C=5MPa 以及內(nèi)摩擦角?=40o。離散模型中平均節(jié)點間距h=2mm。在本算例中，預制裂紋長度為30mm，與水平方向呈45o夾角。為探究預制裂紋深度對結(jié)果的影響，預制裂紋深度分別取10mm、6mm及2mm。

圖7 帶預制裂紋的三維巴西圓盤示意Fig.7 Sketch map of 3B Brazilian disk with prefabricated cracks

當預制裂紋深度為10mm 時，預制裂紋在厚度方向完全貫穿平板，模型在z方向完全一致。該算例的裂紋擴展路徑如圖8所示，可以看到，裂紋擴展路徑與二維情況下相一致，裂紋尖端分別向圓盤上下兩端的壓縮載荷施加處擴展，直至最后裂紋完全貫穿整個圓盤。如圖9 所示，該模型采用RKPM-PD的計算結(jié)果與在三維NMM方法下的模擬結(jié)果［10］及試驗結(jié)果［15］均保持一致。

圖8 三維巴西圓盤含完全貫穿裂紋的裂紋擴展路徑Fig.8 Snapshots of crack path within 3D Brazilian disk with a fully penetrating crack

圖9 含完全貫穿裂紋的三維巴西圓盤在NMM方法及室內(nèi)試驗下的裂紋擴展路徑Fig.9 Crack path of 3D Brazilian disk with a fully penetrating crack obtained by 3D NMM and experiment results

通過對裂紋完全貫穿試件的分析，證明該方法在三維裂紋問題中的有效性。然后考慮裂紋未完全貫穿的情況，設(shè)置裂紋嵌入深度分別為6mm、2mm，其結(jié)果分別顯示于圖10 與圖11。為顯示裂紋未完全貫穿的情況，對于同一時刻的裂紋擴展結(jié)果分別給出其正面圖（z=0）與背面圖（z=10mm）。其中，裂紋從正面嵌入，終止于圓盤內(nèi)部，未穿透背面。

圖10 三維巴西圓盤含嵌入深度6mm未貫穿裂紋的裂紋擴展路徑Fig.10 Snapshots of crack path within 3D Brazilian disk with a penetrating crack (insert depth=6mm)

圖11 三維巴西圓盤含嵌入深度2mm未貫穿裂紋的裂紋擴展路徑Fig.11 Snapshots of crack path within 3D Brazilian disk with a penetrating crack (insert depth=2mm)

通過對未穿透裂紋的結(jié)果分析可知，裂紋首先沿厚度方向擴展，從背面可逐步看到損傷區(qū)域的形成；之后，正面與背面的裂紋均沿豎直方向朝兩端擴展。對比2 個面的裂紋最終擴展軌跡，可以看到兩者大體一致，細微處有所差別。同時可以看到，裂紋切割深度越淺，啟裂所需計算步越多。圖12給出了3 種不同切割深度下圓盤最頂端3 列節(jié)點的軸向平均應力-應變曲線，從圖中可以清晰看出，巴西圓盤在表面裂紋（2mm）時所能承受的載荷遠遠大于深度裂紋（6mm）及完全貫穿裂紋（10mm）時的載荷。

圖12 三維巴西圓盤算例中裂紋擴展過程的應力-應變曲線Fig.12 Axial stress-strain curves of crack propagation in 3D Brazilian disk at different cut depths

5 結(jié)語

針對NOSB-PD 方法面臨的零能模式與邊界誤差問題，提出了一種具備更高精度的RKPM-PD 耦合算法，建立了一套適用于該耦合算法的隱式求解方案，并成功將其應用于三維裂紋擴展問題中。主要結(jié)論如下：

（1）NOSB-PD方法實質(zhì)上等價于采用節(jié)點積分的伽遼金弱形式方法。NOSB-PD 中的零能模式產(chǎn)生的原因是其在弱形式下采用了節(jié)點積分形式。

（2）PDDO 近似與RKPM 近似均滿足一致性條件，但僅有RKPM 近似滿足相容性條件。在提高積分計算階次后，基于RKPM 近似函數(shù)的伽遼金方法計算精度有較為明顯的提升，但基于PDDO 近似函數(shù)的伽遼金方法計算精度并無明顯提高。

（3）RKPM-PD 耦合算法不受零能模式與邊界效應的困擾，具備更高的計算精度，并可有效求解三維動態(tài)裂紋擴展問題。

作者貢獻聲明：

崔 昊：具體研究工作的開展和論文撰寫。

鄭 宏：論文的選題、指導和修改。

李春光：公式推導部分的指導。

韓 月：論文修改。