基于離散元的貨艙內(nèi)鐵精礦運動數(shù)值研究

李筱航，張建偉，吳文鋒，盧金樹

(浙江海洋大學船舶與海運學院，浙江舟山 316022)

國際海事組織于2011 年頒布了IMSBC 規(guī)則，其將固體散裝貨物分為A、B、C 三組，其中A 組為易流態(tài)化貨物。流態(tài)化是指散貨船在航行過程中遭遇海浪產(chǎn)生擺動與振動，貨物內(nèi)部水分滲出，形成自由液面并發(fā)生運動，船舶復(fù)原力矩減小而穩(wěn)性降低，導(dǎo)致事故發(fā)生。針對該類問題，周健等[1]采用PFC3D 對散裝鐵精礦流態(tài)化的宏觀情況進行了研究。SPANDONIDIS，et al[2]基于分子動力學方法對貨艙中固體顆粒的靜止角、材料屬性與運動特性間的聯(lián)系進行了研究。宋喜慶等[3]結(jié)合離散元與有限元方法，對顆粒運動特性進行了研究。李文頡等[4]通過靜傾斜實驗得出含水率和庫倫摩擦間的關(guān)系，并基于離散元法通過改變含水量，提出了船舶穩(wěn)性數(shù)值評估方法。CHEN Wei，et al[5]建立了新的貨物臨界應(yīng)力條件模型，為評估貨物穩(wěn)定性提供了一種手段[5]。DE BLASIO[6]通過增大顆粒間作用力來模擬黏性力，在分子動力學基礎(chǔ)上建立了一套黏性物質(zhì)顆粒的流動方法等。

目前，易流態(tài)貨物在貨艙內(nèi)的運動規(guī)律研究日趨深入。為更清晰、直觀地展示艙內(nèi)顆粒貨物運動情況及主要影響因素，基于離散元法，設(shè)計了共11 種計算工況，研究顆粒間摩擦系數(shù)、粒徑級配、顆粒碰撞恢復(fù)系數(shù)、貨艙橫搖幅值及裝載深度等因素對艙內(nèi)鐵精礦顆粒運動的影響，以期為易流態(tài)貨物海上安全運輸提供參考。

1 數(shù)值模型建立

1.1 理論模型

根據(jù)牛頓第二定律，顆粒i 的運動方程如下：

式中：mi—顆粒i 的質(zhì)量；Ii—顆粒i 的轉(zhuǎn)動慣量；—顆粒i 的加速度；—顆粒的角加速度；F—顆粒在質(zhì)心處所受的合力；M—顆粒在質(zhì)心處所受的合外力矩[7]。

顆粒間接觸模型采用Hertz-Mindlin (no slip)接觸模型[8]。其中，顆粒間接觸法向力Fn可由下式計算：

式中E*為等效彈性模量，R*為當量半徑，α 為接觸半徑。

利用中心差分法對式(1)進行數(shù)值積分，可得到兩次迭代時間中間點表示的更新速度為：

式中：Δt—時間步長；N—對應(yīng)的時間t。

對式(3)進行積分，可得到如式(4)所示的位移等式：

將計算得到的位移帶入位移-力的關(guān)系中可以獲得新的作用，然后反復(fù)迭代計算，可獲取每一時刻顆粒受到的力和所處的位置。計算流程如圖1 所示，首先獲得初始接觸力Fi根據(jù)牛頓第二定律求出ui，然后根據(jù)位移和力的關(guān)系再得出Fi，繼而進行迭代計算。

圖1 離散元方法流程圖Fig.1 Flow chart of discrete element method

1.2 模型驗證

基于文獻[1]進行本文數(shù)值模型驗證，仿真參數(shù)見表1。三維方箱做x 方向的簡諧往復(fù)運動，頻率1 Hz，振幅60 mm。仿真對比結(jié)果如圖2 所示，從顆粒整體運動表面輪廓可見，本文數(shù)值計算結(jié)果與文獻結(jié)果基本一致?？赡艿恼`差在于：顆粒生成方式隨機，顆粒的初始狀態(tài)不能完全一致。

表1 模型驗證仿真參數(shù)設(shè)置Tab.1 Parameter settings of model validation

圖2 結(jié)果對比：(a)本文數(shù)值計算結(jié)果(b)文獻結(jié)果[1]Fig.2 Comparison of results：a：the present results b：results from reference[1]

2 研究方案

2.1 物理模型

以57 000 DWT 散貨船典型貨艙為原型[9]，考慮到數(shù)值模擬的局限性，考慮幾何縮尺比為1：660，實船貨艙及模型尺寸見表2 及圖3 所示。

表2 57 000 DWT 散貨船典型艙段及計算模型主尺度Tab.2 Typical sections of a 57 000 DWT bulk carrier and main dimensions of calculation model

貨艙模型如圖3 所示。

圖3 貨艙數(shù)值模型示意Fig.3 The cargo hold model

2.2 仿真工況

鐵精礦是我國主要的進口礦石貨種，也是易流態(tài)貨物海運事故高發(fā)貨種，與其他A 組貨物相比，其具有較高的比重，顆粒具有多層、微孔結(jié)構(gòu)等特點[10]。鑒于數(shù)值模型驗證的可獲得性，研究選取鐵精礦作為研究對象，其顆粒屬性參數(shù)見表1。考慮到易流態(tài)化貨物的材料屬性及船舶在海上運輸?shù)耐獠織l件等因素，選取橫搖幅值、貨物裝載情況、摩擦系數(shù)、粒徑級配及顆粒碰撞恢復(fù)系數(shù)進行研究。具體參數(shù)如表3 所示。

表3 計算條件設(shè)置Tab.3 Simulation conditions

利用User Defined 模塊配置粒徑級配，級配①為0.45 mm (10%)、0.6 mm (70%)、0.85 mm(20%)；級配②為0.45 mm (20%)、0.6 mm (60%)、0.85 mm (10%)；級配③為0.45 mm (10%)、0.6 mm (60%)、0.85 mm (30%)。運動形式為橫搖，由于實際運輸中顆粒運動比較緩慢，為加速顆粒運動，設(shè)定橫搖周期為1 s，搖擺中心為形心。時間步長20%，網(wǎng)格大小為顆粒的3 倍最小半徑，網(wǎng)格數(shù)為101 250 個。

3 結(jié)果分析

3.1 摩擦系數(shù)對鐵精礦顆粒運動的影響

不同摩擦系數(shù)下顆粒運動情況如圖4，可知工況1 顆粒傾斜程度最大，工況2 次之，工況3 顆粒傾斜程度最小，顆粒幾乎沒有流動?？梢?，顆粒庫倫摩擦系數(shù)越大，顆粒越不易發(fā)生流動。

圖4 9.75 s 時工況1(左)、2(中)、3(右)的顆粒分布情況Fig.4 Particle distribution of condition 1 (left),2 (middle) and 3 (right) at 9.75 s

顆粒法向接觸力為顆粒的法向沖擊力，其越大代表顆粒的碰撞次數(shù)越多、頻率越高、顆粒運動更劇烈，顆粒的動能是衡量顆粒運動速度的重要標志。特殊強調(diào)的一點是，當顆粒靜止不動時，也就是當時間為0時，在重力作用下，顆粒之間的法向接觸力為0.000 88 N。

工況1、2、3 顆粒法向接觸力、動能和重心偏移量對比如圖5 所示?？芍Σ料禂?shù)越小，顆粒間的法向接觸力越大，顆粒的碰撞更劇烈，顆粒的法向接觸力隨著摩擦系數(shù)的增大而減少。摩擦系數(shù)改變后，顆粒的動能隨著摩擦系數(shù)的增大而減小。隨著摩擦力的增大，顆粒的重心偏移量變化越小，且在單個工況下，隨著時間的增加，在慣性的作用下顆粒難以返回初始位置，顆粒的行程會越來越大，導(dǎo)致貨物整體重心偏移量變大。

圖5 工況1、2、3 時法向接觸力(左)、顆粒動能(中)和重心偏移量(右)變化曲線Fig.5 Variation of normal contact force (left),particle kinetic energy (middle) and center of gravity deviation (right) under conditions 1,2 and 3

3.2 橫搖幅值對鐵精礦顆粒運動的影響

工況1、4、5 的仿真結(jié)果對比如圖6 所示，橫搖幅度越大，法向接觸力變化越明顯，峰值越大。三種工況隨著橫搖幅值的減小，顆粒動能遞減，貨物整體重心偏移量明顯減小。

圖6 工況1、4、5 時法向接觸力(左)、顆粒動能(中)和重心偏移量(右)的變化曲線Fig.6 Variation of normal contact force (left),particle kinetic energy (middle) and center of gravity deviation (right) under conditions 1,4 and 5

3.3 粒徑級配對鐵精礦顆粒運動的影響

工況1、6、7 的數(shù)據(jù)對比如圖7。半載情況下，工況1、6 和7 分別裝載顆粒4 450 個、7 810 個和6 450個。從圖7 可以看出顆粒數(shù)量與顆粒碰撞程度成反比。工況6 與工況7 的動能峰值相差不大，工況6 的顆粒動能峰值略大，工況1 的顆粒動能峰值遠遠高于其他兩個工況。裝載深度相同的情況下，顆粒越多重心越高，但重心偏移量相差不明顯。

圖7 工況1、6、7 下法向接觸力(左)、顆粒動能(中)和重心偏移量(右)的變化曲線Fig.7 Variation of normal contact force (left),particle kinetic energy (middle) and center of gravity deviation (right) under conditions 1,6 and 7

3.4 碰撞恢復(fù)系數(shù)對鐵精礦顆粒運動的影響

工況1、8、9 的數(shù)據(jù)對比如圖8，可知顆粒法向接觸力幾乎沒有差別，并且3 種工況下的顆粒動能變化曲線完全一樣。由此可見，不同的碰撞恢復(fù)系數(shù)對顆粒的法向接觸力和動能并沒有較大影響。3 種工況的重心偏移量也基本相同。

圖8 工況1、8、9 時法向接觸力(左)、顆粒動能(中)和重心偏移量(右)的變化曲線Fig.8 Variation of normal contact force (left),particle kinetic energy (middle) and center of gravity deviation (right) under conditions 1,8 and 9

3.5 貨物裝載情況對鐵精礦顆粒運動的影響

工況1、10、11 的數(shù)據(jù)對比如圖9，可知貨物裝載量與顆粒的碰撞程度成正比，工況1 和工況10 的法向接觸力峰值基本相同，但工況1 的法向接觸力變化幅度大，說明工況10 在高裝載率下顆粒碰撞較為劇烈。工況1 和工況10 的顆粒動能幾乎一樣，工況1 的顆粒動能維持在一個較低的水平。工況11 顆粒的重心偏移量變化幅度更大，工況10 的重心偏移量變化幅度最小。

圖9 工況1、10、11 法向接觸力(左)、顆粒動能(中)和重心偏移量(右)的變化曲線Fig.9 Variation of normal contact force (left),particle kinetic energy (middle) and center of gravity deviation (right) under working conditions 1,10 and 11

4 結(jié)論

本文通過與其他文獻結(jié)果進行對比驗證了模型的正確性，確定了11 種工況，對散貨船貨艙內(nèi)的鐵礦石顆粒在不同工況下的運動情況進行了數(shù)值模擬，通過改變摩擦系數(shù)來代替不同含水量的影響，得出以下結(jié)論：

(1)本研究通過改變顆粒間摩擦系數(shù)來替代含水量的影響具備正確性和可行性，本文構(gòu)建的鐵精礦顆粒艙內(nèi)運動數(shù)值仿真模型能夠正確處理顆粒材料在艙內(nèi)的運動仿真。

(2)顆粒摩擦系數(shù)增加即貨物含水量增加，則顆粒運動幅度明顯降低；單個工況下，由于慣性的作用，顆粒的運動幅度隨時間的延長不斷增加。橫搖幅值的增加會導(dǎo)致顆粒碰撞更加劇烈，顆粒的重心偏移量變化幅度會更大；碰撞恢復(fù)系數(shù)的改變對顆粒的運動情況幾乎沒有影響；在裝載體積一定的情況下，顆粒越少的粒徑級配，顆粒的運動幅度越大；隨著顆粒裝載的增多，顆粒的法向接觸力和動能會維持在一個更高的水平。