近場動力學(xué)法頻散特性及其在巖石層裂分析中應(yīng)用

盧志堂， 王志亮, 2

(1.同濟大學(xué) 土木工程學(xué)院，200092 上海； 2.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，230009 合肥)

?

近場動力學(xué)法頻散特性及其在巖石層裂分析中應(yīng)用

盧志堂1， 王志亮1, 2

(1.同濟大學(xué) 土木工程學(xué)院，200092 上海； 2.合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，230009 合肥)

摘要:為了解近場動力學(xué)方法(peridynamics, PD)的計算精度，考察該法用于巖石層裂破壞模擬的效果，對PD進行了頻散分析和算法驗證.首先由頻散分析后發(fā)現(xiàn)：當(dāng)空間步長不變時，隨影響域變大PD法頻散愈嚴重；而空間步長減小時，影響域節(jié)點數(shù)不變，其頻散會變?nèi)?；?dāng)影響域大小不變時，內(nèi)部劃分節(jié)點越密集，頻散越弱.其次，通過該方法與傳統(tǒng)有限差分法的比較表明PD離散方程可看作一系列差分方程的組合，其截斷誤差為影響域半徑δ的二階無窮?。划?dāng)δ為Δx時，PD算法與中心差分法是等價的，且此時計算精度最高.最后，通過PD法應(yīng)用于巖桿一維層裂模擬分析，探討了其空間步長、影響域尺寸對計算結(jié)果的影響，得出層裂時間、層裂位置及損傷分布情況，并與層裂試驗進行對比分析.PD可用于巖石層裂破壞分析，將FDM和PD法兩者結(jié)合進行層裂模擬時，計算時間少、優(yōu)勢明顯.

關(guān)鍵詞:近場動力學(xué)；頻散特性；有限差分法；巖桿；層裂

材料斷裂損傷是工程領(lǐng)域熱點問題之一，除了借助相關(guān)實驗手段外，研究人員也在積極尋求數(shù)值方法上的突破，以期弄清楚材料裂紋產(chǎn)生及其破壞機理，并進而能對斷裂進行預(yù)警.雖然有限元和有限差分等方法在工程分析中得到了廣泛應(yīng)用，但這些傳統(tǒng)數(shù)值方法大多都是基于經(jīng)典連續(xù)介質(zhì)力學(xué)的偏微分方程而來，在解決不連續(xù)問題時會產(chǎn)生奇異值，很難實現(xiàn)對材料內(nèi)部裂縫產(chǎn)生與損傷發(fā)展的模擬.為很好地刻畫裂縫生成和材料損傷演化，革新數(shù)值方法和發(fā)展相關(guān)理論的研究一直未間斷，而近場動力學(xué)(peridynamics，PD)正是其中的理論之一[1].文獻[2]對經(jīng)典連續(xù)介質(zhì)力學(xué)進行了重建，提出了PD 理論，將傳統(tǒng)的偏微分方程替換為積分方程，并可直接對不連續(xù)介質(zhì)力學(xué)問題(比如裂縫擴展與分叉)進行分析.在一定條件下，選擇合適的響應(yīng)函數(shù)，PD方程收斂于經(jīng)典連續(xù)介質(zhì)力學(xué)方程[3].將PD方程離散，結(jié)合損傷模型就可實現(xiàn)對脆性和延性材料的斷裂問題進行模擬.該方法理論清晰、過程簡單，已有模擬結(jié)果與試驗數(shù)據(jù)吻合較好[4-8].在國內(nèi)，文獻[9-10]利用PD理論模擬了混凝土的脆性損傷情況，指出PD法可很好刻畫和模擬混凝土在拉伸情況下的損傷累積與漸進破壞過程.

雖然PD方法能很好反映材料斷裂形式，但是材料斷裂前，通常需要經(jīng)歷一段連續(xù)變形過程，此時該方法的頻散特性還有待深入研究；另一方面，PD的計算效果與基于連續(xù)介質(zhì)力學(xué)發(fā)展起來的數(shù)值算法(如FDM，有限差分法)相比，是否具有優(yōu)勢？能否將其與現(xiàn)有方法聯(lián)合應(yīng)用于材料斷裂分析也值得探討.

1PD法計算原理

1.1基本方程

如圖1所示，材料占據(jù)空間域B，x為B內(nèi)一材料點，Hx為x的影響域，δ為Hx的半徑.假設(shè)x與其影響域Hx內(nèi)的另一材料點x′通過連接鍵產(chǎn)生相互作用，作用力f為

(1)

式中|x′-x|≤δ，x和x′分別為材料點x的位置矢量，u表示位移矢量.

圖1　空間域B與影響域Hx示意

材料點x的PD方程為

(2)

式中:ρ為密度，ü為x點加速度，Vx′為x′點所占體積，b為體力.ξ=x′-x，為相對位置矢量，η=u(x,t)-u(x,t)，為相對位移，見圖2.

相互作用力f可寫為

(3)

圖2　影響域內(nèi)相對位置與相對位移

1.2本構(gòu)關(guān)系與損傷準則

f(η,ξ)可寫成連接鍵伸長率的函數(shù)，即PD本構(gòu)關(guān)系式：

(4)

式中:cm表示微彈模，在一維、二維和三維條件下依次為2E/(Aδ2)、9E/(πδ3)、12E/(πδ4)[11]，E表示彈性模量，A表示材料點占據(jù)的面積，s表示材料點間連接鍵伸長率，其表達式為

(5)

(6)

式中s0為連接鍵臨界伸長率，當(dāng)s ≥ s0時，連接鍵斷裂，且不可恢復(fù).

x點的局部損傷φ可表示為失效的連接作用與全部連接的比值，即

(7)

根據(jù)局部損傷可判斷材料內(nèi)部裂縫的發(fā)育及傳播過程.通常情況下，需定義損傷值φ0作為判斷裂縫是否產(chǎn)生的準則，本文借鑒文獻[4]建議值，即φ0= 0.3.

1.3求解思路

式(2)通常借助數(shù)值積分處理，寫成求和表達式：

(8)

式中：i和p分別表示計算節(jié)點與其影響域內(nèi)的節(jié)點編號，Vp表示點p的在影響域內(nèi)占據(jù)的體積.

圖3給出二維節(jié)點分布，圖中節(jié)點間距均為Δx，每個節(jié)點占據(jù)的面積為(Δx)2.Vp近似[12]為

其中V = (Δx)3.

采用時域中心差分計算加速度：

(10)

式中:上標n表示時步數(shù)，Δt為時間間隔，計算的穩(wěn)定性條件為Δt ≤Δx/c，c表示波速.

PD法數(shù)值計算流程見圖4.

圖3　PD計算的離散域示意

圖4　PD程序流程圖

2頻散特性分析

雖然目前已有研究表明PD在模擬脆性材料破壞過程與試驗吻合較好，但對其在材料破壞前的模擬效果關(guān)注較少.文獻[13]發(fā)現(xiàn)利用PD數(shù)值計算時存在頻散現(xiàn)象，所謂數(shù)值頻散實質(zhì)上是一種因離散化求解波動方程而產(chǎn)生的偽波動[14]，這種頻散不同于波動方程本身引起的物理頻散，而是PD離散方程所固有的本質(zhì)特征.這會導(dǎo)致波在傳播過程中，波前形狀發(fā)生變化，并且逐漸散開，進而引起模擬得出的波速和裂縫擴展速度失真.數(shù)值頻散對研究材料中波傳播過程，尤其是在沖擊荷載(高頻波成分較多，數(shù)值頻散愈嚴重)作用下是不利的.本節(jié)將結(jié)合頻域分析，對PD法的數(shù)值頻散特性進行分析，以求壓制數(shù)值頻散，提高計算精度.以一維方程為例，將式(2)的解寫成一般形式：

(11)

(12)

考慮影響域半徑、空間步長的變化，得出頻散曲線(見圖5，縱坐標表示頻率ω(2π)-1，橫坐標表示波數(shù))，其中參數(shù)按照黑云母花崗巖參數(shù)給出：彈性模量E= 20 GPa，ρ= 2 620 kg·m-3.

圖5　頻散特性影響因素

圖5(a)中空間步長Δx= 0.005 m保持不變，不斷增大影響域，不難發(fā)現(xiàn)δ= Δx時頻散曲線與真實曲線最接近；m越大，得到的曲線與真實曲線相差越大，即影響域增大時，PD方法數(shù)值頻散越嚴重；同時還能看出，頻率越大，PD方法頻散愈加強烈.這表明在利用PD法計算高頻荷載作用時，確定空間步長后，應(yīng)盡量選取較小的影響域，否則會降低其計算精度；圖5(b)給出m= 3時，空間步長Δx對PD方法頻散性質(zhì)的影響，當(dāng)Δx從0.005 m增大到0.02 m時，頻散急劇加重；根據(jù)圖5(c)，當(dāng)保持影響域半徑δ= 0.02 m，m從1增大到4時對應(yīng)的頻散曲線愈接近真實曲線，這表明影響域內(nèi)節(jié)點越密集，PD方法精度越高；圖5(d)給出了Δt從1 μs減小到0.1 μs對PD方法頻散性質(zhì)的影響，發(fā)現(xiàn)頻散曲線對時間步長不太敏感(Δt取值已滿足PD法的穩(wěn)定性條件).通常在模擬裂縫開展時，取影響域半徑為δ= (3~4)Δx[4, 6]，但增大影響域會降低PD方法模擬波傳播問題的精度，這就要求計算時空間步長應(yīng)盡量取小，或者嘗試將預(yù)估裂縫區(qū)處節(jié)點局部進行加密.

3PD算法驗證及應(yīng)用

3.1傳統(tǒng)算法對比

有限差分法(FDM)是一種應(yīng)用廣泛的數(shù)值方法，在巖土力學(xué)問題上具有很好的適用性[15-18].下文將PD方法與FDM進行具體比較，分析兩者之間的差別和聯(lián)系.

不考慮體力和損傷，可將式(8)詳細展開，以嘗試深入分析PD和FDM的差別.一維條件下，影響域δ= Δx時，節(jié)點i的影響域內(nèi)只有兩個節(jié)點(i- 1，i+ 1)，于是式(8)變?yōu)?/p>

(13)

當(dāng)δ=2Δx，3Δx…mΔx時，依次可得：

(14)

(15)

(16)

而一維波動方程經(jīng)過FDM離散后變?yōu)?/p>

(17)

對比發(fā)現(xiàn)，式(17)即式(13)，這說明當(dāng)δ= Δx時，PD法和FDM等價.需要強調(diào)的是，式(13)給出的是一種臨界情況，即節(jié)點i-1和i+1恰好在節(jié)點i影響域邊界上，在PD計算過程中，影響域內(nèi)的節(jié)點位置不斷更新，當(dāng)材料受拉伸時，節(jié)點i-1或i+1可能離開影響域，因此在計算時，δ取值要比節(jié)點間距整數(shù)倍略大來保證計算過程中影響域內(nèi)節(jié)點數(shù)不變，這樣導(dǎo)致FDM和δ=Δx的PD計算結(jié)果會略有差別，下面的算例分析會體現(xiàn)這一點.

根據(jù)泰勒級數(shù)展開可知，式(13)和(14)的截斷誤差分別為O[(Δx)2]和O[(2Δx)2]，這也表明隨著影響域增大，PD方法的誤差增大.一般，當(dāng)影響域半徑δ=mΔx時，得到的PD離散方程可看出一系列差分方程的組合，其截斷誤差變?yōu)镺(δ2)，這樣從理論上說明了影響域變大，該方法頻散嚴重(計算精度降低)的原因.

通過將PD與FDM對比分析，發(fā)現(xiàn)PD離散方程可寫為差分方程組合形式，僅在影響域δ=Δx情況下，PD與FDM計算效果相同，影響域增大之后，其計算精度降低，表明差分法在計算連續(xù)變形方面具有優(yōu)勢.利用PD計算斷裂問題時，通常需要取較大影響域(δ= (3~4)Δx)才能很好模擬材料破壞過程[4, 6]，這反而會降低PD法對波傳播過程的計算精度.為盡可能保證計算精度，在材料破壞之前采用δ=Δx的PD法或直接用FDM，而當(dāng)材料臨近破壞時(文中按連接鍵伸長率達到臨界值考慮)，再采用影響域較大的PD法進行計算，將使得計算結(jié)果更精確.

3.2巖石層裂過程的PD分析

這里以巖石一維層裂試驗為研究對象.根據(jù)一維波動理論，對桿狀試樣入射端施加應(yīng)力脈沖，應(yīng)力波在桿自由端反射形成拉伸波，當(dāng)強度達到巖石的抗拉強度時會引起拉伸斷裂[20].然后根據(jù)層裂位置，就可間接求出層裂強度，并可討論應(yīng)變率效應(yīng)[21-27].根據(jù)PD法頻散性質(zhì)，選擇合適的網(wǎng)格參數(shù)壓制數(shù)值頻散，可有效避免計算求得的應(yīng)力波形在向巖桿自由端傳播過程保持足夠穩(wěn)定，保證計算結(jié)果具有較高精度.

計算中首先施加強度較小的壓力脈沖，采用PD法分析巖桿中波傳播特點，對應(yīng)的定解問題可寫成

(18)

式中p(t)表示桿件所受的沖擊荷載，試驗中荷載時程曲線見圖6，這里將其簡化為上升沿較短的三角波進行分析.三角波峰值pm為8.9 MPa，出現(xiàn)時刻tp=0.04 ms，持續(xù)時間t0=0.12 ms.

首先采用PD法求解該問題，巖桿長l=1 m，取Δx=0.005 m，Δt=1 μs.考慮了4個不同的影響域：δ分別為Δx、2Δx、4Δx和8Δx，同時將FDM(Δx=0.005 m)計算的桿中點(x=0.5 m)應(yīng)力時程曲線給出(見圖7(a)，縱坐標表示軸向應(yīng)力，橫坐標為時間)，規(guī)定拉為正.可看出FDM和δ=Δx的PD方法計算結(jié)果比較接近，但仍存在一定偏差(上文在與FDM對比時已指出)；δ為Δx和2Δx時，圖7(a)得出兩曲線基本吻合；當(dāng)δ增大到4Δx時，t= 0.5 ms后應(yīng)力波波形發(fā)生頻散，計算結(jié)果吻合度變差；當(dāng)δ為8Δx時，從t=0.8 ms后應(yīng)力波因數(shù)值頻散而失真，說明隨著影響域增大，PD法計算結(jié)果變差.圖7(b)中令δ=3Δx，首先保持t0不變，Δx由0.005 m增為0.01 m，約在0.5 ms后出現(xiàn)頻散(帶●標記的曲線)，而帶■標記的曲線無明顯頻散，反映了影響域增大，PD法頻散加重；當(dāng)Δx不變，t0由0.12 ms增大至0.24 ms時(即頻率降低)，波形保持較好，說明隨t0增大，頻散效應(yīng)變?nèi)?通過該對該算例的分析，也可看出PD方法的計算效果與上文頻散分析結(jié)論一致.

采用PD計算斷裂問題時，若Δx，δ等參數(shù)選擇不當(dāng)，產(chǎn)生的數(shù)值頻散會使應(yīng)力波在傳播過程中分散開來，致使尾端加載的壓應(yīng)力波波形偏離原始波形，產(chǎn)生錯誤的計算結(jié)果.

圖6　荷載時程曲線

為說明PD法對層裂模擬分析的效果，方便起見，直接假設(shè)試樣的抗拉強度σt為10 MPa，增大壓力脈沖p(t)，使試樣發(fā)生層裂.計算時考慮m、Δx、pm和t0等參數(shù)的變化，試件斷裂(φ0= 0.3)時尾端損傷分布見圖8(a)所示(縱坐標為損傷，橫坐標表示桿件尾部位置).其中參數(shù)m、Δx、pm、t0、斷裂位置xF，斷裂發(fā)生時間tF列于表1.根據(jù)曲線1~3，可看出當(dāng)影響域半徑由3Δx增大到10Δx時，斷裂的區(qū)域和發(fā)生時間增大，但彼此相差不大.當(dāng)空間步長減小到0.000 5 m時，對比曲線4和5，影響域增大時，斷裂區(qū)域和發(fā)生時間也有變大.同時還發(fā)現(xiàn)，空間步長分別為0.001 和0.000 5 m時，斷裂的區(qū)域和發(fā)生時間相差不大.為了考察沖擊荷載作用時間對層裂位置影響，特將t0增大到0.24 ms并得出斷裂時損傷分布，見圖8(a)中曲線6.對比曲線1和6可看出，層裂位置距自由端變遠，這符合拉應(yīng)力波峰值到時延遲，致使層裂位置遠離自由端.曲線7對應(yīng)的沖擊荷載峰值強度pm為15 MPa，與曲線1相比，層裂位置更加靠近自由端，這是由于沖擊荷載強度提高，反射后拉應(yīng)力峰值達到巖石層裂強度所需時間縮短.圖8(b)給出圖8(a)中曲線1對應(yīng)試件尾部拉應(yīng)力區(qū)隨時間發(fā)展情況，可以發(fā)現(xiàn)隨時間增大，拉應(yīng)力逐漸增大，當(dāng)t=0.439 5 ms時，尾部拉應(yīng)力峰值達到層裂強度，此時巖桿拉伸斷裂瞬間發(fā)生.以上分析表明PD法用于巖石一維層裂模擬可行，效果令人滿意.

圖7　巖桿中點的應(yīng)力時程曲線

圖8　層裂現(xiàn)象PD法模擬

編號mΔx/mpm/MPat0/msΔt/μsxF/mtF/ms130.0010120.120.100.906~0.9090.4395250.0010120.120.100.901~0.9080.43993100.0010120.120.100.900~0.9100.4407430.0005120.120.050.905~0.9070.4389550.0005120.120.050.904~0.9080.4392630.0010120.240.100.811~0.8160.5126730.0010150.120.100.924~0.9280.4327

3.3試驗驗證

試驗中采用的巖桿見圖9(a)，桿長0.99 m，直徑0.07 m，密度2 620 kg·m-3.由于巖石中存在微裂紋，應(yīng)力脈沖在傳播時會發(fā)生衰減.小振幅彈性波在巖石中的衰減是指數(shù)式的，一般情況下可通過彈性波振幅的變化獲得其衰減系數(shù)，關(guān)系式為[26]

(19)

式中:σ0表示初始位置的應(yīng)力峰值，σ(xF)表示傳到距初始位置xF處的應(yīng)力峰值，β為衰減系數(shù).

首先采用低強度應(yīng)力波沖擊巖桿，得到桿中部位置(5#應(yīng)變片，距入射端為0.6 m)應(yīng)變信號，見圖10(縱坐標為電壓，橫坐標為時間).由相鄰拉應(yīng)力波峰值時間差與桿長，求出波速c為2 740 m· s-1，則傳播距離x可由x=ct得出.利用圖10前5個拉應(yīng)力波峰值擬合得出式(19)中的衰減系數(shù)β=0.52 m-1.

加大應(yīng)力波強度，使得巖桿發(fā)生層裂，見圖9(b)，測得其斷裂面距加載端為0.79 m.按照鏡像法，不考慮衰減效應(yīng)，利用5#應(yīng)變片信號，可得到不同時刻自由面附近的應(yīng)力波形，見圖11(縱坐標為軸向應(yīng)力，橫坐標表示桿件尾部位置).然后根據(jù)斷裂位置找出該點的最大拉應(yīng)力，就能確定“名義”層裂強度為σ0= 25.5 MPa.考慮壓應(yīng)力波由5#應(yīng)變片位置傳至自由端形成的拉應(yīng)力波，再到斷裂位置這一過程的衰減(傳播距離x=(0.39+0.20)=0.59 m)，修正后得出層裂強度σt=25.5e-0.52×0.59=18.8 MPa.

圖9　花崗巖桿層裂前后照片

圖10　應(yīng)力波衰減過程

將應(yīng)力脈沖作為邊界條件代入PD計算程序，其中Δx=0.01 m，Δt=2 μs.根據(jù)層裂強度，設(shè)定臨界伸長率s0=σt/E= 0.000 95.考慮3種方法求解：

1)采用δ=3Δx的PD法進行計算；

2)先采用δ=Δx的PD法進行計算，當(dāng)最大連接鍵伸長率smax= 0.9s0時，再令δ=3Δx計算；

3)先采用FDM法進行計算，當(dāng)最大連接鍵伸長率smax= 0.9s0時，再用δ=3Δx的PD法計算.

以上3種方法分別記為：PD3，PD1-PD3和FDM-PD3，得出的斷裂時間和計算耗時見表2，巖桿斷裂時損傷情況見圖12(縱坐標為損傷，橫坐標表示桿件尾部位置).由表2看出，采用FDM-PD3求解耗時最少，PD3耗時最多.PD1-PD3和FDM-PD3得出的斷裂時間接近，而PD3得出的斷裂時間較大，這也反映了影響域增大，頻散加重的特點.由圖12可看出，桿0.79 m處發(fā)生斷裂，與試驗斷裂位置一致，表明該方法能很好模擬桿中拉伸波導(dǎo)致的巖石層裂現(xiàn)象.文章提出的FDM和PD聯(lián)合算法耗時少，精度高，有利于提高PD法對斷裂問題的計算效率.

圖11　應(yīng)力波形疊加

方法斷裂時間/ms計算耗時/sPD30.4782.00PD1-PD30.4660.95FDM-PD30.4640.53

圖12　巖桿斷裂時損傷分布

4結(jié)論

研發(fā)了近場動力學(xué)法計算程序，并詳細分析了該算法頻散特性的影響因素，以及對數(shù)值計算的影響，之后應(yīng)用到一維應(yīng)力波作用下巖桿層裂破壞問題計算分析中，得出如下結(jié)論：

1)當(dāng)空間步長不變時，隨影響域變大，PD法計算精度降低；影響域內(nèi)節(jié)點數(shù)目不變時，空間步長越小，PD法精度會顯著提高；同一影響域內(nèi)節(jié)點分布越密集，精度也有所提高；滿足計算穩(wěn)定性條件的前提下，PD法對時間步長變化不敏感.

2)PD法的離散方程可看作是一系列差分方程的組合，其截斷誤差為O(δ2)，當(dāng)影響域半徑δ等于Δx時，PD算法與中心有限差分法等價，此時計算精度高.當(dāng)影響域增大時，PD法對連續(xù)變形問題計算效果要差于有限差分法.

3)巖桿一維波動層裂模擬結(jié)果與試驗結(jié)果吻合度高.該方法簡單易行，通過選擇合適的計算參數(shù)壓制數(shù)值頻散，能很好地仿真出巖石拉伸波損傷發(fā)展演化直至斷裂破壞的全過程，將FDM和PD法兩者結(jié)合用于脆性材料損傷問題的分析，計算時間少，精度高，能夠發(fā)揮各自優(yōu)勢，提高計算效率.

參考文獻

[1] MADENCI E, OTERKUS E. Peridynamic theory and its applications[M]. New York: Springer-Verlag New York Inc, 2013.

[2] SILLING S A. Reformulation of elasticity theory for discontinuities and long-range forces[J]. Journal of the Mechanics and Physics of Solids, 2000, 48(1):175-209.

[3] SILLING S A, ZIMMERMANN M, ABEYARATNE R. Deformation of a peridynamic bar[J]. Journal of Elasticity, 2003, 73(1/2/3):173-190.

[4] SILLING S A, ASKARI E. A meshfree method based on the peridynamic model of solid mechanics[J]. Computers & Structures, 2005, 83:1526-1535.

[5] HA Y D, BOBARU F. Studies of dynamic crack propagation and crack branching with peridynamics [J]. International Journal of Fracture, 2010, 162(1/2):229-244.

[6] HA Y D, BOBARU F. Characteristics of dynamic brittle fracture captured with peridynamics[J]. Engineering Fracture Mechanics, 2011, 78(6):1156-1168.

[7] WILDMAN R A, GAZONAS G A. A perfectly matched layer for peridynamics in two dimensions[J]. Journal of Mechanics of Materials and Structures 7(8):765-781.

[8] FOSTER J T, SILLING S A, CHEN W W. Viscoplasticity using peridynamics[J]. International Journal for Numerical Methods in Engineering, 2010, 81(10):1242-1258.

[9] HUANG Dan, ZHANG Qing, QIAO Pizhong. Damage and progressive failure of concrete structures using non-local peridynamic modeling[J]. Science China: Technological Sciences, 2011, 54 (3):591-596.

[10]沈峰, 章青, 黃丹，等. 基于近場動力學(xué)理論的混凝土軸拉破壞過程模擬[J].計算力學(xué)學(xué)報, 2013(增1):79-83.

[11] LIU W, HONG J W. Discretized peridynamics for linear elastic solids [J]. Computational Mechanics, 2012, 50(5):579-590.

[12]PARKS M L, PLIMPTON S J, LEHOUCQ R B, et al. Peridynamics with LAMMPS: A user guide, 2008-1035[R]. Albuquerque: Sandia National Laboratories, 2008.

[13]WECKNER O, ABEYARATNE R. The effect of long-range forces on the dynamics of a bar[J]. Journal of the Mechanics and Physics of Solids, 2005, 53(3):705-728.

[14]楊頂輝, 滕吉文. 各向異性介質(zhì)中三分量地震記錄的FCT有限差分模擬[J]. 石油地球物理勘探, 1997, (2):181-190.

[15]吳國忱, 王華忠. 波場模擬中的數(shù)值頻散分析與校正策略[J]. 地球物理學(xué)進展, 2005, 20(1):58-65.

[16]劉東甲. 完整樁瞬態(tài)縱向振動的模擬計算[J]. 合肥工業(yè)大學(xué)學(xué)報(自然科學(xué)版)，2000，23(5): 683-687.

[17]柯宅邦, 劉東甲. 低應(yīng)變反射波法測樁的軸對稱問題數(shù)值計算[J]. 巖土工程學(xué)報, 2006, 28(12):2111-2115.

[18]盧志堂, 劉東甲, 龍麗麗, 等. 基樁低應(yīng)變檢測三維問題的數(shù)值計算[J]. 合肥工業(yè)大學(xué)學(xué)報(自然科學(xué)版)，2011，34(6):905-909.

[19]盧志堂, 王志亮, 劉東甲，等. 管樁低應(yīng)變檢測中的三維效應(yīng)分析[J]. 同濟大學(xué)學(xué)報(自然科學(xué)版), 2012, 40(11):1603-1607.

[20] 王禮立. 應(yīng)力波基礎(chǔ)[M]. 北京: 國防工業(yè)出版社, 2005:1-152.

[21]KLEPACZHO J R, BRARA A. An experiment method for dynamic tensile testing of concrete by spalling[J]. International Journal of Impact Engineering, 2001，25 (4):387-409.

[22]胡時勝, 張磊, 武海軍, 等. 混凝土材料層裂強度的實驗研究[J].工程力學(xué), 2004, 21 (4):128-132.

[23]BRARA A, KLEPACZKO J R. Experimental characterization of concrete in dynamic tension[J]. Mechanics of Materials, 2006, 38(3):253-267.

[24]李夕兵, 陶明, 宮鳳強,等. 沖擊載荷作用下硬巖層裂破壞的理論和試驗研究[J]. 巖石力學(xué)與工程學(xué)報, 2011, 30(6):1081-1088.

[25] WANG Zhiliang, ZHOU Nianqing, WANG Jianguo. Using Hopkinson pressure bar to perform dynamic tensile tests on SFRC at medium strain rates[J]. Magazine of Concrete Research, 2012, 64 (8):657-664.

[26]王志亮, 李洋, 陽棟. C75混凝土動態(tài)層裂強度試驗研究[J]. 建筑材料學(xué)報, 2014, 17 (4):695-699.

[27]WANG ZhiLiang, YANG Dong. Study on dynamic behavior and tensile strength of concrete using 1D wave propagation characteristics[J]. Mechanics of Advanced Materials and Structures, 2015, 22(3):829-836.

(編輯趙麗瑩)

Dispersion characteristics of peridynamics method and its application to spalling analysis of rock

LU Zhitang1, WANG Zhiliang1,2

(1.Department of Geotechnical Engineering，Tongji University，200092 Shanghai，China；2.School of Civil & Hydraulic Engineering, Hefei University of Technology, 230009 Hefei, China)

Abstract:The dispersion characteristics of PD is carried out to figure out the accuracy of peridynamics method (PD) and its simulation effects on rock spalling. Firstly, based on the dispersion analyses of the PD discrete equations, it is found that the numerical dispersion is more apparent with increasing the domain size if the space step is fixed; with the number of nodes in the domain fixed, the numerical dispersion becomes weaker when the space step decreases; the numerical dispersion also gets weaker with the increase of the nodes at a fixed domain. Subsequently, by comparing it to the traditional finite difference method, it is pointed out that the PD discrete equations can be written as the combination of a series of differential equations, whose truncation error is the second order infinitesimal of the domain radius (δ); when δ is set as the space step, PD method is equivalent to the central difference method, owning the highest calculational accuracy. Finally, the PD method is used to explore the one-dimensional spalling phenomenon of rock bar. The effects of space step and domain size on computing results are discussed, and the spalling time, spalling locations and damage distribution are further given. The effectiveness of PD method is also verified by contrast with the spalling test. It is shown that PD can be used for rock spalling analysis and high-precision results can be obtained by PD method coupled with FDM, costing less time and owning a clear advantage.

Keywords:peridynamics；dispersion characteristics；finite differential method；rock bar；spalling

中圖分類號:TU45

文獻標志碼:A

文章編號:0367-6234(2016)02-0131-07

通信作者:王志亮,cvewzL@#edu.cn.

作者簡介:盧志堂(1985—)，男，博士研究生;王志亮(1969—)，男，教授，博士生導(dǎo)師.

基金項目:國家自然科學(xué)基金 (51174145，51379147，51579062)；

收稿日期:2015-06-23.

doi：10.11918/j.issn.0367-6234.2016.02.022

教育部博士點專項資金(20120072110024).