蓮藕主藕體彎曲破壞離散元仿真分析

焦俊，張國忠,2，杜俊,2，劉浩蓬,2，查顯濤,2，邢赫,2

1.華中農業大學工學院，武漢 430070； 2.農業農村部長江中下游農業裝備重點實驗室，武漢 430070

蓮藕(NelumbonuciferaGaertn)，簡稱蓮，別名蓮菜、荷藕等，是一種用途十分廣泛的水生經濟作物，也是我國種植面積最大的水生蔬菜[1]。我國蓮藕種植主要分布于湖北、江蘇、山東、安徽、福建、湖南、浙江等省。近年來，我國蓮藕種植區域逐步擴大，種植面積亦呈穩定增長態勢。蓮藕具有很好的食用和藥用價值，蓮藕和蓮子可供食用，花粉、荷葉、蓮芯等也都可以作為菜肴、飲料及保健食品。蓮藕地下莖為主要食用部位，其由主藕、支藕組成，無序生長于水下淤泥中，斜植交叉，無規律隨機分布。其中，主藕由3～7節膨大節間組成，總長度可達1.2～1.4 m，質量達3.0～4.0 kg[2]。為保證成品合格率，現階段蓮藕采收多用人工挖掘，先由人工腳踏摸索出藕枝大致位置及生長方向，再用手持高壓水槍沿蓮藕外側環繞切割、粉碎淤泥，隨后人工輔助拔出主體蓮藕，使其與淤泥分離[3]。受限于蓮藕生長無規律及子藕旁支錯綜復雜等問題，拔出過程中若部分蓮藕枝節仍埋于淤泥中，易發生彎曲、斷裂，泥水會經斷裂處進入內孔，致使蓮藕品質受損。提高效率并減少收獲過程對蓮藕的損傷與破壞是當前蓮藕機械化采收領域亟待解決的技術難題。

采用三維建模以及數值計算仿真軟件對采挖作業進行仿真分析，是優化改進機械化收獲技術與設備的有效途徑[4-6]。蓮藕作為一種狹長枝節狀的多孔質生物物料，其力學仿真模型缺乏研究。近年來，有關農業物料力學特性的研究如肥料、馬鈴薯[7]、三七種子等散粒以及莖稈類農業物料的離散元建模和仿真參數標定的報道逐漸增多。于慶旭等[8]采用逆向工程技術對三七種子進行了三維建模，通過碰撞試驗、斜面試驗、圓筒提升堆積試驗標定其接觸參數，并開展非規則散粒物料的仿真模擬。馬文鵬等[9]以苜蓿種子休止角與堆積角開展雙指標多目標尋優計算得到其離散元接觸參數，并通過槽輪式排種器進行試驗驗證。張濤等[10]采用抽板試驗測定了玉米秸稈徑向堆積角，以其為目標采用正交試驗標定了接觸參數。廖宜濤等[11]以不同直徑的飼料油菜薹期收獲莖稈為對象，結合圓筒提升堆積試驗以及彎曲破壞仿真試驗分別標定其接觸參數以及粘結參數，建立莖稈的離散元力學模型。

本研究以鄂蓮5號蓮藕為對象，采用物理試驗測定其主藕體本征參數、接觸參數以及內孔尺寸，利用三維反求技術獲取主藕體輪廓以及放樣切割內孔方法完成主藕體三維建模，采用EDEM內置的Hertz-Mindlin with bonding粘結模型建立主藕體離散元模型，進而開展彎曲破壞規律仿真分析，以期為蓮藕機械化收獲過程仿真研究以及蓮藕損傷、破壞途徑與因素的分析提供參考。

1 材料與方法

1.1 試驗材料

試驗材料為湖北省廣泛種植的蓮藕(NelumbonuciferaGaertn)鄂蓮5號。選取第二節主藕體為對象進行仿真分析，取樣區域與部位如圖1所示。

圖1 蓮藕取樣區域Fig.1 Sampling diagram of lotus root

1.2 本征參數測量

1)軸向壓縮試驗。為測定主藕體彈性模量、剪切模量、泊松比等本征參數，選取蓮藕表層和心部等部位制作5 mm×5 mm×5 mm立方形標準試樣，在TMS-Pro質構儀進行平板軸壓縮試驗(圖2)。設置壓縮速度為10 mm/min、加載位移為4 mm，重復10次，通過測量試樣單軸壓縮前后高度與直徑變化量并參考文獻[12]方法計算獲得主藕體彈性模量、剪切模量、泊松比等本征參數。

圖2 單軸壓縮試驗Fig.2 Uniaxial compression test

2)主藕體-鋼摩擦系數測定[10]。采用斜面法獲取主藕體-鋼摩擦系數。試驗前根據前期蓮藕外形尺寸抽樣調查結果，選取橫截面長軸為(70±5) mm、短軸為(60±5) mm的第二節蓮藕，去除兩端，選取中間部位制備成長100 mm筒狀主藕體試樣。靜摩擦系數測量時，為保證試樣在斜面上僅能產生滑動，將2節主藕體沿垂直縱軸方向插入直徑5 mm、長20 mm鋼條，形成并聯結構擺放(圖3A)，緩慢推動斜面至其發生滑動，記錄該時刻斜面與水平面夾角，重復20次取平均值，計算得到主藕體-鋼靜摩擦系數。測量主藕體-鋼滾動摩擦系數時(圖3B)，取1節試樣沿軸向方向放置于斜面上，按照靜摩擦試驗方法試驗，試驗同樣重復20次，最后取平均值計算滾動摩擦系數。

A:滑動試驗 Coefficient of static friction test; B：滾動試驗 Coefficient of rolling friction test.

3)主藕體-主藕體摩擦系數測定。目前尚無對主藕體間摩擦系數測量的文獻報道。本研究參照GB/T 10006―2021《塑料 薄膜和薄片 摩擦系數的測定》標準，采用摩擦系數測試儀進行測量。測量時，沿主藕體周向均勻切取若干個規格為8 cm×20 cm×1 cm切片，并采用粘結劑將切片內表面均勻依次密布粘在滑塊以及試驗臺上，使其外表面相對所采用粘結劑對主藕體外表面無滲透硬化作用(圖4)。粘貼完畢后，驅動滑塊以100 mm/min的速度相對試驗臺直線運動，牽引方向與摩擦方向平行，以模擬主藕體間滑動、滾動等相對運動情形，測量計算主藕體間滾動摩擦、靜摩擦系數。試驗時控制試驗時間，確保試驗環境條件恒溫恒濕，保持主藕體以及切片水分恒定，試驗前后采集同一部位小塊主藕體試樣，采用水分測定儀測定其含水率，均穩定保持在75%±5%。

圖4 主藕體-主藕體摩擦系數測定Fig.4 Measurement of friction coefficient between main lotus root and main lotus root

4)主藕體碰撞系數測定。為測定主藕體-鋼碰撞恢復系數，制備與前述摩擦試驗規格一致主藕體試樣，在其表面涂上粉末以用于標記碰撞點。基于壓縮試驗結果以及防止試驗時損傷主藕體，調整臺架高度H為 20 cm(圖5)。試驗時主藕體蓮藕在高度H處經限位環投放，自由落體至與碰撞板發生碰撞后做拋物線運動落入沙箱。同一試樣多次重復試驗后取其集中落點并采集落點距碰撞點水平距離S1，以及無底座時主藕體落入沙盤的水平距離S2，按照式(1)計算碰撞后水平、垂直方向分解速度vx、vy，按照式(2)計算主藕體-鋼碰撞恢復系數[13]。

(1)

(2)

圖5 碰撞模型Fig.5 Model of collision

為測定主藕體間碰撞恢復系數，開展主藕體間對心碰撞試驗。試驗時，在縱軸高度20 cm處釋放主藕體，其與下方固定狀態主藕體試樣碰撞，以下方主藕體試樣上表面軸向切線方向為橫軸，垂直方向為縱軸，建立如圖6所示坐標系，采用高速攝影儀記錄上方主藕體碰撞后沿縱向彈跳高度，計算二者比值得到主藕體間碰撞恢復系數[14]。

圖6 主藕體對心碰撞試驗Fig.6 Central collision test between mainlotus root and other main lotus root

5)主藕體本征與接觸參數。經以上試驗得到主藕體本征與接觸參數，結果如表1所示。

表1 主藕體本征與接觸參數 Table 1 Intrinsicand contact parameters of main lotus root

1.3 主藕體仿真模型建立

1)三維模型獲取。主藕體結構不規則，其形狀參數直接影響粘結模型中粘結鍵數量以及粘結結構。為此，采用逆向工程技術，通過EinScan-pro手持式3D掃描儀對其掃描，封裝點云數據，從而得到精確主藕體輪廓模型[7]。采用Geomagicstudio 3D軟件，將點云數據轉換為多邊形，刪除尖狀物、填充孔，去除多余特征，針對反映主藕形狀特征的輪廓銳化處理，生成柵格，擬合曲面得到主藕NURBS輪廓曲面模型。由于主藕體內部存在貫通內孔(圖7)，其多孔質結構直接影響物料機械強度，故將主藕體內孔依次標號后統計測量分析其形狀參數(表2)，導入Solidworks三維軟件，放樣切割得到主藕體三維實體模型。

圖7 主藕體截面Fig.7 Cross section of main lotus root

表2 主藕體內孔尺寸 Table 2 Size of inner hole of main lotus root mm

取顆粒半徑3 mm，以外表面輪廓向內縮放2倍顆粒半徑距離建立壁面，將主藕體模型分為內外兩層，用于后續建立外表皮以及心部的顆粒材料(圖8)。

A:點云數據 Point cloud data; B:輪廓曲面模型 Profile surface; C:實體模型 Solid model.

2)顆粒填充與粘結。將主藕體顆粒模型分為心部與外表皮顆粒模型。如圖9A所示，外表皮模型利用Workbench構建單層Prism網格模型后導入Fluent中，利用udf計算得到空心六邊形密排面(圖9B)，構建雙層prism網格模型后分別導出內外表面以及整體節點坐標信息，篩選去除重復坐標信息后得到中層節點信息，將其作為正六邊形中心，最終得到完整密布形外表皮顆粒模型，后處理輸出外表皮XML文件[15]。

設置壁面模型為Hertz-Mindlin無摩擦(no slip)模型，采用落雨法生成心部顆粒，設置box為physical，作為顆粒工廠生成心部顆粒后，利用平板向下壓縮，保證顆粒密實堆積，而后主藕體內層幾何模型由初始Virtual狀態變為physical狀態，沿壁面外形束縛顆粒。壓縮過程中存在部分顆粒進入內孔，設置內孔模型以100 mm/s速率沿軸向直線運動，同時設置接觸幾何與多余顆粒接觸，采用Remove particle API，清理內孔中多余顆粒，后處理輸出心部XML文件。提取2份XML文件中心部與外表皮顆粒坐標信息，輸入XML文件中得到完整填充文件，將完整填充文件輸入EDEM。

若壁面彈性模量過大，會造成顆粒間短時間難以達到平衡狀態，而過小也會造成顆粒在加載壓實過程中穿透壁面使得顆粒無法貼合非球形壁面。為獲取反映蓮藕形狀特征的精準粘結模型，取壁面彈性模量為顆粒彈性模量的1/10，在顆粒間無明顯重疊下重新排列填充直至達到平衡狀態，最終得到主藕體顆粒模型(圖9C)。

A:網格模型 Mesh model; B:六邊形基礎模型 Hexagon base mode; C:完整模型 Complete model.

1.4 主藕體彎曲破壞離散元仿真

1)主藕體彎曲破壞試驗。如圖10A所示，采用TMS-Pro質構儀(FTC公司)進行主藕體彎曲破壞試驗，選用質構儀自帶長(A)×寬(B)×高(H)×厚(R,刃口)為50 mm×10 mm×60 mm×3 mm刀頭(圖10B)，以10 mm/min速率對主藕體施加載荷，記錄加載過程中隨加載位移變化的彎曲破壞力曲線，試驗規格與前述摩擦試驗一致，固定支撐兩點間距為50 mm，參照文獻[11]重復試驗10次，統計顯示最大破壞力平均值為315 N。由圖10C可知：(1)起始時，主藕體在加載過程中經過原點到A間為塑性形變階段，此時局部出現由斜裂紋擴展導致的剪切破壞以及軸向拉裂紋擴展而出現的劈裂破壞[15]；(2)隨后，在A-B階段出現應力集中現象，裂紋擴展迅速，對外表現為局部脆性斷裂；(3)而后，在B-C階段主藕體開始斷裂，裂紋擴張至藕心部，應力增加；(4)最后，壓頭繼續加載，藕體迅速斷裂，表現為典型楔入斷裂。

A:彎曲破壞試驗 Bending failure test; B:彎曲刀頭示意圖 Curved tool head; C:典型主藕體彎曲破壞應力-應變曲線 Typical bending failure curve of main lotus root.

2)粘結模型及參數計算。彎曲破壞試驗顯示主藕體產生劈裂現象，故視其為脆性材料。由于Hertz-Mindlin with bonding模型是1個粘結接觸模型，常用于模擬斷裂問題，故用其來模擬主藕體彎曲破壞過程[17]。Hertz-Mindlin with Bonding模型中顆粒間力學特性滿足梁彈性理論[18]，如式(3)所示：

(3)

其中，Kn、Ks分別為法向接觸剛度、切向接觸剛度，N/m；Fn為法向應力，成對顆粒接觸后產生切向應力增量ΔFs疊加在切向應力Fs上[19]。顆粒間受力以式(4)隨時間步長迭代計算得到[20]：

(4)

離散元模型需要標定的參數：法向粘結剛度kn，N/m3；切向粘結剛度ks，N/m3；切向臨界應力τmax，MPa；法向臨界應力σmax，MPa；粘結半徑R，m。為節約仿真計算時間，球形顆粒物理半徑r為3 mm，接觸半徑取4 mm，粘結半徑一般為接觸半徑的1.2～2.0倍，粘結半徑R取5 mm[21]。根據梁理論，計算得到粘結鍵之間的法向臨界應力和切向臨界應力:

(5)

式(5)中，當拉伸應力、剪切應力超過上述最大值后，粘結鍵斷裂。

3)主藕體彎曲破壞離散元仿真單因素試驗。按照前述測定的結果，開展主藕體彎曲破壞仿真試驗(圖11)。由于Hertz-Mindlin with bonding模型中法向粘結剛度、切向粘結剛度、法向臨界應力、切向臨界應力的粘結參數均對仿真結果存在影響，為探索其影響規律，設置外表皮、心部、外表皮粘結參數一致，采集初始段曲線，對法向粘結剛度x1、切向粘結剛度x2、法向臨界應力x3、切向臨界應力x4，分別設置106、107、1083個水平(表3)，重復10次，開展主藕體彎曲破壞仿真單因素試驗[22]。

圖11 主藕體彎曲破壞仿真試驗Fig.11 Bending failure simulation of main lotus root

表3 單因素試驗 Table 3 Single factor test

試驗時，設置顆粒接觸模型為Hertz-Mindlin with bonding模型，在外表皮顆粒、心部顆粒以及外表皮-心部顆粒間分別添加粘結鍵，輸入粘結參數，設置幾何壁面彈性模量為7.94×1010MPa，基于該模型對時間步長敏感，取低于10%時間步長，取整后為3×10-5s進行運算。

由圖12可知，法向、切向粘結剛度對曲線位移值以及第一峰值影響顯著，法向臨界應力、切向臨界應力對曲線峰值以及位移值影響不顯著。根據單因素試驗獲得法向粘結剛度、切向粘結剛度、法向臨界應力、切向臨界應力參數上下限區間如表4所示。

表4 粘結參數編碼 Table 4 The code of bonding parameter

圖12 單因素試驗Fig.12 Single factor test

4)主藕體彎曲破壞離散元仿真二因子試驗。根據單因素試驗結果，為進一步獲取上述因素對彎曲破壞曲線的影響規律以及二次效應(彎曲性)是否顯著[23]，開展了5個中心點的二水平因子試驗。

5)最速下降法。為快速趨近實際彎曲破壞曲線及關鍵破壞點，從二因子試驗獲取的一階響應曲面方程中心點出發，取沿垂直于擬合曲面等高線的直線為最速下降路徑，取步長與回歸系數成正比，按照最速下降法尋找法向粘結剛度、切向粘結剛度、法向臨界應力、切向臨界應力的最優解。

2 結果與分析

2.1 一階響應曲面分析

由單因素試驗可知，剛度參數同時對曲線第一峰值以及位移值影響顯著，為精確衡量仿真曲線與實測曲線相似性，采用式(6)以仿真以及實測曲線中第一峰值與位移點組成的坐標點，計算坐標間歐拉距離(ED)作為衡量仿真與實測試驗相對誤差的指標。

(6)

式(6)中，Sa、Sb分別為仿真、實測曲線坐標點，Fa、Fb分別為仿真、實測曲線第一峰值，Da、Db分別為仿真、實測曲線位移值。

試驗指標包含峰值以及位移點，兩者度量單位不一，需采用式(7)所示的闕值比較法進行無量綱化處理。

(7)

由式(6)和式(7)可得無量綱化的歐拉距離ED0。

(8)

以表4所示參數上下限數值設計二水平因子試驗(表5)。由表6可知：主效應x1、x2對ED0參數影響顯著(P<0.01)，主效應x3、x4、彎曲性(純二次效應)以及其余交互作用均對ED0參數影響不顯著(P>0.05)，彎曲性不顯著證明一階響應曲面方程具有適合性。使用最小二乘法，利用規范變量以一階模型擬合得到關于ED0參數一階響應曲面方程：

ED0=2.24-7.20×10-9x1-1.06×10-8x2+ 3.75×10-17x1x2

(9)

表5 二水平因子試驗 Table 5 Regular two-level factorial design

表6 一階曲面響應模型方差分析 Table 6 Analysis of variance of first order surface response equation

2.2 最速下降法分析

由一階響應曲面方程可知，x1、x2對ED0參數為負效應。非顯著項x3、x4取中間水平。以一階響應曲面方程中系數最大的x2作為基準值，取0.70 MPa為基本步長，由式(10)所示規范向量x1、x2與自然向量z1、z2關系，設計最速下降試驗(表7)。

(10)

表7 最速下降試驗 Table 7 Steepest ascent

由表7可知，在步長為原點+3Δ處ED0參數最低，該粘結參數下第一彎曲峰值為269.72 N、位移值為7.14 mm，與實測試驗數據相對誤差分別為2.56%、2.00%，總體誤差為2.28%，可以認為此時法向粘結剛度、切向粘結剛度、法向臨界應力、切向臨界應力等粘結參數為最優解，仿真與實測彎曲曲線對比如圖13所示。

圖13 仿真-實測曲線Fig.13 Curve of simulation-test

3 討 論

本研究建立了鄂蓮5號主藕體離散元模型，使用Hertz-Mindlin with bonding模型對其彎曲破壞模型進行了仿真分析，采用單因素試驗、二因子試驗、最速下降法得到法向粘結剛度為5.814×108N/m3、切向粘結剛度為3.450×108N/m3、法向臨界應力為3.80 MPa、切向臨界應力為3.12 MPa。第一彎曲峰值269.72 N、位移7.14 mm，相較實測曲線相差2.56%、2.00%。 對仿真模型各關鍵參數分析顯示，法向剛度粘結系數與切向粘結剛度對于第一彎曲力峰值以及峰值位移點影響顯著，而法向臨界應力與切向臨界應力對第一峰值以及位移點影響不顯著，其中法向剛度大于切向剛度，第一彎曲力峰值大于第二彎曲力峰值。

本研究僅對主藕體彎曲破壞塑形變形階段進行了仿真模擬分析，后期需要進一步對藕節以及主藕、藕節聯合體的受力破壞進行研究。此外，本研究未分析截面藕孔分布以及對藕體破壞性能的影響，后續也需進一步開展研究。