基于剛體動力學(xué)-離散元耦合數(shù)值方法的復(fù)合地層盾構(gòu)滾刀響應(yīng)

徐汪豪， 姚志剛， 方 勇， 王飛陽， 周 健

(1.浙江杭海城際鐵路有限公司， 海寧 314499； 2.西南交通大學(xué)交通隧道工程教育部重點實驗室， 成都 610031； 3.中鐵第四勘察設(shè)計院集團有限公司， 武漢 430063)

近年來，中國隧道及地下工程開始大規(guī)模建設(shè)，各領(lǐng)域的隧道總數(shù)及長度均快速增加[1]。在這些隧道中，較多的工程使用了全斷面機械化掘進設(shè)備，如盾構(gòu)機、TBM設(shè)備。但是，隨著中國盾構(gòu)機應(yīng)用領(lǐng)域及應(yīng)用數(shù)量的不斷擴大，機械設(shè)備施工面臨的困難與日俱增[2]。同時，隨著盾構(gòu)機面臨愈加復(fù)雜的地層環(huán)境，刀盤多同時采用滾刀和刮刀[3]。

在復(fù)合地層中，采用滾刀刀具的盾構(gòu)機將面臨偏磨的挑戰(zhàn)[4]，偏磨的產(chǎn)生與隧道地質(zhì)情況、滾刀響應(yīng)機理有關(guān)[5-11]。目前室內(nèi)實驗?zāi)茌^好得分析滾刀對巖土體的破壞作用。鄒飛等[12]基于實驗分析了滾刀作用下巖體位移與應(yīng)變場。張桂菊等[13]采用數(shù)值分析和實驗相結(jié)合的形式對滾刀切削力學(xué)進行了分析。由于滾刀運動狀態(tài)在實際掘進過程中的不可見性，研究者多采用數(shù)值模擬的手段進行研究。van Wyk等[14]以Paarl花崗巖為宏觀對象，通過離散元數(shù)值分析得到了刀具的切削力和切削深度以及刀具磨損后的刀具切削面大小之間的聯(lián)系。孫偉等[15]使用三維離散元方法，發(fā)現(xiàn)刀刃寬的增大將導(dǎo)致滾刀受力和碎巖體積的增加，而刀刃角對這兩者的影響較?。粡埧萚16]基于二維離散元數(shù)值方法，認為滾刀破巖過程中存在4種基本的破巖模式，這4種破巖模式和圍巖以及刀間距有關(guān)。Cho等[17]使用三維有限元方法，將數(shù)值分析中滾刀的受力行為和室內(nèi)滾刀線性切割試驗進行了對比，最終利用三維有限元方法優(yōu)化TBM(tunnel boring machine)滾刀間距。

目前以數(shù)值模擬為手段對滾刀的研究中，依然存在一定不足。前人研究中，針對硬巖破巖效率的研究較多，針對復(fù)合地層滾刀運動規(guī)律的較少；在三維數(shù)值分析中，多將滾刀轉(zhuǎn)動速度設(shè)定為恒定值[14-20]，認為滾刀的角速度與線速度存在理想匹配，無法實現(xiàn)滾刀轉(zhuǎn)動的動態(tài)平衡，與實際工況存在差異性。為此，使用剛體動力學(xué)(rigid body dynamics，RBD)-離散元(discrete element method,DEM)耦合數(shù)值模擬方法實現(xiàn)滾刀的動態(tài)滾動行為，并以此進行分析。

1 滾刀偏磨與滾刀響應(yīng)機理分析

隨著盾構(gòu)機制造技術(shù)的不斷進步，大直徑盾構(gòu)機的廣泛應(yīng)用，刀盤內(nèi)外側(cè)滾刀的運動速度存在顯著差異，如圖1所示。

滾刀的自轉(zhuǎn)狀態(tài)是由地層摩擦受力的綜合體現(xiàn)，如圖2所示。在圖2中，滾刀迎土面受到的摩擦力f、破巖力F等，綜合體現(xiàn)為滾刀軸承處的法向力Fn和切向力Ft，以及維持滾刀轉(zhuǎn)動的扭矩Tt和阻礙滾刀轉(zhuǎn)動的阻礙扭矩Tci，在大多數(shù)的理想巖體中，Tt一般大于Tci或維持動態(tài)平衡，滾刀能夠形成和線速度相匹配的角速度，導(dǎo)致滾刀接觸面和巖面的相對速度較低。

但是在復(fù)合地層中滾刀面對離散性較大的軟質(zhì)巖土體，滾刀轉(zhuǎn)動扭矩Tt和阻礙扭矩Tci并不能維持較好的動態(tài)平衡狀態(tài)。如滾刀啟動時無法提供足夠的Tt，將導(dǎo)致滾刀角速度和線速度動態(tài)匹配差異性較大，使得滾刀在巖土體的界面處存在較大的相對速度，如圖3所示。

滾刀與巖土體界面處的相對速度大小將決定滾刀刀圈的磨損速率，同時進一步影響滾刀的運動狀態(tài)即角速度和線速度的動態(tài)匹配效果。在單一勻質(zhì)巖體中，滾刀刀圈接觸面處的實際線速度和理論線速度較為接近，可以處于較為理想的工作狀態(tài)；在復(fù)合地層中，滾刀不能達到較好的動態(tài)匹配效果，接觸面處相對速度較高，從而導(dǎo)致滾刀處于非正常工作狀態(tài)。現(xiàn)有工程研究也已表明，當滾刀無法保證正常的自轉(zhuǎn)運動，將導(dǎo)致刀具的異常磨損[5-11]，如圖4所示。

圖1 不同安裝位置下的滾刀線速度方向及大小Fig.1 The direction and size of the hob linear velocity under different installation positions

∝表示角速度與線速度的匹配圖2 滾刀的運動狀態(tài)Fig.2 Movement state of disc cutter

圖3 復(fù)合地層中滾刀界面速度的差異性Fig.3 Difference of interface velocity in composite stratum

圖4 復(fù)合地層中的滾刀偏磨Fig.4 Eccentric wear of disc cutter in composite stratum

2 RBD-DEM耦合數(shù)值模擬方法

常見DEM數(shù)值模擬中，幾何體對象使用較為固定的運動設(shè)定，而實際工況中的幾何體存在較為復(fù)雜的動力學(xué)運動RBD，如被動轉(zhuǎn)動、被動移動等。因此，使用RBD-DEM耦合數(shù)值模擬的方法[21]，將DEM中剛體受到的復(fù)雜受力解傳遞至RBD中進行計算以決定DEM中下一時刻剛體的運動狀態(tài)，最終實現(xiàn)滾刀轉(zhuǎn)動被動平衡的數(shù)值模擬。

圖5 剛體動力學(xué)計算時的物體坐標變化Fig.5 Object coordinate change in RBD numerical simulation

剛體動力學(xué)的計算運算時的物體坐標變化如圖5所示。圖5中，t0為初始計算時刻，此時物體的方向坐標系(local coordinate system，LCS)為初始坐標系，方向坐標系原點的全局坐標(global coordinate system，GCS)為物體初始化時質(zhì)心位置。經(jīng)過一段時間的運動后，物體沿一定運動軌跡運動至新的位置，此時t1時刻的方向坐標系原點的全局坐標、方向角發(fā)生一定改變，此時用向量和方向角方式表示為

(1)

式(1)中：dt為計算時間步，dt=t1-t0；L為運動向量；Lt為下一計算循環(huán)的運動向量；V為運動速度向量；O_AP為方向坐標系原點的全局坐標；APt為下一計算循環(huán)方向坐標系原點的全局坐標。

通過每一計算循環(huán)對運動向量的修改，完成物體的移動計算，物體的轉(zhuǎn)動計算公式為

Ot+Δt=OtMt

(2)

Mt(Aaxis,θ)=

(3)

(4)

式中：Ot+Δt為下一計算循環(huán)步t+Δt時刻的物體方向坐標系的矩陣；Ot為當前計算循環(huán)時物體方向坐標系的矩陣描述；Mt為當前計算循環(huán)時的旋轉(zhuǎn)矩陣；旋轉(zhuǎn)矩陣可通過旋轉(zhuǎn)軸Aaxis和改變角θ表示；ωt為當前循環(huán)的角速度向量，可表示為

(5)

式(5)中：α為角加速度；T為離散元計算得到的滾刀轉(zhuǎn)矩矩陣；ωt+Δt為上一計算循環(huán)步t-Δt時刻的物體交速度向量；I為慣性矩陣。

最終的計算流程圖如圖6所示，其中RBD計算使用基于C++的EDEM二次開發(fā)實現(xiàn)，DEM則由EDEM本身進行。

3 數(shù)值模擬分析

3.1 數(shù)值模擬模型

圖8中為100、150、200、300 kPa圍壓下的實驗室與離散元三軸試驗對比，σ1、σ3為大主應(yīng)力與小主應(yīng)力。

在計算中使用17寸標準的滾刀作為滾刀計算對象，在具體計算中將滾刀詳細尺寸進行簡化，尺寸如圖9所示。

離散元模型如圖10所示，模型顆?？傆?56 131個，模型尺寸為1 m(長)×0.2 m(寬)×0.1 m(高)，計算時間步為1.274 24×10-5s。在計算中，模型除頂部無約束外，其余模型面均使用平面幾何體進行模型的約束固定。墻體、滾刀與顆粒的接觸采用Hertz-Mindlin理論，墻體彈性模擬取60 MPa、泊松比取0.3，滾刀彈性模量取206 GPa、泊松比為

圖6 RBD-DEM耦合計算流程Fig.6 Calculation process of RBD-DEM coupling numerical simulation

表1 離散元計算參數(shù)

圖7 離散元三軸試驗?zāi)MFig.7 Simulation of triaxial test with DEM

σ1為大主應(yīng)力；σ3為小主應(yīng)力圖8 離散元模擬三軸和室內(nèi)試驗數(shù)據(jù)對比Fig.8 Comparison of DEM and practical test

圖9 17寸滾刀簡化尺寸圖Fig.9 Simplified dimension drawing of 17 inch disc cutter

圖10 離散元計算模型Fig.10 computational model in DEM

0.25，墻體-顆粒的摩擦系數(shù)取0.5、滾動摩擦系數(shù)取0.5，滾刀-顆粒的摩擦系數(shù)取0.2、滾動摩擦系數(shù)取0.01。

滾刀從模型A側(cè)向B側(cè)線性移動，移動速度按工況設(shè)定，滾刀轉(zhuǎn)動由RBD耦合計算程序控制。

3.2 不同貫入度下數(shù)值模擬結(jié)果分析

研究表明，在華南地區(qū)的復(fù)合地層盾構(gòu)施工過程中，滾刀貫入度在1～16 mm，隨著硬巖比例的增加而降低貫入度[22]。為探究不同貫入度下滾刀的動態(tài)運動方式，分別設(shè)置2、4、7、8、9、10、12、14 mm貫入度工況進行計算，得到的滾刀轉(zhuǎn)速如圖11所示?？梢钥闯?，隨著貫入度的不斷提高，滾刀響應(yīng)速度逐漸邊塊，滾刀能更快地達到理論轉(zhuǎn)速附近，處于較好的動態(tài)運動狀態(tài)。但是從圖11中可以發(fā)現(xiàn)，在低貫入度時(2、4 mm)，滾刀一直處于較低的轉(zhuǎn)速，沒有進入理想的動態(tài)平衡，滾刀和巖土體間將具有較大的界面相對速度，符合圖3的猜想。但當滾刀的貫入度達到10 mm及以上時，能較好地達到動態(tài)平衡。計算最終各工況模擬得到的滾刀滑移率，如圖12所示，其中滑移率S計算方式為

(6)

式(6)中：S為滑移率；ω′為實際轉(zhuǎn)動角速度；ωo為理論轉(zhuǎn)動角速度。

通過滑移率S可以表征滾刀在運動中滾動運動的占比，即界面相對速度的占比。滑移率S越高，說明界面相對運動速度的占比越高，當滑移率S達到1時說明物體處于完全滑動的運動狀態(tài)。

圖11 滾刀不同貫入度下動態(tài)響應(yīng)的變化規(guī)律Fig.11 The changing law of disc cutter’s dynamic response with different penetration

圖12 各貫入度下最終得到的滾刀滑移率Fig.12 Slip ratio in different simulated conditions

從圖12中可以發(fā)現(xiàn)，隨著貫入度的提高滑移率逐漸下降，當貫入度達到10 mm及以上時滑移率下降趨勢不再明顯，而是趨于緩和，此時滾刀的運動狀態(tài)中滑動運動占比較少，可以認為滾刀處于動態(tài)平衡當中。但是如果繼續(xù)加大貫入度，滾刀的滑移率下降不明顯，說明貫入度對于滾刀進入合適動態(tài)狀態(tài)的貢獻降低。因此，可以認為貫入度對于滾刀進入合適動態(tài)平衡的貢獻只存在于一定區(qū)間內(nèi)，當小于這個區(qū)間時滾刀存在較大的界面滑移速度；當貫入度大于這個區(qū)間時并不能進一步提高動態(tài)響應(yīng)效果。

3.3 不同速度下數(shù)值模擬結(jié)果分析

目前實際工程中，不同盾構(gòu)機型號及地層差異導(dǎo)致盾構(gòu)掘進時刀盤轉(zhuǎn)速存在較大不同。研究表明，在復(fù)合地層中盾構(gòu)機的刀盤轉(zhuǎn)速區(qū)間為1.2～2.0 r/min[23-26]，因此使用1.5 r/min作為數(shù)值分析時的刀盤轉(zhuǎn)速?？紤]目前盾構(gòu)機制造技術(shù)發(fā)展的日新月異、大直徑盾構(gòu)的普遍應(yīng)用，以6 m作為模擬計算半徑的上限，即在1.5 r/min時線速度約為1 m/s。結(jié)合計算機性能的約束，最終設(shè)定0.3、0.5、0.8、1.0 m/s為模擬工況，并采用上述計算得到的10 mm貫入度作為貫入度設(shè)定以保證滾刀出于合適的動態(tài)平衡。

由于滾刀均能進入合適的動態(tài)平衡狀態(tài)，因此控制滾刀運動的主要因素為滾刀受到的扭矩，提取得到各工況中滾刀扭矩時程分布如圖13所示。

圖13 滾刀不同速度時的扭矩時程分布Fig.13 Torque time history distribution of disc cutter at different speeds

從圖13中可以看出，由于使用離散元作為巖土對象的模擬手段，因此扭矩分布存在一定離散性，但是仍具有一定規(guī)律。滾刀受到的扭矩先隨滾刀運動而升高，到達某一峰值點附近后開始下降，隨后在零值附近波動。這規(guī)律與圖11中進入動態(tài)平衡工況顯示的角速度曲線相吻合，即滾刀的扭矩先增大后減小，在滾動階段則在零值附近波動以保持轉(zhuǎn)動的穩(wěn)定性。圖3顯示，滾刀扭矩的提升速度和線性運動速度有關(guān)，隨著線性運動速度的提高，滾刀達到扭矩峰值點附近的時間越短；同樣的，線性運動速度越高滾刀的峰值扭矩越大。

圖13中，各線性速度下的扭矩時程分布曲線規(guī)律與Cho等[17]進行巖石線性切割實驗測得的滾動力分布類似，與其采用理想滾動的數(shù)值模擬結(jié)果有較大不同。因此認為在數(shù)值分析中，使用RBD耦合控制滾刀運動優(yōu)化了數(shù)值分析邊界條件，更為接近現(xiàn)實情況。對圖13的扭矩的時程分布進行擬合，對比擬合函數(shù)峰值點和峰值距離的差異性如圖14所示。

圖14 各工況峰值扭矩和半峰寬值差異性對比Fig.14 Comparison of peak torque and half peak width under different working conditions

圖14顯示峰值扭矩和運動速度間基本成線性規(guī)律，對于線性速度越低的滾刀，其峰值扭矩更低，運動速度為1 m/s的滾刀峰值扭矩是運動速度為0.3 m/s時的1.86倍。半峰寬度的差異性顯示，速度越快的滾刀達到峰值扭矩的時間更低，運動速度為0.3 m/s的滾刀達到峰值扭矩時間是運動速度為1 m/s時的1.72倍。

結(jié)合圖13和圖14的結(jié)果，可以認為滾刀線性運動速度對于滾刀啟動快慢有一定影響。滾刀啟動峰值扭矩和運動速度正相關(guān)，運動速度較低的滾刀峰值扭矩更低；滾刀達到峰值扭矩的時間和運動速度負相關(guān)，運動速度較低的滾刀達到峰值扭矩的時間更長。因此在復(fù)合地層中，近隧道軸線中心處的滾刀對安裝精度、渣土和易性更為敏感，由于理想狀態(tài)下啟動峰值扭矩更低、響應(yīng)時間更長，外界的少許擾動將使?jié)L刀不正常滾動，最終造成異常磨損。

4 結(jié)論

使用RBD-DEM耦合數(shù)值模擬技術(shù)，實現(xiàn)了滾刀在數(shù)值模擬中的動態(tài)轉(zhuǎn)動，進一步使針對滾刀的數(shù)值模擬符合工程實際。通過不同貫入度、不同運動速度，研究了滾刀進入動態(tài)轉(zhuǎn)動的影響原因及響應(yīng)快慢，得出如下結(jié)論。

(1)存在最優(yōu)貫入度區(qū)間，低于該區(qū)間時滾刀不能進入動態(tài)轉(zhuǎn)動狀態(tài)，大于該區(qū)間時改善效果不明顯且增大滾刀受力。

(2)滾刀啟動時受到的扭矩先隨滾刀運動而升高，到達某一峰值點附近后開始下降，隨后在零值附近波動。

(3)滾刀啟動峰值扭矩和運動速度正相關(guān)，運動速度較低的滾刀峰值扭矩更低。

(4)滾刀達到峰值扭矩的時間和運動速度負相關(guān)，運動速度較低的滾刀達到峰值扭矩的時間更長。

因此在復(fù)合地層中，滾刀存在一個最優(yōu)貫入度區(qū)間使得滾刀能夠出于優(yōu)秀的動態(tài)響應(yīng)狀態(tài)且磨耗較低。同時近隧道軸線中心處的滾刀由于運動速度低，因此對安裝精度、渣土和易性更為敏感，因為理想狀態(tài)下啟動峰值扭矩更低、響應(yīng)時間更長，外界的少許擾動將使?jié)L刀不正常滾動，最終造成異常磨損。通過數(shù)值模擬驗證了RBD-DEM的合理性，為后續(xù)相關(guān)滾刀動態(tài)研究提供了參考；初步研究成果為盾構(gòu)刀盤刀具設(shè)計及盾構(gòu)機的操作提供了參考。