盧志堂， 王志亮, 2

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

?

近場動力學法頻散特性及其在巖石層裂分析中應用

盧志堂1， 王志亮1, 2

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

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

關鍵詞:近場動力學；頻散特性；有限差分法；巖桿；層裂

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

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

1PD法計算原理

1.1基本方程

如圖1所示，材料占據空間域B，x為B內一材料點，Hx為x的影響域，δ為Hx的半徑.假設x與其影響域Hx內的另一材料點x′通過連接鍵產生相互作用，作用力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　影響域內相對位置與相對位移

1.2本構關系與損傷準則

f(η,ξ)可寫成連接鍵伸長率的函數，即PD本構關系式：

(4)

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

(5)

(6)

式中s0為連接鍵臨界伸長率，當s ≥ s0時，連接鍵斷裂，且不可恢復.

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

(7)

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

1.3求解思路

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

(8)

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

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

其中V = (Δx)3.

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

(10)

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

PD法數值計算流程見圖4.

圖3　PD計算的離散域示意

圖4　PD程序流程圖

2頻散特性分析

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

(11)

(12)

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

圖5　頻散特性影響因素

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

3PD算法驗證及應用

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

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

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

(13)

當δ=2Δx，3Δx…mΔx時，依次可得：

(14)

(15)

(16)

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

(17)

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

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

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

3.2巖石層裂過程的PD分析

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

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

(18)

式中p(t)表示桿件所受的沖擊荷載，試驗中荷載時程曲線見圖6，這里將其簡化為上升沿較短的三角波進行分析.三角波峰值pm為8.9 MPa，出現時刻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)應力時程曲線給出(見圖7(a)，縱坐標表示軸向應力，橫坐標為時間)，規(guī)定拉為正.可看出FDM和δ=Δx的PD方法計算結果比較接近，但仍存在一定偏差(上文在與FDM對比時已指出)；δ為Δx和2Δx時，圖7(a)得出兩曲線基本吻合；當δ增大到4Δx時，t= 0.5 ms后應力波波形發(fā)生頻散，計算結果吻合度變差；當δ為8Δx時，從t=0.8 ms后應力波因數值頻散而失真，說明隨著影響域增大，PD法計算結果變差.圖7(b)中令δ=3Δx，首先保持t0不變，Δx由0.005 m增為0.01 m，約在0.5 ms后出現頻散(帶●標記的曲線)，而帶■標記的曲線無明顯頻散，反映了影響域增大，PD法頻散加重；當Δx不變，t0由0.12 ms增大至0.24 ms時(即頻率降低)，波形保持較好，說明隨t0增大，頻散效應變弱.通過該對該算例的分析，也可看出PD方法的計算效果與上文頻散分析結論一致.

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

圖6　荷載時程曲線

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

圖7　巖桿中點的應力時程曲線

圖8　層裂現象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.由于巖石中存在微裂紋，應力脈沖在傳播時會發(fā)生衰減.小振幅彈性波在巖石中的衰減是指數式的，一般情況下可通過彈性波振幅的變化獲得其衰減系數，關系式為[26]

(19)

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

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

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

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

圖10　應力波衰減過程

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

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

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

3)先采用FDM法進行計算，當最大連接鍵伸長率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ā)生斷裂，與試驗斷裂位置一致，表明該方法能很好模擬桿中拉伸波導致的巖石層裂現象.文章提出的FDM和PD聯合算法耗時少，精度高，有利于提高PD法對斷裂問題的計算效率.

圖11　應力波形疊加

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

圖12　巖桿斷裂時損傷分布

4結論

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

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

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

3)巖桿一維波動層裂模擬結果與試驗結果吻合度高.該方法簡單易行，通過選擇合適的計算參數壓制數值頻散，能很好地仿真出巖石拉伸波損傷發(fā)展演化直至斷裂破壞的全過程，將FDM和PD法兩者結合用于脆性材料損傷問題的分析，計算時間少，精度高，能夠發(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]沈峰, 章青, 黃丹，等. 基于近場動力學理論的混凝土軸拉破壞過程模擬[J].計算力學學報, 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]楊頂輝, 滕吉文. 各向異性介質中三分量地震記錄的FCT有限差分模擬[J]. 石油地球物理勘探, 1997, (2):181-190.

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

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

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

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

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

[20] 王禮立. 應力波基礎[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].工程力學, 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]. 巖石力學與工程學報, 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]. 建筑材料學報, 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—)，男，教授，博士生導師.

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

收稿日期:2015-06-23.

doi：10.11918/j.issn.0367-6234.2016.02.022

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