魏 炯，高富強(qiáng)

(1.中煤科工開(kāi)采研究院有限公司，北京 100013；2.煤炭科學(xué)研究總院 煤炭資源高效開(kāi)采與潔凈利用國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100013)

巖石力學(xué)是一門(mén)研究巖石材料力學(xué)性能的理論和應(yīng)用學(xué)科，尤其注重于研究巖石的內(nèi)力、變形、空間和時(shí)間之間的相互關(guān)系。隨著計(jì)算機(jī)和計(jì)算數(shù)學(xué)的發(fā)展，科學(xué)計(jì)算繼理論研究和科學(xué)實(shí)驗(yàn)之后，成為現(xiàn)代科學(xué)研究的第三種方法[1]。在過(guò)去數(shù)十年中，數(shù)值計(jì)算巖石力學(xué)蓬勃發(fā)展，采用數(shù)值計(jì)算研究巖石在載荷作用下的變形與破裂等力學(xué)行為成為巖石力學(xué)研究的熱點(diǎn)[2-3]。

數(shù)值計(jì)算可以理解為用計(jì)算機(jī)來(lái)做實(shí)驗(yàn)，數(shù)值計(jì)算結(jié)果應(yīng)與室內(nèi)物理實(shí)驗(yàn)結(jié)果存在一一對(duì)應(yīng)的關(guān)系。但由于計(jì)算能力的限制，一般將三維問(wèn)題簡(jiǎn)化二維問(wèn)題，將準(zhǔn)靜態(tài)問(wèn)題簡(jiǎn)化為靜力問(wèn)題。事實(shí)上，自然界所有的物質(zhì)都處于由空間和時(shí)間構(gòu)成的四維時(shí)空之中，任何力學(xué)過(guò)程都是與空間和時(shí)間相關(guān)的過(guò)程。人們已經(jīng)認(rèn)識(shí)到巖石力學(xué)問(wèn)題是三維的力學(xué)問(wèn)題，將計(jì)算巖石力學(xué)從二維擴(kuò)展到了三維，然而準(zhǔn)靜態(tài)加載條件下的時(shí)間維度依然為人們所忽略[4-9]。一般情況下，高速加載(地震、爆炸、沖擊)條件下采用動(dòng)力計(jì)算，而低速(準(zhǔn)靜態(tài))加載條件下采用靜力計(jì)算。

在進(jìn)行巖石單軸壓縮實(shí)驗(yàn)時(shí)，加載速率是必須輸入的邊界條件，一般采用載荷控制(kN/min)或位移控制(mm/s)。然而，在有限元數(shù)值計(jì)算中，一般將其簡(jiǎn)化為靜力問(wèn)題，按照kN/step或mm/step進(jìn)行加載，加載速率與時(shí)間無(wú)關(guān)[10-16]。這種處理方式雖然可以實(shí)現(xiàn)巖石損傷破裂全過(guò)程模擬并得到全應(yīng)力應(yīng)變曲線，但不能真實(shí)反映物理世界的變化。它與真實(shí)物理實(shí)驗(yàn)在時(shí)間上不存在對(duì)應(yīng)關(guān)系。從嚴(yán)格意義上講，巖石力學(xué)實(shí)驗(yàn)中沒(méi)有靜力問(wèn)題，靜力問(wèn)題只是動(dòng)力問(wèn)題的簡(jiǎn)化。即便在準(zhǔn)靜態(tài)加載條件下，也應(yīng)考慮時(shí)間變量。動(dòng)力問(wèn)題和靜力問(wèn)題(準(zhǔn)靜態(tài)問(wèn)題)應(yīng)具有統(tǒng)一的表達(dá)。

數(shù)值計(jì)算的核心問(wèn)題之一便是建立反映物理世界本質(zhì)的數(shù)學(xué)模型，即建立反映真實(shí)世界各物理量之間關(guān)系的微分方程及相應(yīng)的定解條件。沒(méi)有正確完善的數(shù)學(xué)模型，數(shù)值計(jì)算便無(wú)從談起。針對(duì)動(dòng)力學(xué)問(wèn)題，人們已經(jīng)建立了運(yùn)動(dòng)微分方程，可以用來(lái)解決地震、爆炸、沖擊等問(wèn)題。對(duì)于低速加載的準(zhǔn)靜態(tài)力學(xué)問(wèn)題，運(yùn)動(dòng)微分方程及相應(yīng)的定解條件是否適用，是一個(gè)值得討論的問(wèn)題。因?yàn)閯?dòng)力計(jì)算涉及到應(yīng)力波的傳播、透射和反射，這與靜力計(jì)算完全不同，兩者計(jì)算結(jié)果是否相等鮮有文獻(xiàn)提及。本文采用數(shù)值計(jì)算，探討準(zhǔn)靜態(tài)加載條件下巖石破裂過(guò)程動(dòng)態(tài)分析的可行性，討論其時(shí)間步長(zhǎng)的選取原則。

1 準(zhǔn)靜態(tài)加載的動(dòng)態(tài)分析驗(yàn)證

假設(shè)對(duì)標(biāo)準(zhǔn)方形試件(50 mm×50 mm×100 mm)進(jìn)行單軸壓縮試驗(yàn)，采用軸向載荷控制，加載速率為30 kN/min(0.2 MPa/s)。根據(jù)上述條件，建立一個(gè)50 mm×100 mm的平面應(yīng)力模型，彈性模量取10 GPa，泊松比取0.35，密度取2 400 kg/m3。在試件兩端設(shè)置鋼墊板，墊板尺寸為50 mm×10 mm，彈性模量取200 GPa，泊松比取0.33，密度取7 850 kg/m3，模型示意圖見(jiàn)圖1。采用軸向載荷控制，上邊界施加如圖2所示的載荷，加載速率為30 kN/min(0.2 MPa/s)。下邊界采用固定約束，左右兩側(cè)采用自由邊界。

圖1 巖石單軸壓縮試驗(yàn)示意圖

圖2 邊界應(yīng)力-加載時(shí)間曲線

在瞬態(tài)求解器下，當(dāng)t=5 s時(shí)，其計(jì)算結(jié)果如圖3所示，x方向最大位移為0.89 μm，最小位移為-0.89 μm；y方向最大位移為0 μm，最小位移為-10.00 μm；x方向最大正應(yīng)力為1.18 MPa/s，x方向最小正應(yīng)力為-0.76 MPa/s；y方向最大正應(yīng)力為-0.80 MPa/s，y方向最小正應(yīng)力為-3.00 MPa/s。當(dāng)t=5 s時(shí)，對(duì)應(yīng)的上邊界應(yīng)力為1 MPa。在靜態(tài)分析模式下，其它條件不變，上邊界施加1 MPa靜應(yīng)力，其計(jì)算結(jié)果如圖4所示，x方向最大位移為0.89 μm，最小位移為-0.89 μm；y方向最大位移為0 μm，最小位移為-10.00 μm；x方向最大正應(yīng)力為1.18 MPa/s，x方向最小正應(yīng)力為-0.76 MPa/s；y方向最大正應(yīng)力為-0.80 MPa/s，y方向最小正應(yīng)力為-3.00 MPa/s。

圖3 瞬態(tài)求解器下計(jì)算結(jié)果(t =5 s)

圖5給出了兩種求解器計(jì)算的Mises應(yīng)力對(duì)比結(jié)果(豎向中心線處)，兩者完全重合。以上結(jié)果表明在軸向載荷控制加載下，動(dòng)態(tài)計(jì)算與靜態(tài)計(jì)算的結(jié)果是一樣的?？梢灶A(yù)見(jiàn)，對(duì)于任意時(shí)刻都可以得到相應(yīng)邊界條件下與動(dòng)態(tài)計(jì)算結(jié)果相同的靜態(tài)計(jì)算

圖5 豎向中心線處Mises應(yīng)力對(duì)比(t =5 s)

結(jié)果。因此，在準(zhǔn)靜態(tài)加載條件下選擇動(dòng)態(tài)求解器計(jì)算是可行的。應(yīng)力波的傳播、透射和反射對(duì)計(jì)算結(jié)果沒(méi)有影響，動(dòng)力問(wèn)題和靜力問(wèn)題(準(zhǔn)靜態(tài)問(wèn)題)具有統(tǒng)一性。

2 動(dòng)態(tài)分析的時(shí)間步長(zhǎng)討論

在真實(shí)單軸壓縮試驗(yàn)過(guò)程中，巖石是連續(xù)體，且處于連續(xù)實(shí)時(shí)變化之中。但計(jì)算機(jī)只能處理離散數(shù)據(jù)，在進(jìn)行數(shù)值計(jì)算時(shí)，必須對(duì)求解域進(jìn)行分割，離散為若干單元。涉及動(dòng)態(tài)計(jì)算時(shí)，對(duì)時(shí)間維度也必須進(jìn)行離散，設(shè)置一定的時(shí)間步長(zhǎng)。動(dòng)態(tài)分析過(guò)程中，時(shí)間步長(zhǎng)的選取是十分重要的問(wèn)題。Tedesco等[17]認(rèn)為對(duì)于動(dòng)力問(wèn)題，在數(shù)值積分中，為了保證積分精度，最大時(shí)間步長(zhǎng)必須有所限制。最大時(shí)間步長(zhǎng)的數(shù)值與應(yīng)力波速度和有限單元的尺寸有關(guān)。最大時(shí)間步長(zhǎng)的選擇要確保在此時(shí)間內(nèi)應(yīng)力波能傳過(guò)單元內(nèi)高斯積分點(diǎn)之間的距離，該時(shí)間步長(zhǎng)定義為：

(1)

式中:Lc為單元在波傳播方向的長(zhǎng)度，Cd為縱波波速。當(dāng)Δt≤1/3(Δt)max時(shí)，可以保證計(jì)算精度。

對(duì)于高速加載(高應(yīng)變率)力學(xué)過(guò)程，這種處理方式是合理的。但對(duì)于準(zhǔn)靜態(tài)加載(低應(yīng)變率)力學(xué)過(guò)程，這種處理方式有待商榷。一般巖石單軸壓縮試驗(yàn)持續(xù)時(shí)間為10 min～20 min。如果采用上述時(shí)間步長(zhǎng)，電腦內(nèi)存的消耗是巨大的，而且計(jì)算時(shí)間是漫長(zhǎng)的。以本文單軸壓縮數(shù)值計(jì)算為例，計(jì)算機(jī)處理器為英特爾酷睿i5處理器，內(nèi)存16 G，三角形單元數(shù)目20 992個(gè)，自由度數(shù)目84 802，終止計(jì)算時(shí)間5 s。圖6給出了不同時(shí)間步長(zhǎng)下計(jì)算耗時(shí)。當(dāng)時(shí)間步長(zhǎng)為1 s時(shí)，計(jì)算耗時(shí)為12.33 s；當(dāng)時(shí)間步長(zhǎng)為1×10-1s時(shí)，計(jì)算耗時(shí)為14.06 s；當(dāng)時(shí)間步長(zhǎng)為1×10-2s時(shí)，計(jì)算耗時(shí)為27.96 s；當(dāng)時(shí)間步長(zhǎng)為1×10-3s時(shí)，計(jì)算耗時(shí)為163.03 s；當(dāng)時(shí)間步長(zhǎng)為1×10-4s時(shí)，計(jì)算耗時(shí)為2 424 s；當(dāng)時(shí)間步長(zhǎng)為1×10-5s時(shí)，求解運(yùn)行24 h后，計(jì)算機(jī)內(nèi)存溢出，求解失敗。

圖6 不同時(shí)間步長(zhǎng)下計(jì)算耗時(shí)

可見(jiàn)，時(shí)間步長(zhǎng)對(duì)計(jì)算效率的影響是巨大的。所以應(yīng)該選擇既能保證計(jì)算精度，又使計(jì)算量盡量小的時(shí)間步長(zhǎng)。時(shí)間步長(zhǎng)選取不但與單元尺寸和波速有關(guān)，而且與加載速率相關(guān)。不同加載速率或應(yīng)變率(爆炸、沖擊、普通壓力機(jī)加載、蠕變)應(yīng)該選擇不同的時(shí)間步長(zhǎng)。圖7給出了單軸壓縮數(shù)值計(jì)算中不同時(shí)間步長(zhǎng)下豎向中心線的Mises應(yīng)力。從圖中可以看出，在時(shí)間步長(zhǎng)分別取1 s、0.1 s、0.01 s、0.001 s、0.000 1 s時(shí)，得到的計(jì)算結(jié)果相同。這說(shuō)明針對(duì)具體問(wèn)題，必須具體分析，時(shí)間步長(zhǎng)不能一概而論。在準(zhǔn)靜態(tài)加載條件下，選擇較大的時(shí)間步長(zhǎng)依然可以保證計(jì)算精度。

圖7 不同時(shí)間步長(zhǎng)下豎向中心線處Mises應(yīng)力(t =5 s)

3 準(zhǔn)靜態(tài)加載的高速計(jì)算與觀測(cè)

高速攝像機(jī)能將肉眼無(wú)法分辨的高速世界賦予新生命，顯示肉眼看不見(jiàn)的瞬間動(dòng)作。在巖石動(dòng)力學(xué)試驗(yàn)中，有學(xué)者成功采用高速攝影觀測(cè)到巖石的破裂過(guò)程[18]。但對(duì)于準(zhǔn)靜態(tài)加載力學(xué)過(guò)程，想要在室內(nèi)實(shí)驗(yàn)室成功實(shí)現(xiàn)高速攝影是困難的。巖石的波速一般在2 mm/s～5 mm/s，為了觀察其傳播過(guò)程，數(shù)據(jù)采集頻率必須在微秒級(jí)。否則，就無(wú)法保證有效捕捉到裂紋擴(kuò)展過(guò)程。當(dāng)攝像機(jī)運(yùn)作如此快速時(shí)，大功率的照明設(shè)備是必需的，而且攝影成員必須使用烤箱手套來(lái)操作聚焦透鏡。

如果攝像機(jī)以1幀/μs的速度記錄圖片，假設(shè)一張圖片是500 KB，乘以一百萬(wàn)就是476.84 GB，而且這只是一秒的數(shù)據(jù)，而巖石試驗(yàn)一般持續(xù)時(shí)間為10 min～20 min。因此，它需要有一個(gè)驚人的儲(chǔ)存器，而且每張圖片像素不能太高?？梢韵胂?，如此高性能的相機(jī)，技術(shù)要求頗高，而且一定價(jià)格不菲。最重要的是，攝像機(jī)只能觀測(cè)巖石表面結(jié)構(gòu)，無(wú)法觀測(cè)應(yīng)力和位移的變化。

在數(shù)值模型中，邊界條件類(lèi)似于力學(xué)試驗(yàn)機(jī)，時(shí)間步長(zhǎng)類(lèi)似于攝像機(jī)。這為我們研究準(zhǔn)靜態(tài)加載下巖石破裂問(wèn)題提供了一個(gè)很好的手段，可以得到任意時(shí)刻巖石內(nèi)的應(yīng)力和應(yīng)變分布。圖8給出了在0.2 MPa/s的加載速率下，不同時(shí)刻豎向中心線處米塞斯應(yīng)力，t表示時(shí)間步長(zhǎng)。

圖8 不同時(shí)刻豎向中心線處Mises應(yīng)力

圖8(a)的時(shí)間步長(zhǎng)為1 μs，在1 μs～5 μs時(shí)刻，巖石內(nèi)部應(yīng)力有明顯的動(dòng)態(tài)效應(yīng)，巖石內(nèi)部應(yīng)力隨著時(shí)間的推移而增加；圖8(b)的時(shí)間步長(zhǎng)為10 μs，也可以觀測(cè)到巖石內(nèi)部的動(dòng)態(tài)效應(yīng)。這種動(dòng)態(tài)效應(yīng)在應(yīng)力云圖中展示地更加清楚(見(jiàn)圖9)，應(yīng)力從上邊界逐漸向下傳遞。但圖8(c)—圖8(d)中巖石的動(dòng)力效應(yīng)消失，這是由于加載速率很低，巖石的波速約為3 313 m/s，應(yīng)力波從上邊界傳播到下邊界約需要36 μs，這段時(shí)間只能產(chǎn)生7.2 Pa的應(yīng)力增量，相比于已有的內(nèi)部應(yīng)力，其數(shù)值增量太小。因此，準(zhǔn)靜態(tài)問(wèn)題可以簡(jiǎn)化為靜態(tài)問(wèn)題。但值得注意的是，慣性效應(yīng)是否可以忽略，不但與加載速率相關(guān)，而且與研究對(duì)象的尺寸有關(guān)。如果研究對(duì)象尺寸無(wú)限大，即使在低加載速率下，慣性效應(yīng)也不能忽略。

圖9 不同時(shí)刻Mises應(yīng)力云圖

4 時(shí)間步長(zhǎng)對(duì)巖石損傷破裂的影響

需要注意的是，巖石在加載初期是連續(xù)的，但在中后期會(huì)發(fā)生損傷、破裂直至失穩(wěn)破壞，它是一個(gè)從連續(xù)到非連續(xù)的過(guò)程。如果我們用高速攝影拍攝巖石單軸壓縮試驗(yàn)過(guò)程，無(wú)論我們采用1幀/s，1幀/μs，還是1幀/ns的拍攝速度，試驗(yàn)的最終結(jié)果是一樣的，不會(huì)因拍攝速度的不同而變化。時(shí)間步長(zhǎng)類(lèi)似于拍攝速度，如果計(jì)算程序可以真實(shí)地再現(xiàn)巖石損傷與破裂全過(guò)程，那么計(jì)算結(jié)果也應(yīng)該與時(shí)間步長(zhǎng)無(wú)關(guān)。但事實(shí)上是有關(guān)的，數(shù)值計(jì)算永遠(yuǎn)都不可能做到真實(shí)，這是由計(jì)算機(jī)的特性決定的。因?yàn)楝F(xiàn)實(shí)世界的空間和時(shí)間都是連續(xù)的，而計(jì)算機(jī)只能處理離散數(shù)據(jù)，它所構(gòu)建的虛擬世界必然是非連續(xù)的。任何模擬計(jì)算都只能做到“形似”，無(wú)法做到“神似”。

作者所在的研究團(tuán)隊(duì)一直致力于巖石損傷與破裂過(guò)程的數(shù)值計(jì)算。通過(guò)考慮巖石的非均勻性，將復(fù)雜的宏觀非線性問(wèn)題轉(zhuǎn)化成簡(jiǎn)單的細(xì)觀線性問(wèn)題。通過(guò)引入數(shù)學(xué)連續(xù)物理不連續(xù)概念，將復(fù)雜的非連續(xù)介質(zhì)力學(xué)問(wèn)題轉(zhuǎn)化成簡(jiǎn)單的連續(xù)介質(zhì)力學(xué)問(wèn)題。本節(jié)應(yīng)用此方法研究時(shí)間步長(zhǎng)對(duì)單軸壓縮巖石損傷與破裂的影響。計(jì)算過(guò)程中按照位移加載，加載速率為0.002 mm/s，時(shí)間步長(zhǎng)分別取0.1 s、1 s、10 s、20 s、50 s、100 s、200 s、400 s。圖10給出了其應(yīng)變隨時(shí)間變化曲線。從圖中可以看出，巖石應(yīng)變隨著時(shí)間勻速增加。這也從另一方面說(shuō)明了動(dòng)態(tài)計(jì)算的有效性。

圖10 單軸壓縮巖石應(yīng)變隨時(shí)間變化曲線

圖11和圖12分別給出了單軸壓縮巖石應(yīng)力隨時(shí)間變化和應(yīng)力應(yīng)變曲線。因?yàn)閼?yīng)變隨著時(shí)間勻速變化，所以兩者變化規(guī)律一致。當(dāng)時(shí)間步長(zhǎng)小于10 s時(shí)，時(shí)間步長(zhǎng)對(duì)應(yīng)力曲線影響較小，計(jì)算結(jié)果不隨著時(shí)間步長(zhǎng)變化而變化。當(dāng)時(shí)間步長(zhǎng)大于10 s時(shí)，時(shí)間步長(zhǎng)對(duì)應(yīng)力曲線影響顯著，巖石的峰值應(yīng)力和彈性模量都隨著時(shí)間步長(zhǎng)的增加而減小。但在同一時(shí)刻，不同時(shí)間步長(zhǎng)下的計(jì)算結(jié)果基本重合。這是因?yàn)闀r(shí)間步長(zhǎng)取值過(guò)大，從而跳過(guò)了關(guān)鍵數(shù)據(jù)節(jié)點(diǎn)，造成曲線失真。圖13給出了單軸壓縮巖石損傷隨時(shí)間變化曲線。從圖中也得出類(lèi)似結(jié)論，當(dāng)時(shí)間步長(zhǎng)小于20 s時(shí)，不同時(shí)間步長(zhǎng)下的損傷曲線基本重合。當(dāng)時(shí)間步長(zhǎng)大于20 s時(shí)，不同時(shí)間步長(zhǎng)下的損傷曲線差別顯著。

圖11 單軸壓縮巖石應(yīng)力隨時(shí)間變化曲線

圖12 單軸壓縮巖石應(yīng)力應(yīng)變曲線

圖13 單軸壓縮巖石損傷隨時(shí)間變化曲線

因此，時(shí)間步長(zhǎng)必須小于某一閾值。本文模型峰值應(yīng)力的對(duì)應(yīng)時(shí)刻tp約為280 s，10/280≈0.04，建議選取0.03tp作為時(shí)間步長(zhǎng)。圖14給出了不同時(shí)間步長(zhǎng)下巖石破裂模式。從圖中可以看出，巖石損傷與破裂受時(shí)間步長(zhǎng)的影響較大。當(dāng)時(shí)間步長(zhǎng)小于20 s時(shí)，巖石總是沿一條剪切破裂帶破壞，其破裂模式基本不受時(shí)間步長(zhǎng)的影響。當(dāng)時(shí)間步長(zhǎng)大于20 s時(shí)，巖石破裂模式開(kāi)始發(fā)生變化，巖石呈共軛剪切破壞。時(shí)間步長(zhǎng)為50 s和100 s時(shí)，形成兩條主破裂帶。時(shí)間步長(zhǎng)為200 s、300 s和400 s時(shí)，形成多條主破裂帶。從巖石最終破裂模式也可以看出，時(shí)間步長(zhǎng)取0.04tp比較合理。

圖14 不同時(shí)間步長(zhǎng)下巖石破裂模式

5 結(jié) 論

本文利用數(shù)值計(jì)算，證明了準(zhǔn)靜態(tài)加載條件下的巖石破裂過(guò)程動(dòng)態(tài)分析的可靠性，主要結(jié)論如下：

(1) 在準(zhǔn)靜態(tài)加載條件下選擇動(dòng)態(tài)求解器計(jì)算是正確可行的。應(yīng)力波的傳播、透射和反射，對(duì)計(jì)算結(jié)果沒(méi)有影響，動(dòng)力問(wèn)題和靜力問(wèn)題(準(zhǔn)靜態(tài)問(wèn)題)具有統(tǒng)一性。慣性效應(yīng)是否可以忽略，不但與加載速率相關(guān)，而且與研究對(duì)象的尺寸有關(guān)。

(2) 在動(dòng)力計(jì)算中，時(shí)間步長(zhǎng)選取不但與單元尺寸和波速有關(guān)，而且與加載速率相關(guān)。不同加載速率或應(yīng)變率(爆炸、沖擊、普通壓力機(jī)加載、蠕變)應(yīng)該選擇不同的時(shí)間步長(zhǎng)。在準(zhǔn)靜態(tài)加載條件下，選擇較大的時(shí)間步長(zhǎng)依然可以保證計(jì)算精度。

(3) 時(shí)間步長(zhǎng)對(duì)巖石的應(yīng)力應(yīng)變曲線和損傷與破裂模式影響較大。時(shí)間步長(zhǎng)取值過(guò)大將造成曲線失真和破裂模式的不確定性。當(dāng)時(shí)間步長(zhǎng)小于0.03tp時(shí)，巖石應(yīng)力應(yīng)變曲線和損傷與破裂模式趨于穩(wěn)定。