何 峰，任天嬌，周 剛，卞紅杰，楊 松

(1.遼寧工程技術大學 力學與工程學院，遼寧 阜新 123000；2.陜西火石咀煤礦有限責任公司，陜西 咸陽 713599)

0 引言

地下煤巖體通常受到區(qū)域頂板及底板巖層的擠壓作用，這與單軸壓縮條件下的受力方向只在同一直線上，受單一作用力，產(chǎn)生壓縮變形的載荷作用近似?；诖?，導致煤體具有一定的初始損傷。由于在現(xiàn)場進行實驗的困難，通過在真實巖石和類巖石材料上預先制造裂隙的室內(nèi)模型進行單軸壓縮實驗研究裂隙巖體的破壞過程[1]。朱其志等[2]利用單軸壓縮實驗分析煤體受載內(nèi)部裂紋擴展情況，得到裂紋均是在預制裂縫尖端首先產(chǎn)生，隨后逐漸擴展貫通，最終導致試件的整體失穩(wěn)破壞，煤體整體宏觀力學性能劣化。但在分析煤樣破壞時，采用的微分方程通常會發(fā)生求解不連續(xù)、位置偏導不存在的情況，由此導致分析困難。

近場動力學(Peridynamics，PD理論)作為1種新興的非局部思想建立模型，通過求解空間積分方程來描述物質(zhì)力學行為的方法；同時還兼有分子動力學方法、鍵基理論和無網(wǎng)格方式等多種分析手段，有效避免上述空間微分方程求解不連續(xù)的奇異性，突破在計算尺度上的局限，在宏觀破壞與微觀能量及鍵的斷裂問題上均表現(xiàn)出很高的求解精度[3]。Silling等[4]提出的近場動力學方法通過損傷變量與空間積分方程來描述材料的損傷破壞情況，是1種非局部理論算法，從理論上克服局部計算時對計算網(wǎng)格的依賴，并總結巖體裂紋尖端部分數(shù)值模擬方案和內(nèi)部能量來源及擴散情況；黃丹等[5]概述PD方法的理論基礎、建模思路等，并介紹PD方法在不同尺度上、非連續(xù)力學問題中，均勻與非均勻材料和結構的大變形、損傷、斷裂、沖擊、穿透和失穩(wěn)問題的應用。繼而出現(xiàn)了很多關于函數(shù)的構建、裂紋損傷的數(shù)值模擬以及算法優(yōu)化的研究。PD理論用于分析工程材料斷裂過程中的裂紋擴展問題，構成該理論的本構關系—用鍵的斷裂表示材料斷裂[6]。從理論上克服局部化方法存在的數(shù)值結果對計算網(wǎng)格的依賴性。

基于上述研究，對煤體破壞全過程裂紋擴展情況進行單軸壓縮實驗，根據(jù)近場動力學與損傷變量演繹推導PD損傷本構力函數(shù)，進而基于PD損傷本構力模型分別在宏、微觀2個方面對煤體損傷規(guī)律進行探討，根據(jù)所構建的本構力函數(shù)模型來進行l(wèi)ammps軟件數(shù)值模擬，比對室內(nèi)實驗與模擬結果來研究煤體損傷過程及損傷規(guī)律。

1 實驗

1.1 煤體制備與實驗

樣品取自于陜西火石咀煤礦有限責任公司，共3組煤：1組靠近煤層頂板、1組為煤層中部、1組靠近煤層底板，煤樣加工為50 mm×50 mm×100 mm的長方體。每種煤樣進行5次實驗。通過萬能拉伸實驗機對不同煤樣進行單軸壓縮實驗，以0.05 mm/s的位移速率加載直至煤樣破壞，由于實驗中，部分試件裂紋不明顯或煤樣破碎，無法進行更進一步的研究，所以本次研究主要分析3組中最優(yōu)煤樣。其壓縮前、后對比如圖1所示。

根據(jù)圖1可知，煤樣都會經(jīng)歷壓密階段、彈性階段、屈服階段以及破壞階段。實驗開始時，煤樣處于微壓力初值階段，緊接著煤樣進入稍微壓密階段，隨后煤樣微節(jié)理逐漸穩(wěn)定；壓力逐漸增加，煤樣內(nèi)部微節(jié)理開始擴展；隨著壓力的進一步加大，煤樣在微裂紋的基礎上迅速擴展，導致煤樣開裂；最后壓力增大導致微裂紋相互貫通，煤樣破壞。

卸載后，分叉裂紋開始完全貫通，形成十字交叉的主裂紋，同時有多條新的裂紋出現(xiàn)，試件徹底破壞。煤樣裂紋形式以“X”型裂紋以及“Y”型裂紋為主，同時伴有其他微裂紋的產(chǎn)生及擴展。

1.2 結果分析

煤體材料內(nèi)部存在不同角度、不同尺寸、不同形狀的孔隙、裂隙。原生微裂隙的存在使煤體宏觀力學特性表現(xiàn)為非線性、非均勻的特點。研究表明，應力-應變曲線能夠很好地反應材料損傷狀況，而巖石應力-應變曲線可分為彈性階段、塑性階段、峰后軟化階段和殘余應力階段，這是巖石類材料壓縮過程中常見現(xiàn)象[7]。實驗所得煤樣應力-應變曲線與巖石類材料全應力-應變曲線對比如圖2所示。

圖2 不同煤樣應力-應變曲線

根據(jù)圖2，不同煤體的損傷與巖石類材料的應力-應變曲線有很好的對應性。另外軸向應變的變化過程大致分為：達到峰值→短暫下降→急劇上升→垂直下降。這種變化現(xiàn)象其實是巖石類材料的相繼屈服現(xiàn)象，當煤體試件初始屈服后，隨著應力和變形的增加，屈服應力不斷提高(應變硬化或強化)或者提高到一定程度后降低，這是一般煤體試件壓縮過程中常見的現(xiàn)象。煤樣的垂直應力σ越大垂直應變εy就越大，損傷變量越大煤樣的破碎效果表現(xiàn)得越好。繼而煤體的受力破壞過程，即其內(nèi)部裂隙和節(jié)理的產(chǎn)生、發(fā)展的過程越快。

2 近場動力學(PD)理論

2.1 PD理論

PD理論是在傳統(tǒng)連續(xù)介質(zhì)力學基礎上對固體力學中的運動方程的重新描述，該理論的積分方程代替?zhèn)鹘y(tǒng)力學中微分方程，應用于求解非連續(xù)位移處，避免了微分方程出現(xiàn)奇異解的問題。與傳統(tǒng)的連續(xù)力學相比，PD理論是物質(zhì)點在其鄰域內(nèi)與其他物質(zhì)點之間發(fā)生相互作用?？梢悦枋鑫灰频倪B續(xù)、非連續(xù)變化；模擬非連續(xù)位移處裂紋的損傷擴展過程。同時PD理論具有其本身自己的破壞準則，在求解損傷破壞問題時，無需引入其他的破壞準則，通過求解積分方程的方式描述物體的整個損傷斷裂過程，裂紋可以沿任意方向擴展。因此，PD理論在求解非連續(xù)問題時展示出其獨特優(yōu)勢，能夠精確模擬復雜的損傷和裂紋擴展過程[8-9]。

將空間R離散為點單元，其中某質(zhì)點x在t時刻的基本方程如式(1)所示：

(1)

式中:ρ(x)表示t時刻質(zhì)點x的密度,g/cm3；ü(x,t)表示t時刻質(zhì)點x的加速度，m/s2；u(x,t)與u(x′,t)表示物質(zhì)點x與x′的位移，m；b(x,t)表示物質(zhì)點x受到的外力，MPa；函數(shù)f是點x與x′之間的相互作用力函數(shù)，僅與物質(zhì)點對x與x′之間鍵的變形有關[10]；H為位置點的作用域。

根據(jù)式(1)可知，物質(zhì)點x與x′在t時刻分別產(chǎn)生u(x,t)與u(x′,t)的位移，2質(zhì)點在初始時刻的相對位置向量表示為：x′-x=ξ，并經(jīng)過時間t之后2點間的相對位移向量為：u(x′,t)-u(x,t)=η，所以發(fā)生變形之后的相對位置向量為：ξ+η。此時，物質(zhì)點x在區(qū)域內(nèi)其他點的共同作用下而保持運動狀態(tài)或平衡狀態(tài)，其示意圖如圖3所示。

圖3 近場動力學原理示意

PD理論中使用的積分方程替代傳統(tǒng)力學中的微分方程，避免在非連續(xù)位移處導數(shù)不存在的問題，解決了煤體在損傷失效過程中產(chǎn)生的非連續(xù)問題。PD模型其中一種特例就是鍵基近場動力學模型，在該模型中，假設物質(zhì)點之間的力密度矢量僅與物質(zhì)點之間相對位置有關，而與物質(zhì)點近場區(qū)域內(nèi)的其他物質(zhì)點變形無關。2物質(zhì)點之間存在大小相等、方向相反、滿足牛頓第三定律的力密度矢量，根據(jù)這種特性，鍵基近場動力學模型在保證計算精度的準確性情況下會一定程度上減少計算量，在解決復雜裂隙擴展方面尤為突出，其理論示意如圖4所示。

圖4 鍵基PD模型示意

該理論將模型區(qū)域離散成有限的物質(zhì)點，在以δ為半徑的近場區(qū)域H內(nèi)物質(zhì)點間存在著相互作用力——近場力。

其中，Hx是x的δ領域(近場區(qū)域的大小)，是Hx內(nèi)的點，定義在與物質(zhì)點x距離為δ以內(nèi)的點才會與其中心點x發(fā)生相互作用，與x有相互作用的點可以用集合Hx表示，如式(2)所示：

(2)

通過鍵的斷裂，PD理論能夠完整地模擬煤體的全部損傷失效過程。鍵破壞模式表示在1個物質(zhì)點對內(nèi)，相互作用的2個物質(zhì)點間經(jīng)過一定的相對位移后，當點對之間鍵的伸長率s超過臨界伸長率s0時，該點對破壞，2物質(zhì)點間不再有相互作用，將發(fā)生損傷。鍵的斷裂和裂紋的關系示意如圖5所示。

圖5 鍵的斷裂和裂紋的關系示意

若將物質(zhì)點對鍵看作中心彈簧的鉸接桿單元，可得本構力表達形式如式(3)所示：

(3)

(4)

近場力與鍵伸長率的關系如圖6所示。當鍵的伸長率s超過s0時，鍵發(fā)生斷裂，物質(zhì)點之間不再存在力的作用。

圖6 近場力與鍵伸長率關系

由于鍵型PD模型本構力函數(shù)包含損傷和斷裂條件，不需要引入其他的破壞準則，如傳統(tǒng)的臨界應力強度因子等，運動方程中通過把力密度矢量用不可逆的方式表達損傷，荷載會在物質(zhì)點之間產(chǎn)生重新分配，導致宏觀損傷以自動發(fā)展的方式逐漸擴展。因此，在分析破壞問題時不必進行傳統(tǒng)的開裂判斷和裂紋路徑等分析，裂紋將自然萌生并擴展。修正后的f(η,ξ)函數(shù)很好地體現(xiàn)了“近場”的思想，在近場范圍內(nèi)物質(zhì)點對之間存在相互作用?，F(xiàn)有的損傷模型和方程都是在靜態(tài)實驗條件下完成，得到的大多是宏觀的規(guī)律，缺乏可靠的微觀實驗基礎[11]。該模型中，本構力函數(shù)f與物質(zhì)點對的伸長率s線性相關，這里的伸長率s可以看作是連續(xù)介質(zhì)力學中的“應變”，故該模型可以用于宏觀領域，即能在微觀和宏觀上達到統(tǒng)一。

2.2 脆性材料本構力函數(shù)模型

前蘇聯(lián)學者Kachanov[12]定義的損傷變量如式(5)所示：

(5)

式中：D為損傷變量；A*為煤體受損后損傷面積，m2；A為煤體初始無損橫截面積，m2。

根據(jù)PD本構力函數(shù)f和損傷變量D通過數(shù)學擬合得到式(6)，脆性材料本構力函數(shù)曲線如圖7所示。

圖7 脆性材料本構力函數(shù)曲線

(6)

式中：sc，st分別為斷裂時的壓縮率和伸長率，即當伸長率s超過斷裂臨界值時，材料將發(fā)生破壞；s0c，s0t分別為臨界壓縮率和臨界伸長率，即壓縮或拉伸狀態(tài)下煤體試件不會產(chǎn)生裂紋的最大伸長率；sec，set分別為彈性壓縮率和彈性伸長率，當sec≤s≤set時，材料受拉壓表現(xiàn)為彈性性狀，在不受力時，材料會恢復原狀。

由圖7可知，陰影區(qū)為鍵的損傷演化階段，呈非線性。②，③陰影區(qū)分別表示材料在受拉伸或壓縮時產(chǎn)生部分鍵的斷裂并伴有微裂紋的產(chǎn)生；①，④陰影區(qū)分別表示材料在受拉伸或壓縮時，產(chǎn)生鍵的斷裂并有宏觀裂紋的產(chǎn)生和擴展。在該模型中，物質(zhì)點之間的作用力隨著伸長率增大先增大后減小，力在到達臨界伸長率時鍵發(fā)生破壞，而不是到達峰值作用力時突然發(fā)生破壞，符合巖石類材料先應變硬化再應變軟化的變形機制[13-14]。

在式(6)伸長率s與鍵的損傷變量D關系中，當伸長率小于線彈性伸長率時，鍵不會出現(xiàn)變形斷裂；當伸長量大于彈性伸長率時，鍵的損傷會隨著伸長率呈指數(shù)增長直至斷裂，鍵的斷裂意味著裂紋的形成。具有線性和非線性的力學行為的本構力函數(shù)的基本形式如式(7)所示：

(7)

式中：s可根據(jù)物質(zhì)點的相對位移η和相對位置ξ求出。

由于煤體受到外荷載作用后內(nèi)部微結構會產(chǎn)生不同程度的損傷，損傷產(chǎn)生過程中必然有不同程度上微觀的鍵的斷裂。因此，在PD理論的基礎上，結合損傷變量改進PD模型，進而對煤體損傷規(guī)律進行探究和總結是可行的。

3 數(shù)值模擬

以近場動力學的角度利用lammps原子分子并行模擬器來模擬并建立無網(wǎng)格的PD離散型模型，對于巖石類材料的沖擊模擬，將模型設定為50 mm×50 mm×100 mm的長方體試樣。密度、泊松比、彈性模量等材料參數(shù)均使用數(shù)值軟件所定義的參數(shù)。為便于計算，將總時間步設定為20 000步。將模擬試樣離散為250 000物質(zhì)點，物質(zhì)點間距為Δ=0.1 mm，δ=3.015Δ[15-16]。下邊界固定，上邊界以0.1 mm/min的速度施加均布荷載，為了使模擬結果更為準確，增加馳豫時間以卸載初始應力使模擬試樣從激化的狀態(tài)回復到平衡排列狀態(tài)，強化馳豫過程，觀察其裂紋擴展及內(nèi)部能量的轉換情況。

在模擬過程中采用PD模型來模擬巖體材料的力學特性，采用式(1)作為數(shù)值軟件中的動力學控制方程，式(7)作為本構力函數(shù)。

利用compute stress/atom命令所計算出來的應力是通過對Virial定理所定義的體系應力按原子分解，即將Virial定理中按原子求和的算符刪除(同時應將體系的體積換為單個原子的體積)。但由于單個原子的體積計算太麻煩，所以lammps計算時直接去掉了體積項，導致所算出來的“應力”具有能量的單位。因此，所輸出應力需經(jīng)過換算才能與實際符合，利用“compute v all voronoi/atom”等命令進行后續(xù)換算，對原子應力進行Voronoi體積加權平均即可得到系統(tǒng)瞬時應力，系統(tǒng)瞬時應力的系綜平均值為宏觀測量的系統(tǒng)應力值。換算后煤體內(nèi)部產(chǎn)生破壞的最大應力約為0.1 MPa。經(jīng)過后續(xù)可視化等分析處理可以得到模擬結果，如圖8～11所示。

圖8 壓力閉合階段模擬損傷云圖

第1階段：壓力閉合階段。由圖8可以看出，模擬開始時，模型處于壓密階段，模型內(nèi)部微裂隙、微孔隙逐漸閉合，模型所能承載的最大應力主要集中在模型中部位置，邊緣位置應力較小，煤體模型中部出現(xiàn)了應力集中的現(xiàn)象。此時，模型中部位置發(fā)生輕微變形，宏觀表現(xiàn)不明顯。

圖9 線彈性變形階段損傷云圖

第3階段：非線性變形階段。由圖10可以看出，模型應力分布逐漸穩(wěn)定，呈“X”分布。隨著壓力逐漸增加，模型內(nèi)沿著應力分布區(qū)域大量的鍵發(fā)生斷裂，模型產(chǎn)生損傷。此時模型宏觀表現(xiàn)出有裂紋的產(chǎn)生，伸長、匯合、交叉、貫通，主裂紋呈“X”狀，并在主裂紋附近有很多微小裂紋產(chǎn)生。從煤樣宏觀角度分析，大量裂紋迅速擴展，出現(xiàn)“X”型裂紋，導致模型開裂，裂紋走向及形狀基本穩(wěn)定。

圖10 非線性變形階段損傷云圖

第4階段：宏觀破壞階段。由圖11可以看出，隨著壓力的增大，“X”型裂紋內(nèi)所積聚的應力達到極限狀態(tài)，無法繼續(xù)承載多余應力，微觀上表現(xiàn)為“X”型裂紋周圍鍵全部斷裂，多余應力貫徹遍布整個模型，宏觀上損傷演化加劇，達到屈服應力，裂紋繼續(xù)發(fā)展、交叉、相互匯合、貫通，導致試樣破壞。

圖11 宏觀破壞階段損傷云圖

圖8～11可以直觀反映煤體試樣損傷破壞全過程，試件微觀鍵基變形損傷規(guī)律、應力擴散情況、宏觀裂紋擴展情況能夠較好地對應。宏觀裂紋的萌生與擴展是微觀鍵的斷裂的外在表現(xiàn)，微觀上鍵的斷裂是宏觀裂紋發(fā)生的內(nèi)在原因。根據(jù)模擬方案，模擬煤體在單軸壓縮條件下其內(nèi)部應力-應變變化趨勢，模擬應力-應變曲線如圖12所示。

圖12 模擬應力-應變曲線

與實驗結果對比，實驗結果曲線與模擬曲線走勢大體相同，在受壓閉合階段及線彈性階段均有回彈，表明在受壓過程中均有彈性變形的發(fā)生與恢復，該階段不發(fā)生損傷。在進入非線彈性階段后應力突然下降，即大量鍵產(chǎn)生斷裂，導致試樣產(chǎn)生“X”型宏觀破壞，此時試樣發(fā)生損傷。印證了實驗的準確性與模擬的真實性，說明近場動力學能夠反映煤體的物理力學機制，并且作為一種脆性材料，煤體的本構關系與所構建的脆性本構力函數(shù)關系相同。

4 結論

1)煤體受到外荷載作用后內(nèi)部微結構會產(chǎn)生不同程度的損傷，損傷產(chǎn)生過程中有不同程度上微觀的鍵的斷裂。利用演繹推導出的本構力模型來反演煤體內(nèi)部微結構變化從而得到煤體損傷規(guī)律切實可行。

2)煤體損傷過程可分為：壓密階段、微裂紋的萌生擴展階段、斷裂破壞階段3個階段，除受壓閉合階段其余階段均伴隨鍵的斷裂，至破壞階段時，鍵幾乎完全斷裂，導致煤樣失穩(wěn)破壞。

3)模擬實驗中，模型在受到外載荷作用下，內(nèi)部能量開始聚集，隨著壓力的增大，內(nèi)部能量散開，最終呈“X”型分布，這與實驗加載后煤樣的損傷效果基本一致，可為地下采煤及沖擊地壓防治提供理論支撐。