蘇翠俠，王燕群，蔡宗熙，亢一瀾，黃 田

(天津大學(xué)機械工程學(xué)院，天津 300072)

盾構(gòu)掘進機是目前隧道掘進的專門工程機械，被廣泛用于地下交通、運輸管道等隧道工程建設(shè)．施工過程中刀盤系統(tǒng)與土體相互耦合作用，并受到地質(zhì)等因素的影響，刀盤載荷十分復(fù)雜，因此，研究刀盤載荷分布規(guī)律，分析刀盤系統(tǒng)與土體的適應(yīng)性，對優(yōu)化刀盤拓撲結(jié)構(gòu)、提高施工效率具有指導(dǎo)意義．

盾構(gòu)刀盤的作業(yè)空間處于地表以下，屬于封閉式作業(yè)且工作環(huán)境復(fù)雜多變，無法直接觀測刀盤在掘進過程中的載荷分布，因此目前針對刀盤載荷的研究較少．以往對盾構(gòu)刀盤的研究主要集中在簡化理論與模型實驗等方面：簡化理論方面主要是通過分析各種刀具的切削機理[1]推導(dǎo)盾構(gòu)掘進所需總載荷[2]，進而分析刀盤一般設(shè)計理論及其影響因素[3]；模型實驗方面主要是通過制造等比例小型簡化刀盤，進行土箱實驗分析不同地質(zhì)條件和掘進參數(shù)對掘進推力、扭矩等載荷的影響[4]；另外一種方法是結(jié)合大量的實際工程的施工數(shù)據(jù)建立刀盤掘進參數(shù)的經(jīng)驗公式[5]，總結(jié)各種刀具破壞形式及其檢測方法[6]，對影響刀盤掘進效率的因素進行分析及改進．

近年來隨著數(shù)值技術(shù)的發(fā)展，有限元法為研究盾構(gòu)刀盤載荷提供了一種新工具．盾構(gòu)刀盤的載荷主要由刀盤與土體相互作用產(chǎn)生，刀盤掘進過程本質(zhì)上可以視為動態(tài)的切削問題．采用有限元法從切削角度分析刀盤掘進必須考慮以下幾個方面的問題：有限元模型的建立、土體本構(gòu)的選擇和切屑分離的模擬．

創(chuàng)建適合的模型是數(shù)值仿真的首要問題，在滿足計算精度要求的基礎(chǔ)上應(yīng)對模型進行合理的簡化以節(jié)約計算成本．

土體本構(gòu)的選擇直接影響計算結(jié)果，近年來土力學(xué)的發(fā)展較為成熟，已經(jīng)有越來越多的適合描述不同巖土特性的本構(gòu)模型被提出，目前應(yīng)用較多的巖土模型主要有Duncan-Chang 雙曲線模型、Drucker-Prager非線性彈塑性模型、Cam-Clay 模型、黏彈塑性蓋帽模型等，研究結(jié)果表明它們大都能夠給出滿足適當精度要求的結(jié)果．

有限元模擬切屑分離的方法有很多，根據(jù)分離依據(jù)可以分為兩類：基于應(yīng)力、應(yīng)變能密度或有效塑性應(yīng)變等[7]物理量的物理準則和基于距離等幾何尺寸的幾何準則[8]．具體實現(xiàn)方法包括單元移除、節(jié)點脫粘、節(jié)點分離等．切削仿真中上述切削分離方法選取不當經(jīng)常引起網(wǎng)格畸變導(dǎo)致計算不收斂，采用離散元和有限元結(jié)合[9]、ALE 方法[10]、流固耦合、自適應(yīng)網(wǎng)格劃分技術(shù)[11]等方法在一定程度上能夠描述網(wǎng)格大變形問題，但在應(yīng)用中存在各自的局限性．

以往的研究或著眼于靜力分析或?qū)⒌侗P簡化為簡單的柱形剛體，都忽略了切削過程，相關(guān)的耕耘機刀具與土壤的切削仿真[12]以及金屬切削仿真[8]研究為有限元分析盾構(gòu)刀盤切削問題提供了一定的參考．Shen 等[10]在此基礎(chǔ)上將有限元分析土壤切削的方法應(yīng)用于盾構(gòu)刀盤掘進分析，針對盾構(gòu)刀盤建立了一個簡化三維模型，實現(xiàn)了刀盤掘進仿真．Shen 采用的Ls-Dyna 中的MAT147 材料模型對標準的摩爾-庫倫準則進行了修正，確定了與壓力相關(guān)的峰值剪切強度，并提出了將ALE 方法應(yīng)用于刀盤-土體切削，從流體-固體耦合的角度處理了刀盤-土體的相互作用，解決了切削引起的土體大變形問題．但是一般來說，采用ALE、流固耦合等方法計算時間較長，占用了大量的計算資源；且此類方法研究的刀盤模型都比較簡單，給出的關(guān)于刀盤載荷的信息不是特別豐富．正是由于這些因素的存在，至今未能得到刀盤在掘進過程中的更多的有效的載荷信息．

針對盾構(gòu)刀盤掘進載荷問題，筆者采用有限元軟件ABAQUS/Explicit 進行了數(shù)值模擬仿真計算，建立了模擬盾構(gòu)刀盤掘進過程中連續(xù)切削土體的三維動態(tài)仿真模型．在計算中采用了擴展的 Ducker-Prager 非線性彈塑性的材料本構(gòu)關(guān)系，應(yīng)用了具有單元刪除功能的損傷失效準則，實現(xiàn)了對盾構(gòu)刀盤掘進過程中的土體切屑分離的數(shù)值模擬．最后，本文基于三維仿真結(jié)果，討論了盾構(gòu)刀盤動態(tài)掘進過程中的切削載荷沿刀盤半徑分布、隨時間變化規(guī)律以及刀盤系統(tǒng)切削土體所受到的扭矩．

1 三維動態(tài)仿真模型

隧道掘進是一個復(fù)雜的施工過程，它對周圍環(huán)境擾動較大，且由于盾構(gòu)裝備結(jié)構(gòu)的復(fù)雜性，機械本身和施工環(huán)境影響盾構(gòu)掘進過程的因素也很多．本文重點研究盾構(gòu)刀盤掘進過程中的載荷分布，因此有限元建模包括了盾構(gòu)機刀盤與正前方相互接觸區(qū)域的土體，以此作為盾構(gòu)刀盤系統(tǒng)三維動態(tài)仿真整體模型．

1.1 模型介紹

以天津地鐵工程中的土壓平衡盾構(gòu)刀盤作為原型．刀盤直徑為6.4,m，根據(jù)刀盤的二維設(shè)計圖紙，經(jīng)適當簡化建立相應(yīng)的三維有限元模型．簡化包括刪除不影響結(jié)構(gòu)強度的螺栓孔、泡沫注入口、倒角、拐角等特征．刀盤原型以切刀為主，由于面板相鄰切刀的刀間距較小，因此對切刀進行簡化，在保證刀具前角及后角不變的基礎(chǔ)上在切刀布刀位置以連續(xù)刀刃代替分散的切刀，如圖1 所示．

在刀盤正前方加入圓柱形的土體模型，形成刀盤切削土體的整體有限元分析模型，如圖2 所示．整個有限元模型采用六面體8 節(jié)點實體單元劃分網(wǎng)格，模型節(jié)點總數(shù)為104,364，單元總數(shù)為80,828．為得到較好的計算精度并提高計算效率，將仿真過程中與刀盤直接接觸的被切削的部分土體網(wǎng)格加密，未與刀盤直接作用的土體網(wǎng)格適當加大．

盾構(gòu)刀盤的材料為Q345 鋼．彈性模量E= 210 GPa ，密度 ρ=7 850 kg m3，泊松比ν= 0.26．土體取天津地區(qū)的黏性軟土，主要材料參數(shù)為：彈性模量E= 7.098 MPa ，密度ρ=1 850 kg m3，泊松比ν= 0.3，內(nèi)聚力9.5 kPa ，摩擦角23°．

圖1 盾構(gòu)刀盤模型Fig.1 Model of cutterhead

圖2 刀盤與土體模型Fig.2 Model of cutterhead and soil

在實際施工中，盾構(gòu)刀盤旋轉(zhuǎn)切割土體的同時受到后方液壓千斤頂?shù)捻斶M作用，以緩慢速度向前推進．鑒于本文中只是初步模擬盾構(gòu)刀盤切削土體過程，綜合以上各點，對實際分析模型進行相應(yīng)考慮和簡化，施加的載荷和位移邊界條件可總結(jié)如下．

(1) 初始狀態(tài)盾構(gòu)刀盤與土體即將接觸．

(2) 刀盤的旋轉(zhuǎn)速度為1 r min ．

(3) 刀盤的推進速度為40 mm min ．

(4) 仿真時間為30,s，即盾構(gòu)刀盤環(huán)向切割土體1/2 圈．轉(zhuǎn)速和推進速度在前5,s 以光滑加載方式由零增加到最終值，隨后保持不變．

(5) 約束土體模型外邊界的位移自由度，保持待開挖表面為自由表面．

(6) 利用非光滑接觸條件模擬盾構(gòu)刀盤和土壤之間的相互作用．

1.2 土體模型

掘進過程中土體材料特性直接影響刀盤載荷分布．盾構(gòu)掘進施工面臨的地質(zhì)情況復(fù)雜多變，應(yīng)根據(jù)實際工程采用與之相適應(yīng)的本構(gòu)模型．本文采用了文獻[13]的擴展Drucker-Prager 非線性彈塑性本構(gòu)．

與常用的Tresca 屈服條件、von Mises 屈服條件和雙剪應(yīng)力屈服條件相比，Drucker-Prager 模型不僅在屈服準則中引入了中間主應(yīng)力對材料屈服面的影響，同時考慮了靜水壓力對材料屈服面的影響.Mohr-Coulomb 屈服條件雖然應(yīng)用比較廣泛，但它的屈服面在Π平面上表現(xiàn)為一個不等角的等邊六邊形，具有棱角奇異性，應(yīng)用于數(shù)值計算時會帶來較大的困難．Drucker-Prager 屈服準則克服了這一缺點，并且選擇適當?shù)牟牧铣?shù)可以與 Coulomb 模型相匹配．本文采用的土體模型是對Drucker-Prager 準則加以改進的擴展Drucker-Prager 準則，即

式中：r 為偏應(yīng)力第三應(yīng)力不變量；K為單軸拉伸屈服應(yīng)力與單軸壓縮屈服應(yīng)力比值；d為凝聚力；β為材料摩擦角；q 為Mises 等效應(yīng)力；p為平均壓應(yīng)力；σ1、σ2和σ3分別為3 個方向的主應(yīng)力．實際計算時要考慮土體圍壓的影響．

在傳統(tǒng)的Drucker-Prager 準則基礎(chǔ)上，考慮拉伸、壓縮環(huán)境不同對屈服面的影響引入?yún)?shù)K，如圖3 所示，當材料的拉伸強度與壓縮強度相同時K= 1,t=q，即退化為傳統(tǒng)的Drucker-Prager 準則．流動法則采用關(guān)聯(lián)流動法則．

圖3 偏平面內(nèi)線性擴展Drucker-Prager 屈服準則Fig.3 Typical yield surfaces of Drucker-Prager model in the deviatoric plane

1.3 包含單元刪除功能的單元損傷失效模型

土體從屈服到破壞是一個連續(xù)的過程，以往的仿真研究中對土體切屑的處理一般認為屈服即破壞，即以土體屈服點的各種特征判斷切屑分離、形成準則，而忽略了切削形成到分離這一過程．引用的包含單元刪除功能的單元損傷失效模型則描述了材料從屈服發(fā)展到破壞響應(yīng)過程．

單元刪除功能是為了克服有限元本身缺陷的一種方法．有限元是基于連續(xù)介質(zhì)力學(xué)原理的，即物質(zhì)域在空間中連續(xù)，因此單元本身是不會消失的．而在實際情況中，損傷斷裂的存在必會使得一些單元消失或完全的失效，為了能夠模擬這種情況，有限元中引入了單元刪除功能．單元損傷失效是基于斷裂力學(xué)描述損傷對于材料破壞的影響而提出的，假設(shè)基于特定本構(gòu)關(guān)系的單元材料在達到強度極限以后，材料剛度按照一定的規(guī)律逐漸衰減到零，單元完全喪失承載能力并退出整體模型的計算．如圖4 所示，單元損傷失效過程包含3 個部分：單元失效前的材料響應(yīng)AB段、初始破壞點B點(由初始損傷準則判定)和損傷演變規(guī)律BC段．

圖4 損傷失效模型應(yīng)力-應(yīng)變響應(yīng)曲線Fig.4 Stress-strain curve with progressive damage degradation

如圖4 所示，AB段材料處于彈塑性變形階段，材料達到強度極限后，由單元積分點的等效塑性應(yīng)變建立基于剪切失效準則的初始損傷準則，即假設(shè)材料在開始破壞之前的塑性應(yīng)變?yōu)椋x一個描述塑性變形隨等效塑性應(yīng)變遞增的狀態(tài)變量為

式中：θs為剪應(yīng)力率；τmax為最大剪應(yīng)力；為應(yīng)變率．

當ω s達到1 時，達到初始破壞點B．

材料達到初始破壞點后剛度開始衰減直至喪失承載能力，損傷演變規(guī)律描述了材料剛度衰減規(guī)律(圖4 中BC段)．引入損傷變量D(當D= 1時，材料完全失去承載能力)，初始損傷發(fā)生后任意時刻材料的應(yīng)力張量

式中為當前時刻不考慮失效時(圖4 中BC′段)的有效應(yīng)力張量．

本文采用位移控制損傷定義材料損傷演變規(guī)律．如上所述，當材料達到剪切破壞準則以后，有效塑性位移由式(3)所示的損傷演變方程決定．假設(shè)材料完全失效時(如圖4 中C點)的有效塑性位移為，定義

式中L 為單元特征長度．

單元完全失效后(圖4 中,C點)，隨即被從模型中刪除不再顯示，就像失效單元已經(jīng)消失了一樣．基于連續(xù)介質(zhì)力學(xué)原理，有限元實際處理單元刪除問題時，材料完全失效時單元并未消失，只是失效單元剛度乘了一個極小的數(shù)值后接近于零，使得失效單元退出對有限元整體模型的貢獻．

2 盾構(gòu)刀盤掘進三維仿真結(jié)果分析

2.1 刀盤刀刃切向切削力分布規(guī)律

由刀盤結(jié)構(gòu)可以看出刀具始終高于刀盤面板．在掘進過程中，刀具首先與掌子面接觸并成為整個切削過程的主體，也是整個刀盤切削載荷的主要承力部位．刀具與土體相互作用使得刀具受到垂直于前刃面的作用力，將此作用力按照刀具的切削軌跡沿切向及法向分解即得到切向破壞土體所需要的切削力和法向與土體擠壓的作用力．以下就盾構(gòu)刀盤刀刃部分在掘進過程的數(shù)值計算結(jié)果進行分析．

圖5所示為R= 2.8 m 處刀刃切向切削力在掘進過程中隨掘進時間t的變化曲線，其中R為刀刃所在位置的半徑．從圖中可以看出，在掘進初始階段，刀盤刀刃與土體逐漸接觸，切削力迅速增大，經(jīng)過一段時間進入穩(wěn)定掘進狀態(tài)，切削力也趨于穩(wěn)定．分析刀具切削路徑，由于沿徑向分布的刀刃各部分切削形式相同、切削路徑一致，刀刃各部分的切向切削力隨時間變化規(guī)律一致．

圖5 R＝2.8 m 處刀刃切向切削力隨時間的變化Fig.5 Change of tangential force of simplified cutter when R＝2.8 m

圖6 所示為掘進過程中 3.6 st= 時盾構(gòu)刀盤刀刃切向切削力沿半徑方向的分布．可以看出：隨著半徑增大，刀刃切向切削力逐漸增大．此時為掘進初始階段，刀盤緩慢加載，刀盤刀刃與土體逐漸接觸，土體單元還未出現(xiàn)失效．由于距刀盤中心的距離不同，刀刃各部分的切削進程也不同，距離刀盤中心較遠的刀刃切削力更快地達到穩(wěn)定值．此時盾構(gòu)刀盤刀刃切削力與布刀位置有關(guān)．

圖7給出了掘進過程中 25 st= 時盾構(gòu)刀盤刀刃切向切削力沿半徑方向的分布．可以看出，刀刃切向切削力沿半徑方向分布變化較?。蜻M過程達到穩(wěn)定階段后，刀刃各部分切削形式相同、切削路徑一致，切削力趨于穩(wěn)定，切向切削力沿半徑變化較為平緩；且隨時間變化較小．距離刀盤中心較近的刀刃由于受中心刀的影響，同時切向線速度較低，切削力相對較小．刀盤面板在半徑為2.0,m 處含有環(huán)形筋板結(jié)構(gòu)，由于筋板的影響，刀刃刀角大于實際切刀刀角，該處切向力較?。?/p>

圖6 t＝3.6 s 時刀刃切向切削力沿半徑方向的分布Fig.6 Tangential force distribution of simplified cutter in the radial direction when t＝3.6 s

圖7 t＝25 s 時刀刃切向切削力沿半徑方向的分布Fig.7 Tangential force distribution of simplified cutter in the radial direction when t＝25 s

2.2 刀盤刀刃法向擠壓力分布規(guī)律

圖8給出了掘進過程中刀盤刀刃法向擠壓力合力隨時間變化曲線．可以看出，刀刃法向合力在刀盤與土接觸后迅速增大到穩(wěn)定值．刀盤轉(zhuǎn)動1/4 轉(zhuǎn)后首層單元失效完成1 個切削周期，18 s 時刀刃與重新形成的掌子面接觸開始第2 個切削周期．刀盤刀刃法向合力具有周期穩(wěn)定性．

圖8 掘進過程刀盤刀刃切削力法向合力隨時間的變化Fig.8 Change of total normal force of simplified cutter during excavation

圖9給出了掘進初始階段 3.6 st= 時刀刃法向力沿半徑的分布曲線．可以看出：隨著半徑增大，刀刃法向力逐漸增大．由刀盤結(jié)構(gòu)可以看出掘進初始階段，刀盤面板與掌子面之間存在空隙，凸出面板許多的中心刀將首先與土體接觸并成為面板與掌子面相互作用法向力的主要承重部位，受中心刀的影響，刀盤與掌子面發(fā)生變形．隨著刀盤繼續(xù)推進距中心刀較遠的刀刃先于距中心刀較近的刀刃接觸土體，導(dǎo)致同一時刻刀刃法向力隨半徑增大．

圖9 t＝3.6 s 時刀刃法向力沿半徑的分布Fig.9 Normal force distribution of simplified cutter in the radial direction when t＝3.6 s

圖10給出了掘進過程中 25 st= 時刀刃法向力沿半徑的分布曲線．可以看出：刀刃法向力沿半徑方向分布變化不大．與刀刃切向力同理，掘進穩(wěn)定后刀刃各部分切削環(huán)境相似、切削形式相同，法向切削力在掘進穩(wěn)定后與布刀位置無關(guān)．由于筋板處刀刃刀角增大，導(dǎo)致法向力在半徑為2.0,m 處數(shù)值略有減?。?/p>

圖10 t＝25 s 時刀刃法向力沿半徑的分布Fig.10 Normal force distribution of simplified cutter in the radial direction when t＝25 s

2.3 刀盤掘進過程扭矩變化

將刀盤上與土體相互作用產(chǎn)生的力沿面板面內(nèi)對刀盤中心取矩即得到刀盤掘進過程中克服土體阻力所需的扭矩．圖11 給出了有限元掘進仿真過程中整個盾構(gòu)刀盤受到的扭矩隨時間的變化情況．可以看出：扭矩在掘進初始時刻隨時間迅速增大到穩(wěn)定值，隨后保持不變．盾構(gòu)實際施工過程中，當掘進過程穩(wěn)定后，切削狀態(tài)達到平衡，各掘進參量在整個掘進過程中保持穩(wěn)定不變．圖11 給出的扭矩隨時間的變化與實際施工情況相符，反映了刀盤扭矩的變化規(guī)律．

圖11 掘進過程中刀盤扭矩隨時間的變化Fig.11 Change of torque of cutterhead during excavation

2.4 刀盤面板應(yīng)力分布

圖12給出了掘進過程穩(wěn)定后(25 st = 時)刀盤面板應(yīng)力分布．可以看出：刀盤整體受力均勻，高應(yīng)力區(qū)主要集中在刀刃和開口處以及刀盤筋板與面板連接處.最大應(yīng)力為56,MPa，出現(xiàn)在筋板與面板連接處．計算結(jié)果表明，隨著刀盤掘進仿真趨于穩(wěn)定，刀盤自身的受力情況也逐漸穩(wěn)定，刀盤面板大部分區(qū)域應(yīng)力水平低于20,MPa，遠遠小于強度極限．

以往有限元研究刀盤載荷問題時，通常是直接將均布載荷施加在刀盤面板上，將動態(tài)問題簡化為靜力分析，忽略了刀盤與土體的相互作用．本文模型較以往的相關(guān)研究的有限元簡化模型更接近真實刀盤，并考慮了刀盤切削土體的動態(tài)效應(yīng)，因此本文中給出的盾構(gòu)刀盤在掘進仿真過程中的載荷分布與實際情況更加接近．

圖12 t＝25,s時刀盤有效應(yīng)力分布Fig.12 Nephogram of stress pattern of cutterhead when t＝25,s

3 結(jié) 論

(1) 實現(xiàn)了盾構(gòu)刀盤掘進過程的直接數(shù)值模擬．利用包含單元刪除功能的損傷失效模型模擬土體從屈服到破壞這一連續(xù)性過程，不僅從形態(tài)上反映了盾構(gòu)刀盤切削土體的物理現(xiàn)象，且單元失效后退出計算時的變形相對較小，在一定程度上避免了網(wǎng)格出現(xiàn)極大的扭曲和畸變．與ALE、流固耦合等技術(shù)相比，本文方法求解方程規(guī)模小，節(jié)省了計算資源．

(2) 得到了刀盤切削載荷隨時間變化規(guī)律及沿徑向分布規(guī)律．掘進過程中，盾構(gòu)刀盤的切削載荷主要集中在刀具上．在掘進初始階段，由于中心刀與刀刃切削進程的影響，刀刃法向擠壓力和切向切削力沿徑向遞增且變化較大；而在掘進穩(wěn)定后，則沿半徑變化較?。?/p>

(3) 得到了刀盤切削扭矩隨時間變化規(guī)律．在掘進初始階段，盾構(gòu)刀盤緩慢加載，刀盤與土體逐漸接觸，刀盤切削力矩逐漸增加．掘進穩(wěn)定后，刀盤的刀刃切削及受力狀態(tài)也趨于穩(wěn)定，扭矩增大到最大值后保持穩(wěn)定不變．總體變化趨勢與實際相符．

［1］宋克志，潘愛國. 盾構(gòu)切削刀具的工作原理分析[J].建筑機械，2007(3)：74-76.Song Kezhi，Pan Aiguo. Operation principle analysis of cutting tools on shield[J].Construction Machinery，2007(3)：74-76(in Chinese).

［2］管會生. 盾構(gòu)刀盤扭矩估算的理論模型[J]. 西南交通大學(xué)學(xué)報，2008，43(2)：213-217.Guan Huisheng. Theoretical model for estimation of cutterhead torque in shield tunneling[J].Journal of Southwest Jiaotong University，2008，43(2)：213-217(in Chinese).

［3］[日]土木學(xué)會. 隧道標準規(guī)范(盾構(gòu)篇)及解說[M].宋 偉，譯. 北京：中國建筑工業(yè)出版社，2001.Japan Society of Civil Engineers.Japanese Standard for Shield Tunneling[M]. Song Wei，Trans. Beijing：China Architectural Industry Press，2001(in Chinese).

［4］朱合華，徐前衛(wèi)，鄭七振，等. 軟土地層土壓平衡盾構(gòu)施工參數(shù)的模型試驗研究[J]. 工木工程學(xué)報，2007，40(9)：87-94.Zhu Hehua，Xu Qianwei，Zheng Qizheng，et al. Experimental study on the working parameters of EPB shield tunneling in soft ground[J].China Civil Engineering Journal，2007，40(9)：87-94(in Chinese).

［5］Bernhard Maidl，Martin Herrenknecht，Lothar Nheuser.Mechanized Shield Tunneling[M]. Berlin：Ernst and Sohn，1996.

［6］趙維剛，劉明月，杜彥良，等. 全斷面隧道掘進機刀具異常磨損的識別分析[J]. 中國機械工程，2007，18(2)：150-153.Zhao Weigang，Liu Mingyue，Du Yanliang，et al. Abnormal cutter wear recognition of full face tunnel boring machine(TBM)[J].China Mechanical Engineering，2007，18(2)：150-153(in Chinese).

［7］Carroll J T，Strenkowski J S. Finite element models of orthogonal cutting with application to single point dia-mond turning [J].International Journal of Mechanical Sciences，1988(11)：899-920.

［8］王 霄，盧樹斌，高傳玉. 金屬切削加工有限元模擬技術(shù)的研究[J]. 煤礦機械，2006，27(11)：51-54.Wang Xiao，Lu Shubin，Gao Chuanyu. Research on finite element modeling techniques of metal cutting[J].Coal Mine Machinery，2006，27(11)：51-54(in Chinese).

［9］Onate E，Rojek J. Combination of discrete element and finite element methods for dynamic analysis of geomechanics problems[J].Computer Methods in Applied Mechanics and Engineering，2004，193(27/28/29)：3087-3128.

［10］Shen J Q，Jin X L，Li Y，et al. Numerical simulation of cutterhead and soil interaction in slurry shield tunneling[J].Engineering Computations，2009，26(7/8)：985-1005.

［11］Susila E，Hryciw R D. Large displacement FEM modeling of the cone penetration test(CPT)in normally consolidated sand[J].International Journal for Numerical and Analytical Methods in Geomechanics， 2003 ，27(7)：585-602.

［12］Mootaz Abo-Elnor，Hamilton R，Boyleb J T. 3D dynamic analysis of soil-tool interaction using the finite element method[J].Journal of Terramechanics，2003，40(1)：51-62.

［13］Hibbitt，Karlsson and Sorensen，Inc.ABAQUS Analysis User’s Manual Help Online[M]. USA：SIMULIA，2005.