徐汪豪， 姚志剛， 方 勇， 王飛陽(yáng)， 周 健

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

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

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

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

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

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

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

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

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

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

∝表示角速度與線(xiàn)速度的匹配圖2 滾刀的運(yùn)動(dòng)狀態(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ù)值模擬方法

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

圖5 剛體動(dòng)力學(xué)計(jì)算時(shí)的物體坐標(biāo)變化Fig.5 Object coordinate change in RBD numerical simulation

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

(1)

式(1)中：dt為計(jì)算時(shí)間步，dt=t1-t0；L為運(yùn)動(dòng)向量；Lt為下一計(jì)算循環(huán)的運(yùn)動(dòng)向量；V為運(yùn)動(dòng)速度向量；O_AP為方向坐標(biāo)系原點(diǎn)的全局坐標(biāo)；APt為下一計(jì)算循環(huán)方向坐標(biāo)系原點(diǎn)的全局坐標(biāo)。

通過(guò)每一計(jì)算循環(huán)對(duì)運(yùn)動(dòng)向量的修改，完成物體的移動(dòng)計(jì)算，物體的轉(zhuǎn)動(dòng)計(jì)算公式為

Ot+Δt=OtMt

(2)

Mt(Aaxis,θ)=

(3)

(4)

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

(5)

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

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

3 數(shù)值模擬分析

3.1 數(shù)值模擬模型

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

在計(jì)算中使用17寸標(biāo)準(zhǔn)的滾刀作為滾刀計(jì)算對(duì)象，在具體計(jì)算中將滾刀詳細(xì)尺寸進(jìn)行簡(jiǎn)化，尺寸如圖9所示。

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

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

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

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

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

圖9 17寸滾刀簡(jiǎn)化尺寸圖Fig.9 Simplified dimension drawing of 17 inch disc cutter

圖10 離散元計(jì)算模型Fig.10 computational model in DEM

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

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

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

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

(6)

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

通過(guò)滑移率S可以表征滾刀在運(yùn)動(dòng)中滾動(dòng)運(yùn)動(dòng)的占比，即界面相對(duì)速度的占比。滑移率S越高，說(shuō)明界面相對(duì)運(yùn)動(dòng)速度的占比越高，當(dāng)滑移率S達(dá)到1時(shí)說(shuō)明物體處于完全滑動(dòng)的運(yùn)動(dòng)狀態(tài)。

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

(2)滾刀啟動(dòng)時(shí)受到的扭矩先隨滾刀運(yùn)動(dòng)而升高，到達(dá)某一峰值點(diǎn)附近后開(kāi)始下降，隨后在零值附近波動(dòng)。

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

(4)滾刀達(dá)到峰值扭矩的時(shí)間和運(yùn)動(dòng)速度負(fù)相關(guān)，運(yùn)動(dòng)速度較低的滾刀達(dá)到峰值扭矩的時(shí)間更長(zhǎng)。

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