壓頭作用下巖石破碎過程的細觀模擬

李守巨，李 德，于 申

(1.大連理工大學 工業裝備結構分析國家重點實驗室，遼寧 大連 116024；2.哈爾濱電機廠有限責任公司，黑龍江 哈爾濱 150000)

?

壓頭作用下巖石破碎過程的細觀模擬

李守巨1，李 德2，于 申1

(1.大連理工大學 工業裝備結構分析國家重點實驗室，遼寧 大連 116024；2.哈爾濱電機廠有限責任公司，黑龍江 哈爾濱 150000)

采用平行粘結細觀本構模型，在巖石顆粒尺度上研究了壓頭侵入巖石的破壞過程。分析討論了巖石材料平行粘結模型細觀參數與宏觀力學參數之間的映射關系。為了實現表征顆粒連接破壞模式，應用FISH語言編寫一個伺服程序，檢測每個顆粒連接的受力狀態。依據所建立的顆粒連接破壞準則，并對每個顆粒連接進行破斷與否的評估，實現描述顆粒連接的破壞功能。研究表明，在壓頭垂直壓力作用下，巖石試件內形成了3條主環向拉力鏈，它與宏觀的小主應力方向一致。環向拉力鏈最大影響深度隨著壓頭壓入深度的增加而增加。與此同時，出現了幾條徑向壓力鏈，它與宏觀的大主應力方向基本一致。巖石試件裂紋是由環向拉應力和剪應力大于顆粒平行連接膠結材料的強度引起的。

巖石試件；破碎過程；參數估計;細觀模擬；壓頭侵入；平行粘結模型

0 引 言

離散元方法是當今非連續介質力學數值方法中應用于巖石和混凝土等脆性材料的最廣泛而且最有效的方法之一。其最大特點是可以反映巖石塊體之間的接觸面滑移、分離以及傾覆旋轉，表征巖體內部節理、裂隙和斷層等非連續特性，同時還能夠計算塊體內部的變形和應力，尤其對于描述巖石破壞過程的細觀機制具有內在的優越性。已經廣泛應用于巖石邊坡穩定性分析、地下廠房圍巖穩定性和破壞過程模擬、混凝土結構的裂縫擴展、破壞與倒塌過程分析、盾構機滾刀破巖過程分析和混凝土大壩及其基礎破壞機制等巖石工程領域。Gu[1]提出了剛體-彈簧模型模擬混凝土材料初始缺陷和骨料分布特性，所建立的離散元模型包括粗骨料、細骨料、砂漿和粘結界面，并且考慮了這些材料的隨機分布特性。James[2]研究了非均質巖石細觀裂紋的分布特性，提出了顆粒斷裂的破壞準則。Gray[3]回顧了地質結構體的有限元和離散元建模問題，分析了有限元、離散元和有限差分法的在表征非連續地質體力學特性的適用性。Shmulevich[4]應用離散元方法(PFC)模擬了切刀與土體之間的相互作用特性，研究表明該方法模擬的結果與實驗結果基本一致，其目的在于通過優化切刀的形狀，提高切刀破土的能量利用率。Ono[5]采用三維離散元方法模擬了切土工具與土體之間的作用機制，討論了6種離散元形狀對數值模擬結果的影響，并將實驗結果與模擬結果進行了對比分析。Ucgul[6]采用三維離散元模型模擬了切刀切土的破壞過程，重點分析了無內聚力散體砂土模型參數估計問題，討論了線性和非線性接觸模型的差異性和適用性。Arena[7]提出了一種新方法用于定量計算實驗過程中巖石細觀裂紋，該方法以二維數字圖像為基礎，對于定量分析巖石試件表面裂紋的演化過程具有借鑒和參考價值。Mak[8]介紹了用于表征土體與切刀之間相互作用的離散元模型，討論了細觀模型中某些參數的簡化估計方法，并且進行了切刀切土的實驗和模擬分析研究。Asaf[9]提出了如何根據切刀破土實驗數據進行離散元模型參數校準問題。Innaurato[10]采用實驗和數值模擬方法研究了TBM滾刀的破巖過程，建立了滾刀貫入深度與貫入力之間的關系，研究了滾刀作用下巖石的剪切破壞問題。從以上相關研究工作中可以看出，采用離散元方法研究土體材料破壞過程所進行的工作較多，而研究巖石材料細觀破壞過程的工作相對較少。與此同時，如何較為準確估計巖石材料的細觀模型參數是制約離散元方法解決工程實際問題的瓶頸。巖石的力學行為由其內在的細觀結構和模型參數決定，尤其是巖石細觀裂紋及其演化和發展過程。筆者的目的在于在巖石顆粒尺度上研究壓頭作用下的巖石破壞過程，揭示巖石破碎漏斗的細觀形成機制。

1 巖石材料平行粘結模型細觀參數估計

巖石是多種礦物和膠結物組成的一種高度不均勻材料，在細觀層次上表現在顆粒、粘結物和結構面等力學特性的隨機分布特點。顆粒流(Particle Flow Code)方法是由Cundall和Strack提出的模擬球形顆粒運動與相互作用的細觀離散元方法。該方法已經應用到巖石、土體和混凝土材料的力學特性研究中，著重從細觀力學角度揭示材料的損傷和斷裂機理。顆粒接觸剛度模型包括線性顆粒接觸剛度模型和Herz-Mindlin非線性接觸模型等。如圖1所示，顆粒接觸剛度模型參數包括法向剛度(kn)、切向剛度(ks)和摩擦系數(f)以及阻尼(η)等。

圖1 顆粒接觸剛度模型Fig.1 Contact model of particles

圖2 巖石材料平行粘結模型Fig.2 Parallel bond model for rock materials

平行粘結模型(Parallel bond model)用于模擬巖石和混凝土材料(非散體材料)的開裂問題，Azevedo提出了其簡化模式[11]，如圖2所示，巖石的細觀模型參數還包括：顆粒之間膠結材料的法向剛度(kpn)和切向剛度(kps)、膠結材料切向剪切強度(τs)和法向抗拉強度(σt)，以及膠結圓環的半徑因子(λ)。需要說明的是，顆粒之間膠結材料的切向剪切強度(τs)和法向抗拉強度(σt)與巖石宏觀實驗的剪切強度和抗拉強度并不是一個概念，但二者之間應該具有某種相關性，因為平行粘結模型將顆粒之間的粘結簡化成了一個等效的圓盤。

平行粘結模型法向剛度kn為[12]

kn=2tEc,

(1)

式中t為寬度，對于2維模型，t=1；Ec為顆粒接觸變形模量，近似取為巖石宏觀實驗的彈性模量E，取Ec=43 GPa.切向剛度ks與法向剛度的關系為

ks=Rckn,

(2)

式中Rc為顆粒切向剛度與法向剛度的比，研究表明，這兩個值近似相等[13]，取Rc=1.0.顆粒之間的細觀摩擦系數f與宏觀內摩擦角φ的近似關系為[14]

(3)

顆粒之間膠結材料的法向剛度kpn為

(4)

式中 Epc為平行粘結變形模量；RA和RB分別為顆粒A和B的半徑。平行粘結顆粒之間膠結材料的切向剛度kps為

kps=Rpkpn,

(5)

式中 Rp為平行粘結模型切向剛度與法向剛度的比，取Rp=1.平行粘結模型細觀抗拉強度與宏觀內聚力c和內摩擦角φ的關系為[8]

σt=c×cotφ,

(6)

式中 σt為平行粘結模型細觀法向抗拉強度。平行粘結細觀抗剪強度與宏觀內聚力的關系為[8]

τs=c,

(7)

式中 τs為顆粒之間膠結材料的切向剪切強度；顆粒之間固定平行粘結半徑因子取為λ=0.5.顆粒之間粘結環半徑Rb為

(8)

模擬的巖石試樣為中晶大理巖，其粒徑為1～4 mm，顆粒流模型顆粒的直徑為2 mm，巖石試樣的宏觀力學參數見表1，表2為巖石細觀參數估計值。

表2 巖石試樣的細觀參數估計值

2 壓頭作用下巖石破碎過程細觀模擬

巖石是由多種礦物和膠結物組成的非均勻體，并且含有大量的天然節理和微裂隙。天然巖體中的斷層、節理和裂隙的隨機分布，其空間幾何和力學參數的分布具有統計特征。在長期的地質作用和水文環境下巖體產生的微裂紋、缺陷和孔隙，在外界荷載作用下，在巖石礦物顆粒和膠結物結合面上存在初始微裂紋損傷，在破壞機制上存在局部化現象，并導致應變軟化的本構特性[14]。

在單軸壓縮情況下，巖石試樣的破壞過程是隨機分布強度參數的巖石顆粒的剛度逐漸劣化的結果。從連續介質力學觀點分析，巖石試樣的破壞方式就其力學本質可以歸結為壓縮破壞、拉伸破壞和剪切破壞。為了實現表征顆粒連接力學破壞模式，應用FISH語言編寫一個伺服程序，檢測每個顆粒連接的受力狀態，依據所建立的顆粒連接力學破壞準則，并對每個顆粒連接進行破斷與否的評估，進而實現描述顆粒連接的破壞功能。

圖3 平行粘結模型法向力與法向位移之間的本構關系Fig.3 Constitutive relation between normal force and displacement for parallel bond model

圖4 平行粘結模型切向力與切向位移之間的本構關系Fig.4 Constitutive relation between tangential force and displacement for parallel bond model

如圖3所示，顆粒之間法向接觸力Fn與法向相對位移之間的關系為

Fn=knUn，

(9)

式中Un為法向相對接觸位移，即2個球的重疊量，其為割線模量，kn建立了總體法向力與總體法向位移之間的線性關系。

如圖4所示，切向接觸力以增量形式計算，其計算公式為

ΔFs=-ksΔUs,

(10)

式中ks為切向剛度，其為切線模量；ks建立了切向力增量ΔFs與切向位移增量ΔUs之間的線性關系。顆粒粘結的拉應力破壞準則為

葉總拿過玉墜仔細把玩了一番，又認真咀嚼起王祥的話。錢總則是和王祥一起納起悶來。良久，葉總突然笑了起來，他把玉墜遞還給王祥：

(11)

式中Ftmax為顆粒粘結所承受的極限拉力。顆粒粘結的剪應力破壞準則為

(12)

式中 Fsmax為顆粒粘結所承受的極限剪力。當顆粒粘結剪切破壞后，并且顆粒之間處于受壓狀態時，其滑動阻力Fsf為

Fsf=fFn(Fn>0).

(13)

如圖5所示，金屬壓頭彈性模量為200GPa，剛性壓頭利用墻體模擬，設定壓頭尺寸為20mm×20mm，賦予壓頭的法向剛度取為400GN/m，侵入速度為0.001mm/step，總侵入深度為10.0mm.巖石試樣的幾何尺寸為200mm×100mm，巖石顆粒密度為2 700kg/m3.模型邊界利用墻體模擬[15]，令其法向剛度等于顆粒法向剛度。設定巖石顆粒粒徑為2.0mm，顆粒數目為5 000個。隨著巖石試樣的加載過程，巖石材料將不斷產生新的裂紋，部分顆粒之間的粘結將發生斷裂力學形式的破斷，某些顆粒配位數也將隨之不斷變化。

圖5 壓頭侵入深度2 mm時的力鏈分布Fig.5 Force chains in rock specimen under 2 mm penetration depth

圖6 壓頭侵入深度4 mm時的力鏈分布Fig.6 Force chains in rock specimen under 4 mm penetration dept

在平行粘結力鏈圖中，深色部分表示顆粒之間處于受拉狀態，黑色部分表示顆粒之間處于受壓狀態。從圖5中可以看出，在壓頭垂直壓力作用下，在巖石試件內，形成了3條主環向拉力鏈，它與宏觀的小主應力方向一致，環向拉力鏈最大影響深度約為30mm.同時，出現了幾條徑向壓力鏈，它與宏觀的大主應力方向基本一致。

從圖6中可以看出，隨著壓頭侵入深度的增加，環向拉力鏈最大影響深度擴展到約為42mm；主環向拉力鏈發展得更加明顯。

圖7 壓頭侵入深度6 mm時的力鏈分布Fig.7 Force chains in rock specimen under 6 mm penetration depth

圖8 壓頭侵入深度8 mm時的力鏈分布Fig.8 Force chains in rock specimen under 8 mm penetration depth

從圖7中可以看出，隨著壓頭侵入深度的增加，環向拉力鏈最大影響深度擴展到約為56mm；在壓頭下部附近出現了一部分破碎區。

對比圖7，8和9可以發現，壓頭下面的破碎區域(淺色區域)隨著壓頭侵入深度的增加而逐漸增加，直至形成破碎漏斗。

圖9 壓頭侵入深入10 mm時的力鏈分布Fig.9 Force chains in rock specimen under 10 mm penetration depth

從圖9中可以發現，在壓頭正下方形成了一個明顯的破碎漏斗(淺色區域)，破碎漏斗的角度為70左右。巖石是一種典型的非均質和準脆性材料[16]，巖石晶粒和初始缺陷的分布具有隨機性和不確定性，在荷載作用下，各組分成分的剛度差異很大，使得對力的傳遞響應和自身變形的差異，導致在巖石內部應力場和應變場分布極不均勻，產生局部應力集中，在初始缺陷或弱結合面區域產生微裂紋，且裂紋的擴展路徑具有明顯的非對稱性。巖石中裂紋擴展路徑強烈依賴于巖石材料的非均勻性，材料的非均勻性對裂紋的貫通和抑制有著重要的影響。在細觀尺度上，巖石內部所含有的晶粒和膠結物的力學性質差異較大，巖石材料具有明顯的非均勻性。

3 結 論

1)巖石材料在彈塑性變形直至破壞過程中，其宏觀體積響應只是各種微結構之間相互作用的平均結果，變形的主要機制是微裂紋成核和增長及其聚集而生成細觀裂紋、擴展成宏觀裂紋、直至宏觀破壞。由于巖石(巖體)內含有大量的初始缺陷(節理、層理、微裂紋等)，連續介質力學方法對于表征巖石的這些非連續特性存在著固有的缺點，尤其是分析裂紋群擴展過程時遇到的困難更加突出，離散元方法為解決這類難題提供了新的途徑；

2)采用非連續介質力學數值方法，發展巖石的細觀力學模型有助于從材料破壞機理的層面上對宏觀巖石破壞過程進行理解和解釋。研究表明，在壓頭垂直壓力作用下，在巖石試件內，形成了3條主環向拉力鏈，它與宏觀的小主應力方向一致。環向拉力鏈最大影響深度約隨著壓頭壓入深度的增加而增加。與此同時，出現了幾條徑向壓力鏈，它與宏觀的大主應力方向基本一致；

3)巖石試件裂紋是由環向拉應力和剪應力大于顆粒平行連接膠結材料的強度引起的。壓頭作用下巖石的破壞特性細觀模擬分析成果，為TBM滾刀破巖分析、滾刀優化布置和滾刀形狀設計提供了一種可以借鑒的方法；

4)如何精確定量確定巖石細觀模型中的參數是制約離散元方法發展和應用的瓶頸問題之一。由于巖石顆粒大多在mm量級，在該尺度上直接測量力學參數是非常困難的。如果能夠在巖石試樣宏觀尺度上進行實驗，并且基于宏觀實驗數據進行力學細觀模型參數反演，將是一項十分有意義的工作。

References

[1]GUXiang-lin，LIHong，WANGZhuo-lin，etal.Amodifiedrigid-body-springconcretemodelforpredictionofinitialdefectsandaggregatesdistributioneffectonbehaviorofconcrete[J].ComputationalMaterialsScience，2013，77(3)：355-365.

[2]JamesPM，FordM，JivkovAP.Anovelparticlefailurecriterionforcleavagefracturemodelingallowingmeasuredbrittleparticledistributions[J].EngineeringFractureMechanics，2014，121-122(1)：98-115.

[3]GrayGG，MorganJK，SanzPF.Overviewofcontinuumandparticledynamicsmethodsformechanicalmodelingofcontractionalgeologicstructures[J].JournalofStructuralGeology，2014，59(2)：19-36.

[4]ShmulevichI，AsafZ，RubinsteinD，Interactionbetweensoilandawidecuttingbladeusingthediscreteelementmethod[J].Soil&TillageResearch，2007，97(1)：37-50.

[5]OnoI，NakashimaH，ShimizuH，etal.Investigationofelementalshapefor3DDEMmodelingofinteractionbetweensoilandanarrowcuttingtool[J].JournalofTerramechanics，2013，50(4)：265-276.

[6]UcgulM，FielkeJM，Saunders.Three-dimensionaldiscreteelementmodelingoftillage：Determinationofasuitablecontactmodelandparametersforacohesionlesssoil[J].BiosystemsEngineering，2014，121(2)：105-117.

[7]ArenaA，PianeCD，SaroutJ.Anewcomputationalapproachtocracksquantificationfrom2Dimageanalysis：Applicationtomicro-cracksdescriptioninrocks[J].Computers&Geoscience，2014，66(1)：106-120.

[8]MakJ，ChenY，SadekMA.Determiningparametersofadiscreteelementmodelforsoil-toolinteraction[J].Soil&TillageResearch，2012，118(5)：117-122.

[9]AsafZ，RubinsteinD，ShmulevichI.Determinationofdiscreteelementmodelparametersrequiredforsoiltillage[J].Soil&TillageResearch，2007，92(1-2)：227-242.

[10]InnauratoN，OggeriC，OrestePP，etal.ExperimentalandnumericalstudiesonrockbreakingwithTBMtoolsunderhighstressconfinementrockmech[J].RockMechanicsandRockEngineering，2007，40(5)：429-451

[11]AzevedoNM，LemosJV，AlmeidaJR.Influenceofaggregatedeformationandcontactbehaviorondiscreteparticlemodelingoffractureofconcrete[J].EngineeringFractureMechanics，2008，75(6)：1 569-1 586.

[12]PotyondyDO，CundallPA.Abonded-particlemodelforrock[J].InternationalJournalofRockMechanicsandMiningSciences，2004，41(8)：1 329-1 364.

[13]李守巨，李 德，于 申.基于宏觀實驗數據堆石料細觀本構模型參數反演[J].山東科技大學學報，2015，34(5)：20-26.

LIShou-ju，LIDe，YUShen.Meso-parametersInversionofconstitutivemodelforrock-fillmaterialsbasedonmicroexperimentaldata[J].JournalofShandongUniversityofScienceandTechnology，2015，34(5)：20-26.

[14]孫其誠，厚美瑛，金 峰.顆粒物質物理與力學[M].北京：科學出版社，2011.

SUNQi-cheng，HOUMei-ying，JINFeng.Physicalandmechanicalsofgranularmaterial[M].Beijing：SciencePress，2011.

[15]喬衛國，楊 麟，李大勇.高頻振動沉樁的離散元模擬分析[J].西安科技大學學報，2012，32(5)：604-609.

QIAOWei-guo，YANGLin，LIDa-yong.Numericalanalysisofpilesdrivingbyhighfrequencyvibrationusingdiscreteelementmethod[J].JournalofXi’anUniversityofScienceandTechnology，2012，32(5)：604-609.

[16]馬文偉，趙光明，孟祥瑞.巖石SHPB劈裂分析理論研究[J].西安科技大學學報，2013，33(5)：283-290.

MAWen-wei，ZHAOGuang-ming，MENGXiang-rui.TheoreticalstudyforsplittinganalysisofrockusingSHPB[J].JournalofXi’anUniversityofScienceandTechnology，2013，33(5)：283-290.

Meso-simulation for fracturing process of rock specimen under action of indenter

LI Shou-ju1，LI De2，YU Shen1

(1.StateKeyLaboratoryofStructuralAnalysisforIndustrialEquipment，DalianUniversityofTechnology，Dalian116024，China；2.HarbinElectricMachineryCompanyLimited，Harbin15000，China)

The parallel bond model in discrete element method is used to simulate the fracturing process of rock specimen under action of indenter.The nonlinear mapping relations between meso-parameters of constitutive model and macro mechanical parameters of rock specimen are analytically bridged.In order to characterize the fracture mode between the bonds of parcels each other，the serving program is coded with FISH language，and the stress states for every bound is checked.Based on the proposed fracturing criterion for the bonds of parcels，the evaluation for whether exists fracture between the bonds of parcels is performed.The investigation shows that under the action of vertical indenter，the tree tensile force chains are produced in circumferential direction，and their directions are same as the directions of the minor principal stress.The directions of radial compressive force chains are same as the directions of the major principal stress.The cracks in rock specimen are induced by tensile force in circumferential direction and shear force when the tensile stress or shear stress is larger than material strength.

rock specimen；fracturing process；parameter evaluation；meso-simulation；indenter penetration；parallel bond model

10.13800/j.cnki.xakjdxxb.2016.0602

1672-9315(2016)06-0769-06

國家重點基礎研究發展計劃(2015CB057804)；國家自然科學基金(11572079)；工業裝備結構分析國家重點實驗室開放基金(S14206)

李守巨(1960-)，男，遼寧沈陽人，教授，博士，E-mail:lishouju@dlut.edu.cn

2016-05-10 責任編輯：楊泉林

TU 443

A