SPH法切削破巖數(shù)值仿真

歐陽(yáng)義平，楊啟，2*，程旭東

（1.上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240；2.上海交大海洋水下工程科學(xué)研究院有限公司，上海 200231；3.中國(guó)鐵建港航局集團(tuán)有限公司，廣東 珠海 519070）

SPH法切削破巖數(shù)值仿真

歐陽(yáng)義平1，楊啟1，2*，程旭東3

（1.上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240；2.上海交大海洋水下工程科學(xué)研究院有限公司，上海 200231；3.中國(guó)鐵建港航局集團(tuán)有限公司，廣東 珠海 519070）

借助LS-DYNA計(jì)算軟件，用SPH法求解偏微分方程進(jìn)行切削破巖數(shù)值仿真，并以不同切削參數(shù)的切削破巖實(shí)驗(yàn)為原型進(jìn)行計(jì)算，得出切削力與實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。數(shù)值仿真結(jié)果與實(shí)際相符，刀齒附近巖樣單元應(yīng)力和塑性應(yīng)變均最大，切削過后巖樣上形成切削槽，兩側(cè)存在破碎角，比以往數(shù)值仿真更接近實(shí)際情況，且破碎角大小與剪脹角有密切關(guān)系。大量仿真結(jié)果顯示，計(jì)算所得切向力無論變化規(guī)律還是數(shù)值均與實(shí)驗(yàn)吻合較好，尤其是平均切向力，證實(shí)了SPH法仿真切削破巖的可行性，為切削破巖研究提供了一種新方法。

切削破巖；SPH法；LS-DYNA；數(shù)值仿真；切削力；破碎角

0 引言

切削破巖在隧道挖掘、航道和港口建設(shè)中有廣泛應(yīng)用。但切削破巖中巖石破碎機(jī)理復(fù)雜，切削力計(jì)算困難。通過實(shí)驗(yàn)方法研究其破碎機(jī)理雖然可靠，但成本太高。目前國(guó)內(nèi)外越來越多學(xué)者通過數(shù)值仿真研究切削破巖。在切削破巖中使用的數(shù)值方法主要為有限元法和離散元法。Yu[1]和夏毅敏等[2]基于有限元法用LS-DYNA對(duì)切削破巖進(jìn)行數(shù)值仿真。Su等[3]，Huang等[4]和 Van 等[5]用離散元法模擬切削破巖，并求得切削力。上述研究仍無法準(zhǔn)確估計(jì)切削力和再現(xiàn)刀齒切削巖石破碎情況。近來SPH（Smoothed Particle Hydrodynamics）法在巖石破碎研究中有不少應(yīng)用，且在研究受壓作用下巖塊破碎機(jī)理和破碎形態(tài)[6-8]，沖擊破巖[9]，巖石坍塌[10]取得較好效果。Zhang等[11]和Li等[12]用SPH法仿真鎬齒切削巖石，探討了刀齒中應(yīng)力峰值和巖石變形、失效情況，但均未涉及切削力。本研究借助LS-DYNA用SPH法對(duì)刀齒切削破巖進(jìn)行數(shù)值仿真，重點(diǎn)分析該方法對(duì)切削力和破碎情況的預(yù)測(cè)能力，結(jié)合實(shí)驗(yàn)結(jié)果，探討SPH法用于切削破巖研究的可行性。

1 計(jì)算方法

SPH法在天體物理學(xué)問題研究中誕生，目前越來越多應(yīng)用于固體力學(xué)研究。該方法是一種無網(wǎng)格拉格朗日法，用于數(shù)值求解偏微分方程，其基礎(chǔ)是插值理論[6-7，10]。無網(wǎng)格的特點(diǎn)使其能夠處理許多不規(guī)則結(jié)構(gòu)和復(fù)雜變形問題。其基本原理是：將幾何體離散為許多具有一定體積和質(zhì)量的粒子，以粒子取代單元；每個(gè)粒子攜帶了必要的物理量（如：密度、速度、壓力和內(nèi)能等），任意時(shí)刻，空間任意一點(diǎn)的物理量通過對(duì)該時(shí)刻臨域內(nèi)所有粒子插值求和得到，從而得到連續(xù)場(chǎng)的近似；最后用離散化的場(chǎng)函數(shù)求解偏微分方程。SPH法引入了核函數(shù)對(duì)場(chǎng)函數(shù)進(jìn)行光滑，所得這種場(chǎng)的“近似”也叫場(chǎng)的核估計(jì)。場(chǎng)函數(shù)f在任意一點(diǎn)i的核估計(jì)＜f（xi）＞等于以點(diǎn)i為球心（圓心），2h為半徑的球形（圓心）區(qū)域內(nèi)所有粒子的場(chǎng)函數(shù)值插值的和（h為光滑長(zhǎng)度），即，

式中：xi為點(diǎn)i的位置坐標(biāo)；xj為粒子j的位置坐標(biāo)，均在空間固定坐標(biāo)系下；f（xj）為粒子j處場(chǎng)函數(shù)的值；W為核函數(shù)，反映點(diǎn)的距離與貢獻(xiàn)程度大小；h為光滑長(zhǎng)度，通常與點(diǎn)i距離超過2h的粒子對(duì)點(diǎn)i的物理量無貢獻(xiàn)，以xi為球心（圓心），2h為半徑的球形（圓形）區(qū)域也稱為點(diǎn)i的影響區(qū)域；j為落在點(diǎn)i影響區(qū)域內(nèi)粒子的編號(hào)；mj為粒子j的質(zhì)量；ρj為粒子j的密度。f在點(diǎn)i的微分形式的核估計(jì)為：

核函數(shù)W是光滑函數(shù)，可采用Monaghan[13]給出的三次樣條函數(shù)。

假定巖石為連續(xù)介質(zhì)，且在切削問題中不考慮重力，則質(zhì)量和動(dòng)量守恒方程的離散形式為：

式中：α、β表示坐標(biāo)方向；vα為速度分量；σαβ為應(yīng)力張量；∏ij為人工黏性，是為了使數(shù)值計(jì)算穩(wěn)定而人為加入的項(xiàng)；δαβ是克羅內(nèi)克函數(shù)，即α=β，δαβ=1，否則為0。為消除邊界處單元缺失帶來的計(jì)算誤差，需要對(duì)方程（4）采用重整化修正，詳見文獻(xiàn)[8]和[14]。

2 數(shù)值仿真模型

本研究對(duì)刀齒切削硬石膏[15]進(jìn)行數(shù)值仿真。原模型中巖樣尺寸為0.5 m×0.5 m×0.3 m，巖樣正方形面用膠水固定于切削平臺(tái)上，切削速度為0.15 m/s，切削厚度1.5～7.5 mm。將巖樣離散為間距1 mm的SPH單元，如果按照實(shí)際模型尺寸建模，并按實(shí)際速度切削，完成1 s切削，至少需要用1臺(tái)普通4核臺(tái)式計(jì)算機(jī)計(jì)算125萬h。受限于計(jì)算機(jī)內(nèi)存和計(jì)算能力，所取數(shù)值仿真模型尺寸比原模型小，如圖1，為100 mm×50 mm×30 mm，切削速度提高為10 m/s，仿真時(shí)間縮減為0.003 s，該模型下，計(jì)算耗時(shí)7.5 h。刀齒切削方向平行于巖樣上表面，且與巖樣端面垂直。在LS-PrePost中建立計(jì)算模型，巖樣單元間距為1 mm，刀齒劃分為用六面體和八面體立體實(shí)體單元，單元最大邊長(zhǎng)為2 mm。巖樣兩個(gè)側(cè)面和底面節(jié)點(diǎn)施加固定約束，限制6個(gè)自由度，其余面自由。刀齒做剛體平移，除沿著切削方向的平移運(yùn)動(dòng)外，其余自由度均被約束。

圖1 數(shù)值仿真模型Fig.1 Numerical simulation model

巖樣選用#173-Mohr-Coulomb材料模型，材料參數(shù)見表1[15]，其中巖樣剪脹角大小根據(jù)文獻(xiàn)[16]的研究，約為內(nèi)摩擦角的一半。因巖石動(dòng)力參數(shù)較難測(cè)得，因此本文所取均為準(zhǔn)靜態(tài)參數(shù)。刀齒選用#020-Rigid材料模型，密度為7 800 kg/m3，彈性模量為210 GPa，泊松比為0.3。

表1 巖樣材料參數(shù)Table 1 Rock material parameters

在LS-DYNA的*Control中設(shè)置體積黏性和沙漏能控制，均按默認(rèn)值設(shè)定。能量計(jì)算中將所有能量均包含在內(nèi)，包括沙漏能和耗散能。*Control_SPH中設(shè)置計(jì)算區(qū)域?yàn)檎麄€(gè)巖樣區(qū)域和巖樣正上方10 mm以內(nèi)的區(qū)域。仿真時(shí)間為0.003 s，每隔5×10-6s進(jìn)行1次采樣，輸出刀齒和巖樣的相互作用力，巖樣內(nèi)部應(yīng)力和應(yīng)變。

3 仿真結(jié)果與分析

將LS-PrePost生成計(jì)算K文件導(dǎo)入LS-DYNA進(jìn)行計(jì)算，再將結(jié)果導(dǎo)入LS-PrePost查看結(jié)果。得到的RCFORC文件可直接用MATLAB進(jìn)行處理和繪圖。

3.1 巖樣的應(yīng)力和塑性變形

用寬度為10 mm的刀齒切削巖樣，切削厚度為4.5 mm，切削前角為20°。圖2給出了數(shù)值仿真中t=0.001 5 s時(shí)刀齒齒尖附近巖樣的等效應(yīng)力和塑性變形，應(yīng)力單位為MPa，圖中未顯示刀齒。圖2（a）為過齒尖且垂直于切削方向所在截面的應(yīng)力和刀齒對(duì)稱面所在截面的應(yīng)力。從圖知，與刀齒刀刃接觸的巖樣單元應(yīng)力最大，尤其是與刀齒兩個(gè)齒尖接觸的單元，隨著與刀具距離的增加，巖樣單元受到的應(yīng)力值逐漸減小。隨著切削推進(jìn)，最大應(yīng)力區(qū)域也不斷向前移動(dòng)，但刀齒附近巖樣的應(yīng)力分布情況基本相同。在邊界附近的巖樣單元應(yīng)力值非常小，這說明巖樣模型尺寸變小后與原模型差異不大，可以近似視為相同。圖2（b）為與圖2（a）相同截面上的塑性變形。切削中，與刀齒接觸的巖樣塑性變形很大，且這些單元的相鄰單元的塑性變形也很大。刀齒向前切削，使其前方的巖樣單元逐漸失效并將其從巖樣上剝離，刀齒兩側(cè)和底部附近的巖樣單元均有較大塑性變形，巖樣被切削破壞。

圖2 t=0.001 5 s時(shí)刀齒附近巖樣的范式等效應(yīng)力和塑性變形Fig.2 Von Mises stresses and plastic displacements of rock elements near the cutter in different moments at t=0.001 5 s

3.2 巖樣破碎情況

刀齒切削過后，巖樣上會(huì)留下一條切削槽，切削槽各垂直于切削方向的截面形狀差異不大，且均有破碎角（breakout angle[15]）存在。圖3給出了不同剪脹角下，仿真所得距離起始端面距離為5 mm處切削槽截面形狀。從圖可知，剪脹角越大，切削槽破碎角越明顯（越大）。實(shí)際切削中，切削槽兩側(cè)均存在破碎角，但在此前的數(shù)值仿真均未提及。仿真中破碎角的出現(xiàn)能在現(xiàn)象上更接近實(shí)際，另外，它也改變了切削槽形狀，當(dāng)?shù)洱X再次切削經(jīng)過切削槽附近時(shí)，刀巖相互作用情況將與不存在破碎角時(shí)不同，進(jìn)而影響后續(xù)切削，這恰好是實(shí)際應(yīng)用中經(jīng)常碰到的問題，因此將會(huì)有重要意義。

圖3 不同剪脹角下切削槽形狀Fig.3 Shape of grooves at different dilation angles

3.3 切削力

切削中刀齒與巖石的相互作用力為切削力，可視為切向力、法向力和側(cè)向力的合力。由于仿真模型對(duì)稱，側(cè)向力理論上為零，實(shí)際仿真中側(cè)向力相對(duì)于切向力和法向力也很小，不考慮。圖4為切向力和法向力隨時(shí)間變化曲線。從刀齒與巖樣接觸開始，切削力迅速增加，而后在一定范圍內(nèi)波動(dòng)，與實(shí)際相符，但切削力波動(dòng)幅度要比各自平均值小很多，這與實(shí)際有差異。

為便于將仿真結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比，求出切削力最大值和平均值?？紤]到切削破巖中功耗主要受切向力影響，只對(duì)切向力進(jìn)行分析。按照表2中的切削參數(shù)調(diào)整仿真模型，其余參數(shù)保持不變，得到仿真結(jié)果，并求得最大切向力和平均切向力，見表2。從表可知，仿真所得切向力與實(shí)驗(yàn)結(jié)果比較接近，平均切向力差異更小。

圖4 切削力時(shí)歷曲線Fig.4 Cutting force vs time

表2 切削力實(shí)驗(yàn)結(jié)果與數(shù)值仿真結(jié)果對(duì)比Table 2 Comparison of the experimental results and numerical simulation results of cutting forces

用最小二乘法將仿真與實(shí)驗(yàn)中所得切向力進(jìn)行擬合，如圖5，圖中符號(hào)含義與表2相同。結(jié)果顯示，數(shù)值仿真結(jié)果與實(shí)驗(yàn)有較高相關(guān)性，尤其是平均切向力。數(shù)值仿真結(jié)果普遍比實(shí)驗(yàn)結(jié)果大，結(jié)合表2可知，切削前角為-10°和0°時(shí)，數(shù)值仿真所得最大切向力要比實(shí)驗(yàn)值大；切削前角為10°和20°時(shí)吻合得非常好；切削前角為30°時(shí)，仿真結(jié)果比實(shí)驗(yàn)值略小。平均切向力隨切削前角變化亦有類似規(guī)律，但不如前者明顯。這可能是切削前角較小時(shí)，巖樣在發(fā)生剪切失效之前已經(jīng)存在其他形式失效，因此切向力無需達(dá)到剪切失效所需值。

圖5 數(shù)值與實(shí)驗(yàn)切向力擬合Fig.5 Fitting of numerical and experimental tangential cutting forces

4 結(jié)語(yǔ)

1）借助LS-DYNA計(jì)算軟件，用SPH法求解偏微分方程，進(jìn)行刀齒切削破巖數(shù)值仿真，所得結(jié)果與實(shí)驗(yàn)結(jié)果較吻合。

2）在切削過程中，與刀齒接觸的巖樣單元受力和塑性應(yīng)變均最大，隨著離刀齒距離增加，應(yīng)力和塑性應(yīng)變均減小。最高應(yīng)力區(qū)始終出現(xiàn)在刀齒前刃面附近，隨著刀齒切削前進(jìn)而移動(dòng)。

3）切削過后，巖樣上形成的切削槽兩側(cè)出現(xiàn)破碎角，比此前的類似仿真更接近實(shí)際情況。仿真結(jié)果顯示，破碎角與剪脹角關(guān)系密切，剪脹角越大，破碎角越大。

4）數(shù)值仿真所得切削力時(shí)歷曲線與實(shí)際比較吻合。所得平均切向力大小與實(shí)驗(yàn)值接近。

[1]YU B,KHAIR A W.Numerical modeling of rock ridge breakage in rotary cutting[C]//Rock mechanics:meeting society's challenges and demands.Proceedings of the 1 st Canada-US rock mechanics symposium.Vancouver,Canada:American Rock Mechanics Association,2007:519-526.

[2] 夏毅敏，薛靜，周喜溫.軟巖切削過程的三維數(shù)值模擬[J].長(zhǎng)安大學(xué)學(xué)報(bào)：自然科學(xué)版，2010，30（4）：102-106.XIA Yi-min,XUE Jing,ZHOU Xi-wen.3-D dimensional numerical simulation of soft rock cutting process[J].Journal of Chang′an University:Natural Science Edition,2010,30(4):102-106.

[3] SU O,AKCIN N.Numerical simulation of rock cutting using the discreteelementmethod[J].International Journal of Rock Mechanics and Mining Sciences,2011,48(3):434-442.

[4]HUANG H,LECAMPION B,DETOURNAY E,et al.Discrete element modeling of tool-rock interaction I:rock cutting[J].International Journal for Numerical and Analytical Methods in Geomechanics,2013,37(13):1 913-1 929.

[5] VAN WYK G,ELS D N J,AKDOGAN G,et al.Discrete element simulation of tribological interactions in rock cutting[J].International Journal of Rock Mechanics and Mining Sciences,2014,65(1):8-19.

[6]DAS R,CLEARY P W.Effect of rock shapes on brittle fracture using Smoothed Particle Hydrodynamics[J].Theoretical and Applied Fracture Mechanics,2010,53(1):47-60.

[7]MA G W,WANG X J,REN F.Numerical simulation of compressive failure of heterogeneous rock-like materials using SPH method[J].International Journal of Rock Mechanics and Mining Sciences,2011,48(3):353-363.

[8] TANG X W,ZHOU Y D,LIU Y L.Factors influencing quasistatic modeling of deformation and failure in rock-like solids by the Smoothed Particle Hydrodynamics method[J].Mathematical Problems in Engineering,2013(1):681-703.

[9]SHAH Q H,HAMDANI A.The damage of unconfined granite edge due to the impact of varying stiffness projectiles[J].International Journal of Impact Engineering,2013,59(1):11-17.

[10]KAREKAL S,DAS R,MOSSE L,et al.Application of a mesh-free continuum method for simulation of rock caving processes[J].International Journal of Rock Mechanics and Mining Sciences,2011,48(5):703-711.

[11]ZHANG Q,LI S J.Simulation study of the strength of pick based on LS-DYNA in different hardness of rock[J].Advanced Materials Research,2012,567(1):232-235.

[12]ZHAO L J,ZHOU Z H,GUAN Q Z,et al.Application of SPH in numerical simulation of roadheader hard rock cutting[J].Advanced Materials Research,2013,619(1):203-206.

[13]MONAGHAN J J,GINGOLD R A.Shock simulation by the particle method SPH[J].Journal of Computational Physics,1983,52(2):374-389.

[14]RANDLESPW,LIBERSKYLD.SmoothedParticleHydrodynamics:some recent improvements and applications[J].Computer Methods in Applied Mechanics and Engineering,1996,139(1-4):375.

[15]NUH Bilgin.Investigations into the mechanical cutting characteristics of some medium and high strength rocks[D].Newcastle Upon Tyne:University of Newcastle Upon Tyne,1977.

[16]ZHAO X G,CAI M.Influence of plastic shear strain and confinement-dependent rock dilation on rock failure and displacement near an excavation boundary[J].International Journal of Rock Mechanics and Mining Sciences,2010,47(5):723-738.

Numerical simulation of rock cutting using SPH method

OUYANG Yi-ping1,YANG Qi1,2*,CHENG Xu-dong3
（1.State Key Laboratory of Ocean Engineering,Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration,Shanghai Jiao Tong University,Shanghai 200240,China;2.Shanghai Jiao Tong University Underwater Engineering Institute Co.,Ltd.,Shanghai 200231,China;3.CRCC Harbour&Channel Engineering Bureau Group Co.,Ltd.,Zhuhai,Guangdong 519070,China）

Rock cutting of different cutting parameters were simulated by using SPH method to solve the partial differential equation with the help of LS-DYNA software,the rock cutting experiment with different cutting parameters is carried out as a prototype,and the cutting force is compared with the experimental results.The numerical results were in accord with the experimental results.The Von Mises stresses and plastic strains of those elements which contacted with the cutter were maximum.The cutting groove was formed on the rock sample after cutting,and there are breakout angles on both sides,which was much similar to the actual situation than former simulation results,and the size of the breakout angle is closely related to the dilatancy angle.A series simulation results show that the change law and the numerical value of the cutting force is in accord with the experiments,especially the mean tangential cutting forces.It is proved to be a suitable method that using SPH method to simulate rock cutting,and provides a new method for rock cutting research.

rock cutting;SPH mothed;LS-DYNA;numerical simulation;cutting force;breakout angle

U615.31；TU45；V448.253

A

2095-7874（2017）10-0013-05

10.7640/zggwjs201710003

2017-03-17

2017-05-23

國(guó)家自然科學(xué)基金項(xiàng)目（51179104）

歐陽(yáng)義平（1987— ），男，湖南永州人，博士在讀，主要從事船舶與海洋和疏浚設(shè)備研究。

*通訊作者：楊啟，E-mail：yangqi110@sjtu.edu.cn