離散元法用于絞刀齒切削珊瑚礁巖過程的三維模擬

林挺 ，馮沛洪 ，馬源 ，郭濤 *

（1．上海交通大學(xué)船舶海洋與建筑工程學(xué)院，上海 200030；2．中交疏浚技術(shù)裝備國家工程研究中心有限公司，上海 201208）

絞吸挖泥船是一種典型的疏浚工程船，作業(yè)時(shí)絞刀通過橋架下放到水底，絞刀旋轉(zhuǎn)并松散、破碎水底的密實(shí)砂土或巖體，再由管路、泥泵等輸送系統(tǒng)抽吸、排送砂土、碎巖與水的混合漿體，實(shí)現(xiàn)疏?；虼堤畹哪康摹＝陙?，國內(nèi)疏浚行業(yè)發(fā)展迅速且不斷承接國際的港口或海洋工程，沿海疏?；驆u礁吹填項(xiàng)目中不乏挖掘珊瑚礁石的工程。未來，水下挖巖的工程量可能進(jìn)一步增加，設(shè)計(jì)專用挖掘礁巖、風(fēng)化巖等巖土的絞刀、刀齒很有必要。

當(dāng)前，挖巖絞刀的功率預(yù)估、結(jié)構(gòu)設(shè)計(jì)、強(qiáng)度校核主要依靠同類施工經(jīng)驗(yàn)，這很可能造成過量設(shè)計(jì)與匹配浪費(fèi)，也對(duì)新型挖巖絞刀的優(yōu)化與設(shè)計(jì)造成了困擾。但是絞刀齒與巖石接觸并破碎巖石的具體受力難以通過絞刀軸上宏觀表現(xiàn)來考察，即使通過物理實(shí)驗(yàn)也難以直接測(cè)量，而且物理實(shí)驗(yàn)具有耗時(shí)長、費(fèi)用大等缺陷。

離散元方法（Discrete element method，DEM）源起上世紀(jì)70年代，可將物質(zhì)處理為由不連續(xù)的顆粒組成的集合，通過求解各顆粒的運(yùn)動(dòng)方程實(shí)現(xiàn)模擬物質(zhì)的運(yùn)動(dòng)與變形，是巖土力學(xué)領(lǐng)域中的典型數(shù)值研究方法[1]，已在海底巖體破碎等研究中得到良好應(yīng)用[2]。本文將采用DEM三維求解絞刀挖掘珊瑚礁巖時(shí)刀齒上的受力情況，為挖巖絞刀的設(shè)計(jì)提供參考。

1DEM物理模型

處理彈性顆粒球形面的變形時(shí)，DEM廣泛采用考慮切向庫倫摩擦力[3]以及回復(fù)阻尼系數(shù)[4]的Hertz-mindlin模型。該模型中，相互接觸的球面1、球面2因彈性變形而產(chǎn)生的相互作用力為：

式中：Fn、Ft分別為法向力和切向力；E*、G*分別為等效楊氏模量、等效剪切模量；R*為等效半徑；δn、δt分別為兩球面的法向重疊距離與切向重疊距離，即兩球面的球心距離相較于初接觸時(shí)刻（此時(shí)兩球面點(diǎn)接觸，球面無變形，亦無彈性接觸力）的徑向變化量以及兩球接觸面中點(diǎn)相對(duì)于初接觸時(shí)的切向變化量；μs為顆粒間的靜摩擦系數(shù)。并有：

式中：υ為泊松比；其他各字母含義同上，下標(biāo)1、2分別代表不同球面?？梢?，兩球面間剪切力的

式中：ρ為顆粒密度。實(shí)際上對(duì)于高配位數(shù)的準(zhǔn)靜態(tài)顆粒堆積（比如三維球形顆粒），通常時(shí)間步小于 0．2tRay是合適的。

Hertz-Mindlin模型沒考慮顆粒間的黏聚力，并不適于表征具有黏聚力的巖石類物質(zhì)。Potyondy與Cundall在此基礎(chǔ)上提出了一種用于巖石的顆粒粘合模型[6]，該模型類似于在相近顆粒之間填充了“固體粘合劑”，可以較大程度上模仿巖石類物質(zhì)。本文使用該模型來表征珊瑚礁巖中各顆粒間的相互作用。該模型假定在初始時(shí)刻，若相鄰離散顆粒的球面中心距離小于接觸半徑Rc，則球面之間將存在一種粘合作用（bonding），本文稱為鍵。鍵的兩端分別為兩球面中心，顆粒發(fā)生相對(duì)運(yùn)動(dòng)時(shí)，作用鍵使得球面間存在類似線性彈簧（包括力與轉(zhuǎn)矩）的相互影響，直到鍵被破壞。鍵對(duì)顆粒的作用力變化關(guān)系滿足式（4）：

式中：v為顆粒間相對(duì)速度；S為單位面積鍵剛度，其量綱為N·m-3；A為鍵的橫截面積；t為時(shí)間間隔；下標(biāo)n、t分別代表法向、切向分量。作用鍵在特定條件下會(huì)被破壞而導(dǎo)致“粘合效應(yīng)”失效，這種破壞不可逆。作用鍵的失效準(zhǔn)則包括：鍵上承受了極限法向應(yīng)力（CNS）、極限切應(yīng)力（CTS）或者鍵達(dá)到極限長度。當(dāng)存在作用鍵時(shí)，時(shí)間步長需要考慮鍵力的傳播。

2 珊瑚礁巖的力學(xué)性能及DEM模擬參數(shù)

珊瑚礁是由珊瑚蟲鈣質(zhì)遺骸及其他造礁生物的石灰質(zhì)長期累積而成的巖土體，主體為礁灰?guī)r，在中國南部沿海地區(qū)分布廣泛。由于珊瑚蟲種類繁多，珊瑚礁石形成過程中環(huán)境不斷變化，礁灰?guī)r的巖體結(jié)構(gòu)復(fù)雜，其力學(xué)性能也差異明顯。中科院王新志等人[7]曾對(duì)永暑礁及渚碧礁的礁灰?guī)r進(jìn)行了采樣及力學(xué)研究，得到礁灰?guī)r力學(xué)特性如表1所示[7]。研究中還發(fā)現(xiàn)礁灰石內(nèi)部孔隙率較高，極限不能超過最大靜摩擦力，否則球面間將發(fā)生相對(duì)滑動(dòng)或轉(zhuǎn)動(dòng)。

為保證計(jì)算的準(zhǔn)確性，需要確保時(shí)間步明顯小于瑞利時(shí)間步——剪切波穿過一個(gè)顆粒的時(shí)間。瑞利時(shí)間步 tRay由式（3）[5]確定：干燥礁巖的天然密度約為（1．25～1．35）× 103kg/m3，孔隙率約為45%～55%；礁石壓裂后會(huì)存在殘余強(qiáng)度，隨著壓裂破壞的增加強(qiáng)度進(jìn)一步降低，這與脆性灰?guī)r不同。

表1 采樣礁灰?guī)r力學(xué)特性Table 1 The mechanical property of the sampled coral reef rocks

為使DEM模型能夠準(zhǔn)確模擬水下飽和珊瑚礁灰?guī)r的力學(xué)性能，按DL/T 5368—2007《水電水利工程巖石實(shí)驗(yàn)規(guī)程》中的實(shí)驗(yàn)規(guī)程，對(duì)礁石的巴西劈裂法測(cè)抗拉強(qiáng)度實(shí)驗(yàn)以及單軸抗壓強(qiáng)度實(shí)驗(yàn)進(jìn)行相同的DEM數(shù)值實(shí)驗(yàn)，并根據(jù)實(shí)驗(yàn)結(jié)果調(diào)整計(jì)算參數(shù)。計(jì)算中，基本離散粒子選為球形；為保證顆粒直徑總是小于試件直徑的1/10，球形粒徑設(shè)為4 mm；接觸半徑為3 mm；為克服球形顆粒易滾動(dòng)的缺陷，限制顆粒的旋轉(zhuǎn)自由度；基于Mohr-Coulomb強(qiáng)度理論，假定破碎的主要原因是巖體抗剪切強(qiáng)度不足。

經(jīng)過多次“選擇參數(shù)-計(jì)算-修改參數(shù)”的循環(huán)后，選用表2所示的參數(shù)，得到與表1中實(shí)驗(yàn)結(jié)果均值相近的礁巖模型。

表2 模擬礁灰?guī)r顆粒參數(shù)Table 2 Particle parameters of the simulated coral reef rocks

巴西劈裂法測(cè)抗拉強(qiáng)度的實(shí)驗(yàn)中，圓柱形試件的周面上、下兩側(cè)各有寬約4 mm的墊條緩慢向內(nèi)擠壓，圖1給出了實(shí)驗(yàn)前后的鍵變化。圖2給出了實(shí)驗(yàn)中試件承受拉力隨墊條位移的變化過程。圖中可見，礁巖在不足1 mm的變形量下即發(fā)生突然的張性斷裂，斷裂面基本經(jīng)過上下的墊條，其抗拉強(qiáng)度約為1．2 MPa。

圖1 巴西劈裂實(shí)驗(yàn)前、后鍵狀態(tài)Fig.1 Statues of the bonds before and after Brazilian splitting test

圖2 巴西劈裂數(shù)值實(shí)驗(yàn)計(jì)算結(jié)果Fig.2 Numerical results of the Brazilian splitting test

圖3給出了礁石單軸抗壓強(qiáng)度實(shí)驗(yàn)前、后的鍵變化。圖4給出了單軸抗壓強(qiáng)度的數(shù)值計(jì)算結(jié)果。從圖中可見，數(shù)值模擬的礁巖圓柱在軸向形變約1 mm時(shí)達(dá)到最大單軸抗壓強(qiáng)度，約為6 MPa。此后壓力隨著軸向形變?cè)龃蠖杆贉p少，這點(diǎn)與實(shí)驗(yàn)中壓力階梯式減小的真實(shí)情況不符，這可能是數(shù)值模擬中忽略了真實(shí)礁巖材料內(nèi)不均勻分布的孔隙所造成的。盡管數(shù)值結(jié)果在應(yīng)力衰減上存在偏差，但在極限應(yīng)力、極限形變上與實(shí)驗(yàn)情況非常相近，仍認(rèn)為當(dāng)前計(jì)算參數(shù)能夠表征珊瑚礁巖的主要力學(xué)性能。

圖3 單軸壓裂實(shí)驗(yàn)前、后鍵狀態(tài)對(duì)比Fig.3 Statues of the bonds before and after uniaxial compression test

圖4 單軸壓裂數(shù)值實(shí)驗(yàn)計(jì)算結(jié)果Fig.4 Numerical results of the uniaxial compression test

3 絞刀齒切珊瑚礁巖的三維計(jì)算與分析

絞吸船工作時(shí)，橋架上的絞刀探入水下并繞軸旋轉(zhuǎn)，帶動(dòng)刀齒、刀臂切削巖土。船體一邊由橫移鋼纜牽動(dòng)繞工作鋼樁轉(zhuǎn)動(dòng)，一邊由臺(tái)車推動(dòng)改變船體與工作鋼樁的距離，同時(shí)橋架受巖土的反作用力變化，導(dǎo)致其下探角度略微波動(dòng)。因考察局部刀齒挖掘過程，絞刀的變半徑、變深度的弧形運(yùn)動(dòng)可以近似簡(jiǎn)化為繞軸旋轉(zhuǎn)以及水平橫移的疊加，如圖5所示。

圖5 簡(jiǎn)化后絞刀運(yùn)動(dòng)示意圖Fig.5 Schematic diagram for the simplified cutter movement

DEM的計(jì)算量隨粒子數(shù)增加而指數(shù)級(jí)上升。受限于計(jì)算資源，本文僅計(jì)算了某大型六臂絞刀的第4排相鄰刀齒以大挖巖橫移速度挖掘礁巖時(shí)的受力情況。該絞刀額定功率2 000 kW，最大旋轉(zhuǎn)半徑約1．74 m。由于安裝半徑、安裝角度不同，不同刀齒的挖掘面積也不同。以本絞刀第4排刀齒為例，其相鄰刀齒錯(cuò)開排列，齒尖半徑分別為1．232 m與1．309 m，尖端半徑較小的齒（后稱齒A）在軸向上較鄰齒（后稱齒B）前突約0．1 m，相鄰齒尖的周向夾角分別為67°及53°。絞刀旋轉(zhuǎn)速度22 r/min、橫移速度15 m/min，圖6給出向上挖掘礁巖時(shí)相鄰的齒A、齒B受力的計(jì)算模型。圖中模型的顆粒數(shù)約為85 000個(gè)，虛線為齒尖的運(yùn)動(dòng)軌跡。由圖可見，受刀齒錯(cuò)開排列的影響，齒A在橫移面上的挖掘面明顯小于齒B，但齒A在軸向上突前，軸向挖掘深度明顯大于齒B。

圖6 計(jì)算刀齒切削過程模型Fig.6 The model for simulating the cutting processes

計(jì)算時(shí)，先由齒A破環(huán)礁巖表面并留下挖掘痕跡，考察后續(xù)相鄰的齒B、齒A繼續(xù)切削礁巖的受力情況。珊瑚礁巖按圖6中的樣式構(gòu)建以減少離散粒子數(shù)量，礁巖軸向厚度0．28 m，絞刀軸位于礁巖上沿0．7 m。

為了提高計(jì)算效率，根據(jù)式（1）、（3）、（4）以及相似準(zhǔn)則，調(diào)整計(jì)算參數(shù)如表3所示。其他計(jì)算設(shè)置與前文一致，計(jì)算結(jié)果與選用表2參數(shù)的結(jié)果相似，相似比例為1∶1 000，計(jì)算耗時(shí)為后者的1/30。計(jì)算的時(shí)間步長為2×10-6s，總時(shí)長為1．2 s，使用Intel I7-4790處理器滿負(fù)荷計(jì)算共耗時(shí)約30 h。圖7給出了計(jì)算后礁巖的破壞情況。

表3 等效計(jì)算時(shí)選擇的參數(shù)設(shè)置Table 3 Parameters used in the equivalent calculations

圖7 刀齒切削后礁巖模型Fig.7 The reef model after cutting with cutter teeth

圖7中，礁巖上刀齒所過處礁巖破裂，部分碎裂成粒子，部分仍保持團(tuán)塊狀，刀齒未經(jīng)過處的內(nèi)部也存在著明顯裂痕。在頂部，齒A在軸向上基本挖穿了礁巖；在底部，齒B的挖掘深度約為齒A的一半。

設(shè)定齒B在挖掘起點(diǎn)處為絞刀的0°相位角，相位角將隨著絞刀旋轉(zhuǎn)而增加。圖8給出了計(jì)算得到的齒A與其相鄰的齒B在挖掘過程中所承受的礁巖作用力隨絞刀旋轉(zhuǎn)的變化。

從圖中可見，齒A與齒B上承受的峰值力分別約為5×104N與9×104N。齒B上承受了較高的作用力，這是由于齒B旋轉(zhuǎn)半徑較大，齒A切削范圍內(nèi)的部分礁巖顆粒已被齒B破壞了內(nèi)部鍵。

圖8 刀齒上切削合力隨相位角變化Fig.8 Cutting forces on the teeth vary with rotating anlges

假設(shè)絞刀第4排的3對(duì)A、B齒在挖掘過程的受力情況相同，根據(jù)刀齒受力矢量和相對(duì)位置可以折算出該排刀齒整體所受的挖掘力矩，該值以120°為周期往復(fù)性變化。圖9給出了2個(gè)周期內(nèi)的估算力矩及擬合曲線?？梢姡g刀第4排刀齒承受的力矩大致在（1～9）×104N·m的范圍內(nèi)，在齒B開始切削后又旋轉(zhuǎn)了約42°時(shí)達(dá)到最大值，在齒A完成切削時(shí)達(dá)到最小值。

圖9 第4排刀齒挖掘力矩隨相位角變化Fig.9 The cutting torque of the fourth row teeth vary with rotating anlges

圖8中，刀齒B在相位角約為40°～50°時(shí)受力最大，計(jì)算此時(shí)間段內(nèi)的刀齒表面承受壓力的時(shí)均值，采用插值映射的方法將該值作為刀齒載荷施加在有限元軟件上，考察作用力可能引起的刀齒應(yīng)力分布及變形。刀齒材料選擇為某低合金耐磨鋼，計(jì)算結(jié)果如圖10、圖11所示。圖中，礁巖對(duì)刀齒的壓力主要集中在齒尖上，在齒B承受最大受力時(shí)，齒尖壓力峰值高達(dá)1．14×108Pa。此時(shí)，齒尖處將發(fā)生最大約0．54 mm的變形，齒背上加強(qiáng)筋處的von-Mises等效應(yīng)力最大，約1．68×108Pa，同理此處的安全系數(shù)最小，約為7．1。

圖10 刀齒最大受力時(shí)的變形Fig.10 The deformation of the cutter teeth under maximum force

圖11 刀齒承受的最大等效應(yīng)力Fig.11 The maximum equivalent stress on the cutter teeth

4 總結(jié)與討論

使用DEM三維模擬了絞刀第4排刀齒挖掘珊瑚礁巖的過程，量化求解了挖巖機(jī)具的受力情況，并對(duì)計(jì)算結(jié)果進(jìn)行了有限元分析，為挖巖絞刀或類似機(jī)械的設(shè)計(jì)、校核與優(yōu)化提供了一種數(shù)值解決方案，得到以下結(jié)論：

1）Hertz-Mindlin接觸模型與顆粒粘合模型，可以更有效地模擬珊瑚礁巖等巖體的內(nèi)聚力效應(yīng)，模擬結(jié)果更接近真實(shí)情況。

2）本文確定的DEM參數(shù)可以較準(zhǔn)確地表征珊瑚礁巖的力學(xué)性能，并能較高效地處理絞刀挖掘礁巖的受力、破碎問題。

3）計(jì)算的刀齒挖掘礁巖時(shí)的最大受力約為9×104N，最大等效應(yīng)力約為 1．68×108Pa，佐證了該刀齒的結(jié)構(gòu)設(shè)計(jì)與安裝角度滿足挖掘礁巖的強(qiáng)度要求。

受限于資源耗費(fèi)與計(jì)算效率，本文在顆粒形狀、粒子尺度、粒子參數(shù)、模型范圍等方面尚有不足，例如未使用不同形狀、非均勻粒徑的粒子以及非均勻的鍵值等方式來模擬礁巖內(nèi)的裂隙現(xiàn)象，進(jìn)一步研究工作可在這方面做出改善，實(shí)現(xiàn)更高精度的計(jì)算，并推廣到整個(gè)絞刀挖巖過程的模擬。