摘 要：流形元是新出現(xiàn)的一種被較為廣泛地應(yīng)用于材料破壞模擬的數(shù)值分析方法。為了驗(yàn)證該方法在材料破壞模擬中的有效性，分別利用流形元與有限元兩種不同的數(shù)值方法對(duì)巖石沖擊破壞過程進(jìn)行了模擬分析。模擬結(jié)果表明，流形元法對(duì)材料破壞的模擬結(jié)果更符合實(shí)際情況，這主要是由于該方法采用了兩種不同的覆蓋系統(tǒng)——數(shù)學(xué)覆蓋和物理覆蓋，并引入能夠客觀反映材料裂紋的產(chǎn)生與擴(kuò)展準(zhǔn)則，以及相應(yīng)的塊體運(yùn)動(dòng)理論。流形元的出現(xiàn)有望對(duì)材料破壞模擬開創(chuàng)出一條新的途徑。

關(guān)鍵詞：巖石沖擊破壞；數(shù)值模擬；流形元；有限元；對(duì)比分析

自1960年Clough最先在進(jìn)行飛機(jī)結(jié)構(gòu)分析時(shí)提出有限元以來^[1]，數(shù)值計(jì)算方法就開始在工程結(jié)構(gòu)分析、流體分析及場分析等眾多領(lǐng)域中得到廣泛的應(yīng)用，同時(shí)數(shù)值方法本身也得到不斷的發(fā)展與完善。目前從種類上來說，主要包括有限元法、離散元法和流形元三類。這三類數(shù)值方法在實(shí)際應(yīng)用中各有一定的優(yōu)缺點(diǎn)，如有限元法在求解連續(xù)介質(zhì)發(fā)生連續(xù)變形時(shí)非常有效，而離散元法在求解非連續(xù)材料發(fā)生非連續(xù)變形時(shí)非常有效。然而隨著生產(chǎn)實(shí)踐的發(fā)展，又提出了一些更新的問題如初始完整性很好的材料在外力作用下產(chǎn)生破壞的過程模擬等，由于有限元法是以連續(xù)性為基礎(chǔ)，它在模擬材料破壞方面有著不可克服的困難，于是就采用一種近似的方法來模擬材料的破壞，即以某一物理量作為材料的破壞準(zhǔn)則，如應(yīng)力、應(yīng)變或損傷等，然后確定一個(gè)閾值，當(dāng)某單元的該物理量的值超過閾值時(shí)就認(rèn)為該單元已經(jīng)破壞，即把該單元從整個(gè)結(jié)構(gòu)中刪除，實(shí)際上是把單元的質(zhì)量、剛度等乘以一個(gè)很小的數(shù)值，使該單元不參加以后的計(jì)算。這種方法在目前的計(jì)算中也得到了一定的應(yīng)用^[2-4]，但其顯然是一種近似的方法，也可以說是有限元為彌補(bǔ)其在材料破壞模擬方面的不足而采用的一個(gè)補(bǔ)救措施，與實(shí)際情況存在著一定的差距。同樣盡管離散元法在模擬材料的破壞方面有著巨大的優(yōu)勢，但是它卻不能模擬新裂紋的產(chǎn)生，盡管目前也有不少學(xué)者采用虛擬裂紋的方法來模擬完整材料或部分完整材料破壞過程，并取得了良好的效果^[5]，但顯然這也是一種近似的方法。

于是流形元即數(shù)值流形方法(NMM)就應(yīng)運(yùn)而生，該方法是由美籍華人石根華博士繼提出不連續(xù)變形分析方法(DDA)之后，又提出的一種新的數(shù)值方法^[6]。雖然該方法是在20世紀(jì)90年代初才被提出，但由于其在模擬材料破壞方面具有明顯的優(yōu)越性，一開始就吸引了眾多學(xué)者的關(guān)注。到目前為止，該方法已在材料破壞模擬分析中得到了成功應(yīng)用^[7,8]；同時(shí)該方法在塊體運(yùn)動(dòng)模擬方面，完全吸收了DDA中關(guān)于塊體運(yùn)動(dòng)的理論，能夠很好地模擬塊體破壞后的飛散過程，這一重大進(jìn)展對(duì)以連續(xù)介質(zhì)為基礎(chǔ)的有限元法來說是一個(gè)重大的突破，克服了在利用有限元法計(jì)算時(shí)僅能給出待分析域內(nèi)應(yīng)力分布，而不能模擬其破碎及破碎后塊體運(yùn)動(dòng)等情況的不足。目前關(guān)于流形元模擬材料破壞、運(yùn)動(dòng)等方面已有不少的相關(guān)研究成果發(fā)表。但是至今還沒有看到關(guān)于有限元與流形元的相關(guān)對(duì)比研究成果發(fā)表，如果能夠采用這兩種方法對(duì)同一問題進(jìn)行分析研究，則更能夠很好地反映出它們之間的相互關(guān)系及各自的優(yōu)缺點(diǎn)，因此本文就采用這兩種方法對(duì)巖石沖擊破壞過程進(jìn)行模擬研究，并對(duì)其模擬結(jié)果進(jìn)行了對(duì)比分析。

 

巖石在沖擊載荷作用下的破壞問題屬于明顯的斷裂動(dòng)力學(xué)問題。斷裂動(dòng)力學(xué)就是研究那些慣性效應(yīng)不能忽略的斷裂力學(xué)問題^[9]，其求解方法明顯不同于斷裂靜力學(xué)問題，它們的一個(gè)最重要的區(qū)別就是材料的慣性效應(yīng)不能忽略，考慮了慣性效應(yīng)的斷裂力學(xué)就是斷裂動(dòng)力學(xué)或動(dòng)態(tài)斷裂力學(xué)。在動(dòng)態(tài)加載時(shí)，試件除產(chǎn)生彈塑性變形外，內(nèi)部各質(zhì)點(diǎn)的自由振動(dòng)要獲得一定的加速度，從而產(chǎn)生慣性力，這就是所謂的動(dòng)態(tài)加載時(shí)的慣性效應(yīng)。

在動(dòng)力問題的分析中，數(shù)值流形方法引入了慣性矩陣，它相當(dāng)于有限元方法中的質(zhì)量矩陣，以充分考慮動(dòng)力學(xué)問題中的慣性效應(yīng)。數(shù)值流形方法在處理動(dòng)力學(xué)問題時(shí)，與處理靜力學(xué)問題相比的一個(gè)最大差別就是在當(dāng)前步的計(jì)算中，各單元繼承了前一時(shí)間步的速度，而不是像處理靜力學(xué)問題一樣置當(dāng)前時(shí)間步的速度為零。而這種處理方法對(duì)一般動(dòng)力學(xué)問題來說過于簡單，并且這種處理方法也僅僅考慮了單元的慣性效應(yīng)，而沒有考慮單元的阻尼，所以還不夠精確。因此這里采用的數(shù)值流形方法計(jì)算程序是借用動(dòng)力有限元法的求解思想，利用動(dòng)力有限元中的Newmark方法改進(jìn)后的數(shù)值流形方法計(jì)算程序，改進(jìn)的原理及改進(jìn)后的程序優(yōu)越性已在相關(guān)文獻(xiàn)中進(jìn)行了論述^[7]。有限元的動(dòng)力計(jì)算方法同樣是采用瞬態(tài)動(dòng)力分析Newmark方法。

 

由于目前數(shù)值流形方法在實(shí)際中的應(yīng)用還不是太廣泛，所以有必要首先對(duì)其計(jì)算原理進(jìn)行簡單介紹。

流形方法在求解動(dòng)力學(xué)問題中所使用的方程為結(jié)構(gòu)動(dòng)力學(xué)方程：

 式中：M為質(zhì)量矩陣；C為阻尼矩陣；△δ為位移增量；分別是位移速度和加速度；K=K_e+K_cn+K_cs+K_f，其中K_e是剛度矩陣，K_cn、K_cs分別為塊體及不連續(xù)面之間的接觸矩陣，K_f是約束矩陣。△F為總載荷增量矩陣，△F=F_p+F_b+F_f-F₀+F_cn+F_cs+F_fr，其中F_p是外載荷向量，F(xiàn)_b是體積力向量，F(xiàn)_f是已知約束位移引起的等效載荷向量，F(xiàn)₀是初應(yīng)力向量；F_cn、F_cs分別為法向和切向接觸引起的等效載荷向量；F_fr，為接觸面之間的摩擦力引起的等效載荷向量。

在外力作用下，材料的破壞過程一方面是由于材料內(nèi)部已存在的裂紋被激活擴(kuò)展，而另一方面就是在外載作用下，材料內(nèi)部的應(yīng)力超過了材料的抗拉或抗剪強(qiáng)度而導(dǎo)致材料被拉斷或剪壞。而對(duì)于巖石等抗壓強(qiáng)度遠(yuǎn)遠(yuǎn)大于其抗拉和抗剪強(qiáng)度的脆性材料來說，其破壞的主要形式為拉伸破壞和剪切破壞。

材料的破壞過程是裂紋在外力作用下不斷產(chǎn)生、聚集和匯合的過程。在流形程序判斷新裂紋的起裂和已存在裂紋的擴(kuò)展問題中，采用了不同的標(biāo)準(zhǔn)。對(duì)于已存在的裂紋是否擴(kuò)展的問題采用斷裂力學(xué)的應(yīng)力強(qiáng)度因子作為判斷標(biāo)準(zhǔn)，根據(jù)線彈性斷裂力學(xué)來判斷裂紋的起裂與擴(kuò)展^[7]。對(duì)于新裂紋的產(chǎn)生，一般考慮采用以應(yīng)力為判據(jù)。這里采用摩爾-庫侖準(zhǔn)則作為裂紋起裂的判據(jù)，其判斷方法如下：①當(dāng)?shù)?主應(yīng)力大于材料的抗拉強(qiáng)度時(shí)，材料被拉伸破壞，產(chǎn)生新裂紋；②當(dāng)某點(diǎn)的最大剪應(yīng)力大于材料的抗剪強(qiáng)度時(shí)，材料被剪切破壞，產(chǎn)生新裂紋。

在外力作用下，巖體內(nèi)的初始裂紋將發(fā)生擴(kuò)展并最終把巖體切割成相應(yīng)的塊體，這些塊體在自重或其他外力作用下產(chǎn)生運(yùn)動(dòng)。數(shù)值流形方法對(duì)塊體的運(yùn)動(dòng)模擬中，一定要保證塊體之間無嵌入、無拉伸，它完全借鑒了DDA中塊體運(yùn)動(dòng)學(xué)的理論。而這一過程是通過反復(fù)的開一合迭代，即在可能的接觸邊界之間施加和移去彈簧來實(shí)現(xiàn)的。在邊界處塊體之間有張開、鎖定和滑動(dòng)三種狀態(tài)，當(dāng)邊界從張開發(fā)展到鎖定時(shí)施加切向和法向彈簧；當(dāng)邊界從滑動(dòng)到鎖定時(shí)加切向彈簧；當(dāng)邊界之間滑動(dòng)時(shí)在邊界兩側(cè)加摩擦力，從而對(duì)剛度矩陣作出修正，計(jì)算接觸力并疊加于本次計(jì)算的平衡力列陣中，參加下一次的迭代計(jì)算。

 

在本算例中，取爆炸載荷作為施加的沖擊載荷。相應(yīng)的實(shí)際問題是對(duì)應(yīng)于二維平面問題的圓形裝藥在有限域內(nèi)爆炸后形成爆破漏斗的過程模擬。炸藥在炮孔中起爆后，巖石將發(fā)生如下的破碎過程：

(1)強(qiáng)大的沖擊波壓應(yīng)力使炮孔周圍巖石受壓破碎，在瞬時(shí)形成壓縮破碎和初始裂隙。

(2)環(huán)向拉應(yīng)力及應(yīng)力波反射拉應(yīng)力使巖石中的裂隙擴(kuò)展，引起巖石進(jìn)一步破裂，包括初始裂隙的擴(kuò)展和二次裂隙的形成。

(3)爆生氣體膨脹作用使巖石中的裂隙貫穿形成碎塊，巖塊運(yùn)動(dòng)形成爆破漏斗。

由于實(shí)際的爆破過程十分復(fù)雜，因此在滿足工程要求的條件下，可以對(duì)巖石爆破破壞的過程進(jìn)行簡化處理，通過在炮孔內(nèi)壁上施加均布的沖擊三角波載荷來模擬爆炸荷載，沖擊載荷的升壓時(shí)間為80ms，降壓時(shí)間為220ms，整個(gè)荷載作用時(shí)間為300ms，載荷的最大峰值壓力為5GPa。

所建立的計(jì)算模型是二維分析域內(nèi)的圓形裝藥在有限域巖體內(nèi)爆炸后，形成爆破漏斗的過程。取計(jì)算模型的尺寸為4m×1.5m，圓形炮孔直徑為O.5m，炮孔中心距上邊界的距離為0.5m，距左右邊界的距離均為2m。假定巖石為各向同性均質(zhì)巖體，其中模型的上邊界為自由邊界，其余三個(gè)邊界均為固定邊界，計(jì)算模型示意圖如圖1所示。

在本算例中，所分析問題的類型為動(dòng)態(tài)類型，因此所采用的材料計(jì)算參數(shù)均應(yīng)為動(dòng)態(tài)參數(shù)。取大理巖的參數(shù)作為本算例中的巖石參數(shù)：動(dòng)彈性模量80GPa，動(dòng)泊松比0.20，密度2400kg／m³，動(dòng)態(tài)斷裂韌性0.5MN／m^3/2，節(jié)理面的摩擦角45°，黏結(jié)力10MPa，動(dòng)抗拉強(qiáng)度30MPa。

4.2.1流形元計(jì)算結(jié)果

取計(jì)算過程中的3個(gè)瞬時(shí)狀態(tài)如圖2(a)～(c)所示。

4.2.2有限元計(jì)算結(jié)果

采用目前應(yīng)用較為廣泛的有限元數(shù)值計(jì)算軟件ANSYS對(duì)上述模型進(jìn)行計(jì)算分析，模擬結(jié)果如圖3（略）所示。由于有限元軟件還不能很好地考慮裂紋的產(chǎn)生及塊體的運(yùn)動(dòng)，所以采用了有限元法中對(duì)材料破壞模擬的近似方法——單元生死法，即采用某一準(zhǔn)則對(duì)計(jì)算結(jié)果中的單元進(jìn)行刪除，即刪除失效的單元。這里采用Mises有效應(yīng)力失效準(zhǔn)則，當(dāng)某單元的Mises應(yīng)力超過1MPa時(shí)，就認(rèn)為該單元已經(jīng)失效，并把該單元從計(jì)算模型中刪除，用于表示巖石的破壞；同時(shí)還給出了模型中的應(yīng)力分布情況。

4.2.3兩種計(jì)算結(jié)果的對(duì)比分析

由上述流形元和有限元模擬的結(jié)果可以看出：

(1)由于流形元采用的是材料的拉、剪破壞準(zhǔn)則和線彈性斷裂力學(xué)理論，因此它可以很容易地模擬完整巖石在沖擊荷載作用下新裂紋的起裂及擴(kuò)展；對(duì)于塊體的運(yùn)動(dòng)采用的足DDA中的塊體運(yùn)動(dòng)理論，所以它也能夠很好地模擬塊體的運(yùn)動(dòng)過程。而相比之下，由于有限元法是以連續(xù)介質(zhì)為基礎(chǔ)，它不具備考慮巖石中新裂紋的產(chǎn)生及擴(kuò)展的功能，所以也很難得出符合巖石實(shí)際破壞過程的計(jì)算結(jié)果。但有限元為了克服其在模擬材料破壞方面的不足，引入了單元?jiǎng)h除的方法，而圖3(B)所給出的結(jié)果就是采用單元?jiǎng)h除的方法得到的巖石破壞情況計(jì)算結(jié)果，這在一定程度上彌補(bǔ)了有限元在這方面的不足。

(2)從對(duì)巖石破壞過程模擬的準(zhǔn)確性來看，巖石的破壞是一種漸進(jìn)的過程，也就是說巖石的破壞是從無裂紋到一條裂紋再到多條裂紋，以至于最后的完全破壞。圖2的流形元模擬結(jié)果非�？陀^地反映了這種破壞過程，即材料首先產(chǎn)生一條裂紋，然后再到多條裂紋，最后切割成塊體。相比之下，有限元的計(jì)算結(jié)果則無法反映這種客觀的破壞過程，在O.7ms以前巖石還處于少數(shù)幾條裂紋產(chǎn)生與擴(kuò)展階段時(shí)，有限元的計(jì)算結(jié)果則幾乎無法反映這種少量、局部的破壞情況；而當(dāng)破壞的范圍逐漸增大，即在3ms以后有限元?jiǎng)t以應(yīng)力的形式體現(xiàn)出了巖石中的應(yīng)力不同，進(jìn)而可以根據(jù)應(yīng)力失效準(zhǔn)則得出材料可能的破壞區(qū)域，但有限元仍無法模擬出巖石破壞后的塊體運(yùn)動(dòng)情況。因此，從對(duì)巖石破壞過程模擬的準(zhǔn)確性來看，流形元法則顯得更為優(yōu)越。

(3)從對(duì)爆炸空腔的模擬結(jié)果來看，流形元法得到的爆炸空腔是一個(gè)倒三角形的形狀，這與巖石在近地表處的炸藥爆炸破壞的實(shí)際結(jié)果是符合的，這是因?yàn)樵诮�地表處由于自由面的影響，裂紋很容易擴(kuò)展達(dá)到地表，形成爆破空腔。有限元模擬的爆破空腔則是在炮孔周圍較大，而在地表則沒有形成一個(gè)倒三角形的爆破漏斗。造成兩者差別的根本原因是兩種數(shù)值方法采用的算法不同，由于流形元采用的是斷裂力學(xué)方法，當(dāng)巖石中出現(xiàn)初始裂紋以后，由外力作用產(chǎn)生的能量將很容易使已有裂紋繼續(xù)擴(kuò)展，這樣就很容易形成圖2(c)中的兩條優(yōu)勢裂紋。當(dāng)這兩條優(yōu)勢裂紋達(dá)到地表以后，就形成了所謂的爆破漏斗，即使這兩條優(yōu)勢裂紋中間的巖石塊體沒有得到破碎也絲毫不影響漏斗區(qū)的形成。而有限元方法則不同，由于其沒有裂紋的形成，所以也不存在能量的優(yōu)先分配，所有的能量幾乎都是按照單元距離或位置對(duì)稱的原則進(jìn)行分配，所以其形成的爆破空腔也是非常對(duì)稱的，顯然這種數(shù)學(xué)上的對(duì)稱是不符合巖石破壞的實(shí)際情況的。

(4)從裂紋的發(fā)展情況來看，流形元在模擬巖石的沖擊破壞過程中裂紋的發(fā)展，存在著明顯的分岔現(xiàn)象。即裂紋并不是沿著最初出現(xiàn)的裂紋一直向前發(fā)展，而是在發(fā)展過程中會(huì)偏離原來的路徑而分為兩條或多條分支裂紋，甚至在原裂紋沒有到達(dá)的地方出現(xiàn)了新的小裂紋，這是裂紋動(dòng)態(tài)擴(kuò)展中一個(gè)最常見的現(xiàn)象之一^[10]。而有限元的模擬中卻沒有發(fā)現(xiàn)這一點(diǎn)，這也是由于有限元的局限性所致。

 

本文對(duì)目前應(yīng)用最普遍的有限元和最新出現(xiàn)的流形元兩種數(shù)值方法，在計(jì)算巖石沖擊破壞中的優(yōu)缺點(diǎn)進(jìn)行了對(duì)比分析，說明了流形元在材料破壞模擬中的優(yōu)勢，即由于流形元是根據(jù)材料的實(shí)際破壞特點(diǎn)采用相應(yīng)的裂紋起裂及擴(kuò)展準(zhǔn)則進(jìn)行計(jì)算的，因此對(duì)材料破壞過程的模擬與實(shí)際情況是比較吻合的。而有限元法則是以連續(xù)介質(zhì)為基礎(chǔ)，沒有考慮到材料的實(shí)際破壞過程，而是以應(yīng)力或應(yīng)變等失效準(zhǔn)則作為判斷材料破壞的判據(jù)，因而計(jì)算出的材料破壞情況與實(shí)際有較大的差距。因此，以連續(xù)介質(zhì)為基礎(chǔ)發(fā)展起來的有限元法在計(jì)算連續(xù)材料在外力作用下產(chǎn)生連續(xù)變形方面是具有很大優(yōu)勢的，而在計(jì)算不連續(xù)破壞方面則具有一定的局限性。由于流形元法是以物理覆蓋和數(shù)學(xué)覆蓋雙重覆蓋為基礎(chǔ)建立的，因而能夠方便地考慮材料在外力作用下的破壞過程，破壞后單元的重新構(gòu)建也十分方便，同時(shí)由于其在塊體運(yùn)動(dòng)方面完全借鑒了DDA中的塊體運(yùn)動(dòng)理論，因此對(duì)材料破壞后的運(yùn)動(dòng)模擬也是十分有效，有望對(duì)材料的破壞模擬研究開創(chuàng)出一條新途徑。

摘自《工程爆破》總第61期

[1]Clough R W．The Finite Element in Plane Stress Analysis [A]．In：Proc．2^ndASCE Conf．on Electronic Computation [C]．1960．

[2]Zhi-liang Wang，Yong-chi Li，R．F．Shen．Numerical Simulation of Tensile Damage and Blast Crater in Brittle Rock Due to underground Explosion[J]．International

Journal of Rock Mechanics & Mining Sciences。2007，44：730—738．

[3]尹潘．巖石破壞機(jī)理的LCA模型數(shù)值模擬[J]．遼寧工程技術(shù)大學(xué)學(xué)報(bào)，2006，25(3)：364—366．

[4]王來貴，趙娜，周永發(fā)，等．巖石受拉破壞的數(shù)值模擬方法[J]．返寧工程技術(shù)大學(xué)學(xué)報(bào)，2007，26(2)：198—200．

[5]張國新，武曉峰．裂隙滲流對(duì)巖石邊坡穩(wěn)定的影響——滲流、變形耦合作用的DDA法[J]．巖石力學(xué)與工程學(xué)報(bào)，2003，22(8)：1269—1275．

[6]Shi Genhua．Manifold Method of Material Analysis[A]．In：Transactions of the Ninth Army Conference on Applied Mathematics and Computing[C]．Minneapolish，Minncsoda，USA：1992，51—76．

[7]劉紅巖，秦四清，楊軍．爆炸荷載下巖石破壞的數(shù)值流形方法模擬[J]．爆炸與沖擊，2007，27(1)：50一56．

[8]張國新，趙妍，石根華，等．模擬巖石邊坡傾倒破壞的數(shù)值流形方法[J]．巖土工程學(xué)報(bào)，2007，29(6)：800—805．

[9]范天佑．?dāng)嗔牙碚摶�礎(chǔ)[M]．北京：科學(xué)出版社，2003．

[10]劉再華，解德，王元漢，等．工程斷裂動(dòng)力學(xué)[M]．武漢：華中理工大學(xué)出版社，1996．