基于離散元的礦砂船貨物載荷計算方法

宋喜慶 張少雄 胡豐梁

(武漢理工大學交通學院1) 武漢 430063) (中國船級社2) 北京 100007)

0 引 言

超大型礦砂船(VLOC)與普通散貨船相比在船體結構型式上存在較大差異，近年來礦砂船諸多的傾覆事故及其不斷大型化的發展勢頭使其結構安全性備受關注，一般需采用船舶直接計算方法進行強度校核.載荷的計算是決定直接計算準確與否的關鍵，礦砂貨物載荷是其重要組成.

干散貨對于壁面的作用的研究屬于土力學范疇，擋土墻壓力理論其重要參考.在船舶工程領域， Koichi等[1-2]多次開展實船航行過程中貨物壓力進行監測，然而數據顯示礦砂船航行過程中受船體的振動、晃蕩等因素的影響，貨物壓力的實測結果與按庫倫土理論整理的貨物壓力存在較大差異，且無明顯規律.

目前船舶工程中所應用的載荷計算公式通常以規范公式為指導，如文獻[3]中對于干散貨貨物載荷的計算要求，通過側壓力系數計算非水平板的側向壓力.這種用于船舶強度直接計算的內部載荷計算方法不能考慮干散貨顆粒的形狀參數、艙壁的傾斜程度、滑裂面假設[4]等，也忽略了干散貨顆粒裝載過程中的隨機分布、顆粒間的摩擦系數、顆粒與艙壁的摩擦系數及顆粒的剪切模量，是一種近似的干散貨內部載荷的描述方法.因此有必要對船舶直接計算中傳統計算方法的可靠性進行重新考慮.

1 考慮方法

干散貨尤其是礦砂類貨物的的力學特征體現為“散”“動”，使用連續體力學分析干散貨壓力問題依賴了太多假設.通過離散元法(DEM)描述干散貨堆積狀態可以使貨物載荷與實際更加貼合.隨著計算機運算性能的快速提升，DEM計算方法的可行性將逐漸顯現出來.本文基于離散元分析程序EDEM對這種貨物載荷的分析方法展開探索.

EDEM是一款具有高工程應用價值的離散元分析軟件[5]，它可以對用戶自定義形狀的顆粒添加力學參數并設置顆粒與顆粒間及顆粒與幾何體間的接觸關系，同時可仿真顆粒隨預定義幾何體的勻速和變速運動，在降低實驗成本的前提下能夠為實際工程問題提供更好的解決方案.本文采用MSC.Patran建立250 000 DWT VLOC三維有限元模型，提取貨艙的結構有限元網格，通過網格處理程序，為EDEM分析提供網格模型，在EDEM中模擬干散貨的裝載情況，分析貨物載荷產生的堆積壓力場.最后通過基于python語言開發的EDEM與Patran的數據接口，實現離散元壓力場向有限元計算模型的映射.

2 實現方法

選取服務于澳大利亞至中國的250K DWT VLOC為目標船.該船總長315.29 m、型寬57.00 m，型深25.00 m，設計吃水18.00 m，貨艙體積19 544 m3，方形系數0.863，最大貨物密度2.5 t/m3，貨艙結構模型見圖1.

圖1 貨艙結構模型

2.1 計算流程設計

船舶直接計算是一項復雜、高難度的船舶結構安全性評估方法，涉及多項技術的綜合應用[6]，同時載荷形式多樣，為保證不對現有直接計算流程產生影響，離散元分析得到的靜壓力加載過程設置在建立船舶有限元計算模型并完成網格附屬性之后和靜平衡調整之前完成.

離散元分析程序EDEM和有限元分析工具Patran具有完全不同的結構網格和仿真環境.其中EDEM僅允許三角形網格，承受載荷并限制顆粒運動，而Patran建立的有限元網格則以體現結構力學特性和傳遞外載荷為主要功能，為提高計算精度，船舶計算中有限元網格以四邊形網格為主，并帶有少量三角形網格.結合貨物壓力及程序間結構網格的功能特點，本文設計形成的靜載荷加載流程見圖2.

圖2 離散元分析流程圖

2.2 有限元網格及網格處理

通過離散元仿真實現靜壓力場計算首先要由目標船的有限元模型提取貨艙模型，該模型為由四邊形和三角形網格組成的五面封閉的網格模型，本文直接采用直接計算中800 mm的網格尺寸即一個肋距，提取得到的貨艙有限元網格模型見圖3，包含了內底板、貨艙底凳前后側板、槽型艙壁、底邊艙斜板、內舷板及艙口圍板.

圖3 有限元分析模型A

在此基礎上，通過網格處理程序對原模型文件進行轉化形成了用于離散元分析的模型數據文件.離散元分析中的結構模型保留了模型A的節點位置信息，建立的適用于EDEM分析的貨艙模型見圖4.該模型具有7 070個節點，14 050個三角形單元.

圖4 離散元分析模型B

2.3 DEM仿真參數設計

DEM仿真參數包括計算模型的選擇、材料參數、接觸關系模型、裝載方式等.其中在分析模型的選擇上，本文采用的分析模型[7]為Hertz-Mindlin(No Slip)模型，見圖5.

圖5 球形顆粒碰撞接觸模型

當量彈性模量滿足：

當量半徑滿足：

切向力Ft由切向重疊量δt和切向剛度St決定，滿足

Ft=-Stδt

式中：G*為當量剪切模量.離散顆粒法向力和切向力具有阻尼分量，其中阻尼系數和恢復系數相關，切向摩擦力遵守庫倫摩擦定律，滾動摩擦力通過接觸獨立定向恒轉矩模型實現.

材料參數和接觸關系模型的設置中，本文根據實船裝載特點及裝載的貨物類型，材料參數根據實際裝載貨物類型進行設置，采用CSR規定的休止角和材料特性設置接觸參數，同時，為提高仿真效率，鋼質圍壁與鐵礦砂的滾動系數設置為0.由此建立表1～2的材料參數及接觸關系模型.

表1 材料參數

表2 接觸關系

在裝載方式上[8]，采用在以艙口中心為原點的6 m邊長的8邊形范圍內隨機生成顆粒并在重力加速度作用下裝載堆積的方式.由此建立的貨艙重載裝載模型見圖6.

圖6 EDEM裝載模型

2.4 壓力場映射

離散元分析中考慮了顆粒的彈塑性及運動、碰撞，在顆粒裝載完畢后通過總體坐標系下三個方向上的合力隨時間的變化曲線可取出相對穩定態下的壓力場分布.裝載量為93.4萬個顆粒的貨艙裝載模型X方向合力隨時間的變化曲線見圖7.

圖7 貨艙X方向合力隨時間的變化曲線

離散元分析中干散貨顆粒與三角形網格在接觸平衡模型可求得三角網格節點上的集中力.由于輸出格式及節點重編號的問題，網格節點上的集中力難以直接用于船體強度的直接計算.利用python語言開發的EDEM與Patran的數據接口，可從EDEM數據庫讀取并完成節點位置及載荷信息的篩選，在將有效信息加工為Patran的命令流字符串后輸出創建PCL腳本文件[9-10]，以python的OS程序模塊建立程序控制關系，使節點位置及壓力場在模型間建立映射關系，實現離散元分析得到的靜載荷在Patran中的程序化快速施加.其中，由于模型A與模型B節點編號的改變，程序對重建節點的位置關系進行搜索判斷，消除舊節點，將重建節點上的載荷映射到原節點所在位置上，保留了原網格的形狀及質量.詳細實現流程見圖8.

圖8 壓力映射程序流程設計

3 建立對比條件

為驗證通過離散元方法開展貨物載荷計算的可行性，本文加載至貨艙的貨物載荷按本文第二部分描述的方法完成，并以本節定義的載荷總量控制模型和載荷輸出形式為標準與CSR法貨物壓力施加方法建立對比條件.

3.1 載荷總量控制模型

圖9 貨物堆積截面參數設置

由此建立堆積上表面的高度曲面方程為

將Hsur代入得到側壓力方程

根據確定的貨物上表面，由此建立控制關系模型為

式中：N為顆粒的裝載數量；P為球形顆粒的孔隙率；ρ為初始設定的貨物密度；Mt為貨物裝載總重；Fz為貨艙重載工況下貨物的垂向總壓力.

3.2 載荷形式

DEM法通過節點力的方式考慮顆粒與結構網格的接觸關系，輸出結構模型各節點上壓力、摩擦力的合力在總體坐標系三個坐標軸方向上的等效分量.CSR將貨物靜載荷分解為兩部分載荷，即垂直于貨艙圍壁的壓力、傾斜板沿板面向下的剪切力.兩種計算方法下貨物靜載荷的體現見圖10.

圖10 兩種計算方法下貨物靜載荷

本文將CSR計算得到的兩部分載荷按總體坐標系三個方向等效分解到節點上，與DEM法得到的節點力進行對比.

4 分析結果

過程中選取半徑0.95，0.85，0.75，0.65，0.55，0.45，0.35，0.25和0.15 m的顆粒按2 000個/s的裝填速度分別開展了裝載實驗，提取得到相對穩定態下作用于整個貨艙垂向合力隨顆粒半徑的變化關系見圖11.

圖11 貨艙Z向合力隨顆粒半徑的變化關系

由圖11可知，當半徑達到0.15 m，網徑比(網格尺寸/半徑)達到5.3時，顆粒與單個網格最大接觸數量達到8.4個，DEM法計算得到的總體載荷與CSR法計算得到的總體載荷已基本一致，差距僅為0.36%.在此條件下形成的重載狀態下CSR法和DEM法貨物靜載荷等效為節點載荷后的壓力云圖見圖12.

圖12 兩種計算方法下貨物靜載荷等效節點力云圖

由圖12可知，貨物顆粒在內底板等構件上出現了明顯的隨機壓力分布，與CSR法中貨物載荷理想化的均勻的載荷分布有較大出入，這與實際礦砂裝載過程中因貨物顆粒間的隨機性碰撞后的堆積行為相符.

在整體上，DEM計算方法在作用于貨艙結構網格上的X,Z方向合力與CSR法差距較小，Y向合力雖然與CSR法有量級上的差距，但與X方向合力同量級，這與粒徑半徑的選擇相關.貨艙重載工況下三個方向的合力表現見表3.

在該離散元仿真參數下，本文提取了半徑在0.15～0.95 m范圍內，9種顆粒半徑裝填完成在內舷板上產生的側壓力隨高度的變化曲線，并與CSR法進行對比，離散元分析方法得到的側向壓力隨裝填顆粒半徑的縮小，在內舷板上的側向壓力逐漸增大，側壓力變化逐漸平滑,見圖13.

表3 兩種計算方法下貨物載荷合力對比

圖13 內舷板Y方向沿高度方向壓力變化

5 結 論

1) 本文開展了離散元與有限元計算模型結構網格間壓力場映射方法研究，建立了基于Python語言和Patran PCL命令流的壓力場映射方法的可行的流程設計，可在不改變原模型節點位置及網格質量的情況下，可通過場函數的形式進行大量節點力的高效加載.

2) 從離散元分析結果看出，使用網徑比達到5.3的分析模型進行散貨裝載仿真還是會在內底板等構件上出現壓力場的概率性分布，這與CSR理想化的規律的載荷分布模型不符.

3) 本文建立了一套基于離散元的貨物壓力計算流程，并以實際貨艙模型為例進行了裝載測試，為貨物載荷的研究提供了思路.