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

盧志堂， 王志亮, 2

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

?

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

盧志堂1， 王志亮1, 2

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

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

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

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

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

1PD法計(jì)算原理

1.1基本方程

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

(1)

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

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

材料點(diǎn)x的PD方程為

(2)

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

相互作用力f可寫為

(3)

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

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

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

(4)

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

(5)

(6)

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

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

(7)

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

1.3求解思路

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

(8)

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

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

其中V = (Δx)3.

采用時(shí)域中心差分計(jì)算加速度：

(10)

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

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

圖3　PD計(jì)算的離散域示意

圖4　PD程序流程圖

2頻散特性分析

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

(11)

(12)

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

圖5　頻散特性影響因素

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

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

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

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

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

(13)

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

(14)

(15)

(16)

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

(17)

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

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

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

3.2巖石層裂過程的PD分析

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

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

(18)

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

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

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

圖6　荷載時(shí)程曲線

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

圖7　巖桿中點(diǎn)的應(yīng)力時(shí)程曲線

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

編號(hào)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試驗(yàn)驗(yàn)證

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

(19)

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

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

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

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

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

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

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

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

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

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

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

方法斷裂時(shí)間/ms計(jì)算耗時(shí)/sPD30.4782.00PD1-PD30.4660.95FDM-PD30.4640.53

圖12　巖桿斷裂時(shí)損傷分布

4結(jié)論

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

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

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

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

參考文獻(xiàn)

[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]沈峰, 章青, 黃丹，等. 基于近場(chǎng)動(dòng)力學(xué)理論的混凝土軸拉破壞過程模擬[J].計(jì)算力學(xué)學(xué)報(bào), 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]吳國(guó)忱, 王華忠. 波場(chǎng)模擬中的數(shù)值頻散分析與校正策略[J]. 地球物理學(xué)進(jìn)展, 2005, 20(1):58-65.

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

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

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

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

[20] 王禮立. 應(yīng)力波基礎(chǔ)[M]. 北京: 國(guó)防工業(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]胡時(shí)勝, 張磊, 武海軍, 等. 混凝土材料層裂強(qiáng)度的實(shí)驗(yàn)研究[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]李夕兵, 陶明, 宮鳳強(qiáng),等. 沖擊載荷作用下硬巖層裂破壞的理論和試驗(yàn)研究[J]. 巖石力學(xué)與工程學(xué)報(bào), 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混凝土動(dòng)態(tài)層裂強(qiáng)度試驗(yàn)研究[J]. 建筑材料學(xué)報(bào), 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

中圖分類號(hào):TU45

文獻(xiàn)標(biāo)志碼:A

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

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

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

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

收稿日期:2015-06-23.

doi：10.11918/j.issn.0367-6234.2016.02.022

教育部博士點(diǎn)專項(xiàng)資金(20120072110024).