基于离散元的矿砂船货物载荷计算方法

宋喜庆 张少雄 胡丰梁

(武汉理工大学交通学院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) 本文建立了一套基于离散元的货物压力计算流程，并以实际货舱模型为例进行了装载测试，为货物载荷的研究提供了思路.