基于ALE的矿粉货物液化晃荡问题并行数值模拟

丁峻宏，金允龙，王 惠（上海超级计算中心，上海003；上海船舶运输科学研究所，上海0035）

基于ALE的矿粉货物液化晃荡问题并行数值模拟

丁峻宏1，金允龙2，王 惠1
（1上海超级计算中心，上海201203；2上海船舶运输科学研究所，上海200135）

含水量较高的矿粉货物在海上运输过程中易出现液化，形成自由表面并使晃荡现象加剧，严重威胁船舶运输安全。针对船载液化矿粉晃动和舱壁冲击问题，采用ALE有限元方法对其进行了细致建模和计算模拟，从三维角度考察了在船舱一定装载率和运动状态下液化矿粉的晃荡现象和特性；同时，借助不同仿真软件，对计算结果的合理性和准确性进行了相互比对和分析。模型求解借助了高性能计算资源，以解决问题求解时间长和多组计算工况带来的大规模计算需求；结合所建计算模型特点和流固耦合特性，研究了多核环境下两种不同区域分解策略和实现方式，通过并行计算性能数据比较分析，以探求更为合理的并行加速策略。

矿粉液化；ALE；晃荡；流固耦合；并行计算；区域分解

0 引 言

随着中国经济社会的快速发展，对矿产资源的需求一直处于高位，许多重要能源和矿物原材料持续表现出高峰需求。近几年来，从国外进口的铁矿和有色金属精矿供应日趋增多，以镍矿为例，一方面有限的国内资源供应难以满足市场需求，另一方面镍又是高品质不锈钢生产过程中的重要添加剂，工业需求催生了国外镍矿的大量进口。然而，镍矿货物在运输过程中，时常由于种种原因导致含水量超标，从而在遇到摇摆和振动时发生水分析出矿体表面的现象（亦称矿粉"液化流动"），形成自由液面，导致船体稳性损失及局部强度超载，严重威胁船舶运输安全[1]。

仅以2010年底到2011年底为例，在近一年时间内，在中国周边海域就有四艘船舶相继由于矿粉货物液化造成船体稳定性不足最终导致船体倾覆沉没，包括“建富星”号、“南远钻石”号、“宏伟”号和“皇后”号，事故的高度频发引发业界对这种货物的运输安全性日趋重视，各国行业组织[2]也开始致力于研究制定更加统一实用的安全规范。

不同矿源地和不同含水量条件下的矿粉物性存在差别，可以在实验室条件下开展相应的物性研究和小比例液化矿粉晃荡试验模拟，但是实验模型力学行为和实际船舶运输条件下相比仍有一定的区别，因此借助数值仿真手段进行液化矿粉晃荡和安全性研究具有其独特和重要的作用。

由于船舱晃荡引起液化矿粉或矿砂的流动，涉及到非线性、大变形和流固耦合等复杂物理特性。过去有不少针对LNG液舱晃荡开展的研究，可以从技术手段上提供借鉴，比如说有的研究采用欧拉法结合多相流来模拟，并选用流体体积法（VOF法）跟踪自由表面[3]；有的研究采用任意拉格朗日—欧拉法（ALE）来求解含自由液面的液体晃荡问题[4]；有的采用无网格的拉格朗日法来模拟，对大变形区域的模拟和自由表面的捕捉也有其独到之处[5]。但是，从现有文献综合来看，目前国内外使用数值仿真方法来开展粘稠泥沙或矿砂流动模拟的先例不多[6-7]；虽然国内近年来也开始有一些从事液化矿物运输安全性方面的研究出现，但多专注于船舶倾覆机理分析和预防监控应急措施[8-9]，鲜见开展液化矿粉船舱晃荡数值模拟的研究案例。

此外，以往类似研究都基于二维分析或者小比例研究模型，如若拟实建立三维船舱模型和细致的空间网格，则可以对液化矿粉船舱晃荡过程进行更为立体真实的仿真和考察。但大规模的流固耦合计算和深入比较分析，却并非一般计算机的计算能力能够企及。

本文针对某型运输船舶中一定装载率下的液化矿粉，建立了大规模全三维流固耦合有限元模型，计算则基于ALE有限元显式算法，并借助高性能计算机和多核并行求解工具LS-DYNA软件，研究了液化矿粉在晃动过程中的力学行为和周期性变化规律，并与其他仿真结果进行了对比。

1 原理与方法

1.1 ALE算法和VOF法

对于呈液化流动状的泥沙、矿砂或矿粉计算模型来说，拉格朗日算法将网格处理为随材料一起运动，界面跟踪和边界设置比较容易，但很难处理由此引起的网格大变形和大位移；欧拉算法将材料处理为从固定网格中穿越，避免了网格畸变的同时却不利于追踪界面和边界，难点还有多相流的模拟[10]。ALE算法系统计算所有单元的质量、动量和能量输运量、单元密度和速度等参量，并随时间步进行更新，最后采用更新后的密度和状态方程中单位内能来计算单元中的压力值。

ALE算法在拉格朗日和欧拉两种坐标系基础上，又引入任意参考坐标系，使得物质导数描述为：

式中：Xi是拉格朗日坐标，xi是欧拉坐标，wi是物质速度v和网格速度u之间形成的相对速度差。

公式（2）至公式（4）表述了ALE算法的守恒控制方程：

a.质量守恒方程

b.动量守恒方程

c.动量守恒方程

式中：ρ为流体密度，σi j为应力张量，b为单位质量的体积力，E为单位质量的总能量。

有限差分法被用于实现每个节点任意时刻的速度u和位移x更新：

式中：M为对角质量矩阵，Fint为内力矢量，Fext为与体积力和边界条件等有关的外力矢量。

VOF法通过研究网格单元中流体和网格体积比函数Vr来确定自由面，追踪流体的变化，被广泛用于求解流体和固体力学中的很多问题。在晃荡计算模型求解时，Vr取值为1或0时分别代表全部为指定相流体所占据单元和无指定相流体单元两种状态，部分填充单元（0＜Vr＜1）构成了自由液面，界面追踪算法在网格更新和物质输送之前执行。

1.2 多CPU/核并行求解

流固耦合工程问题非线性特点强，计算规模大而耗时久，最好能借助并行计算手段[11]。目前有限元程序并行求解时，多采用大规模并行处理（MPP）模式。如图1所示，一般通过区域分解方法（Domain decomposition method，简称DDM）来处理分析对象，即把整个模拟问题区域划分成很多相对小的求解子区域,把每个子区域分配给不同的CPU处理器（或核心，Core）分别进行求解，区域分解允许每个核单独求解一部分问题，而核之间通过交互机制进行数据交换，最后将所有子区域的解汇总得到整个区域的全局解。

图1 并行求解的模型分区和CPU布局Fig.1 Model of domain decomposition for parallel computing and CPU distributing

在对区域分解并行求解的计算效果进行比较和评价时，常见的考核性能评价指标包括模型计算总耗费时间TC，加速比SN和并行效率EN。其中SN和EN分别用来考量并行加速处理效果和并行计算时整体系统的资源利用率，即：

2 应用实例

2.1 计算模型

本文开展液化矿粉船舱晃荡数值模拟研究的基本计算模型见图2（隐去部分舱壳单元），图3显示的是后续为开展压力时程对比而选定的舱壁中下部监测点位置。整体模型基本尺寸为29.5 m（长）× 32.3 m（宽）×18.4 m（高），所装载的矿粉高度距内底板8.0 m，矿粉液化后假定上层还出现了深0.3 m的稀泥浆，船舱重心距内底板5.3 m。矿粉、泥浆和空气各对象有限元模型均采用ALE六面体体网格划分，其中泥浆网格还进行了局部加密；船舱采用四面体壳单元来模拟，整个计算模型单元和节点总数均在24万左右。

图2 矿粉货物晃荡计算有限元模型Fig.2 FEM model for sloshing of ore fines cargo

图3 舱壁压力监测单元位置示意图Fig.3 Diagram of pressure monitoring elements at bulkhead

各对象的基本物理属性见表1。其中船舱舱壳用MAT_RIGID刚性材料模型来模拟，而矿粉、泥浆和空气三者都采用MAT_NULL材料模型，同时使用EOS_GRUNEISEN状态方程[12]定义可压缩流体材料的压力：

式中：p为介质压力，ρ和ρ0分别为当前和初始密度，μ=ρ/ρ0-1，C为冲击波速，γ0为GRUNEISEN常数，S1、S2和S3是vs-vp曲线的斜率系数，a是压缩时γ0的一阶体积修正，E为内能。以本文液化矿粉为例：可取波速C=1 647 m/s，S1=1.921，S2=-0.096，S3= 0，a=0，γ0=0.35。

在装载液化矿粉船舱晃荡的模拟过程中，设定船舱绕重心横向往复摆动，最大摇摆角度为10°，取晃荡激励周期T为10 s，每个工况的计算时间取为3T。

表1 各建模对象的基本材料参数表Tab.1 Basic material parameters of each modeling object

2.2 并行求解

本文建立的计算模型规模不小，每个工况问题物理时间久，而且ALE算法求解时间也比拉格朗日算法要长得多。如果使用普通工作站（4CPU Cores；内存8～16 GB），计算一个工况至少需要近7天，如果计算工况多的话，则对计算效率、磁盘空间直至研究的顺利开展都提出了挑战，因此必须要借助高性能计算资源。上海超级计算中心搭建有专门针对工业和工程领域应用需求的“蜂鸟”集群，该集群包括65台HS23刀片计算节点，每个节点含两个Intel八核Intel处理器（主频2.6 GHz,共享64 G内存）。

利用LS-DYNA中ALE算法进行问题求解时，采用坐标递归对分方法（Recursive Coordinate Bisection，RCB）来进行计算模型的多核区域分解处理。该方法将网格沿垂直于模型最长坐标方向将模型对分，将相邻边界节点复制分布到两侧区域，持续递归分割直至满足分割收敛条件。通常采用程序推荐的基于ALE单元均分策略来组织进行计算区域的剖分方案A（见图4（1））；进一步考虑到船舱尺寸特征、晃荡方向性特点以及流固耦合界面的均衡分配需求，本文在此基础上又提出了新的子区域剖分策略方案B（见图4（2）），并就这两种不同区域分解方案对应并行加速效果进行测试对比。

图4 不同的区域分解策略及分区结果Fig.4 Different domain decomposition strategies and results

3 计算分析

3.1 船载液化矿粉晃荡分析

一般摇摆船舱内液体的晃动可以认为是驻波、行波、水跃和涡旋等形式的组合，如何从三维的角度去再现舱内装载矿粉液化后随船摇摆下的晃荡现象非常重要，由于篇幅关系，本文在此仅给出了几个特定时刻舱内装载矿粉晃荡时自由表面变化和密度分布云图（见图5）。可以发现，与一般舱载均质低密度低粘度液体晃荡现象和特点有所不同，绝大部分液化矿粉的晃动形式偏向驻波，而其上层部分矿粉和表层稀泥浆在两侧舱壁之间往复流动现象更似行波，从三维视角来看自由表面流动平缓且不规则（见图6）。另外，液化矿粉在舱内晃荡流动时，与船舶本身横摇运动之间有一定相位差，该现象在船舶从摇摆至两侧最大倾角继而向平衡位置运动时最为明显，若在海上持续受到不利外载作用并导致矿粉逐渐向一侧倾斜聚集，终将影响船舶的运输安全。

图5 不同时刻矿粉货物的自由表面形状Fig.5 Free surface shapes of ore fines cargo at different time

图6 矿粉货物的流体密度云图（t=25 s）Fig.6 Fluid density contour plot of ore fines cargo(t=25 s)

3.2 仿真数据对比和应急措施验证

船载矿粉发生液化产生晃荡时多发生在海上运输途中，一般很难直接获取来自现场的数据和资料。因此仿真结果数据的对比，只能借助于相似比例模型试验或者其他仿真算法。将本文利用LS-DYNA软件计算得到的舱壁监测点（见图3）压力时程曲线，与利用CFD软件FLUENT计算得到的结果进行了对比，见图7。两组结果数据初始时刻压力值有所不同是因为FLUENT模拟时考虑了静压，但曲线后续变化趋势表明随着船舶开始晃荡后静压的影响逐渐降低，曲线整体周期性规律和最大峰值基本趋于一致，但相位和幅值略有差异，这是受不同流体模型参数设置以及不同软件具体算法的影响。初步计算表明两种软件所得结果还是比较接近的，这可为借助不同方法开展晃荡仿真数值模拟时的合理性和准确性评判提供相互验证。

图7 不同模拟软件计算得到的压力时程对比Fig.7 Comparison of pressure time history due to different simulation softwares

图8 安装隔板后矿粉货物的流体密度云图（t=25 s）Fig.8 Fluid density contour plot of l ore fines cargo after installing middle baffle(t=25 s)

在易液化矿粉装船以及航行途中发生横倾时，现场都有一些规范和应急举措，基于现有三维仿真模型，还可以从中选取一些建模方便的方案开展验证和评估，充分发挥数值计算模型的优势。如图8所示，有种措施是在舱内中部用木板纵向埋于矿粉表面并与舱前后两侧固定，与图6相比，计算表明该措施的确可以降低晃荡幅度，从而减小自由液面的影响。

3.3 并行求解性能分析

以船载液化矿粉晃荡计算模型计算三个摇摆周期（30 s）为测试对象，从图9（1）来看，1核和2核由于两种区域分区结果无异，因此计算耗时也一致。4核时方案B开始略有不同，使用8核和16核时该方案节省计算时间现象最为明显，到32核和64核时能继续保持其比较优势，表明其仍具有一定的可扩展性。从计算耗时具体数值来看，单核计算需要25天，双核也要13天，不只是计算时间上令人难以接受，更重要的是对软硬件计算环境工作稳定性也提出了很高要求，一般很难绝对保证问题能在预期时间顺利算完；使用16核后，计算耗时可显著降低到3天之内，给多工况求解和深入研究提供了时间上的可能性。

图9 不同区域分解策略下并行性能数据对比和分析Fig.9 Parallel performance data comparison and analysis with different decomposition strategies

进一步从图9（2）分析来看，两种分区方案对应的加速比效果均随核数增加呈非线性上升趋势，其中方案B自4核之后始终保持领先并略呈扩大趋势，表明其负载均衡得到了进一步改善。但从并行效率数据来看，32核之后两种方案并行效率都已开始回落至50%甚至更低，从计算资源的有效配置和使用角度来看，不建议再使用过多的计算核心数。

因此对于本文以及类似同等规模计算模型而言，求解时使用16～32核应该是比较合理的选择，单工况求解时所耗时间将会大大降低，2～3天就可完成计算；而在此基础上若采用方案B执行计算模型的区域分解，则可将单工况总计算耗时在常用方案A的基础上再进一步缩减3～10小时。

4 结 论

本文详细阐述了如何利用高性能计算软硬件资源，与流固耦合算法以及并行计算技术紧密结合，就船舶远洋运输过程中矿粉液化后产生的晃动冲击问题开展深入研究，给出了具体实现方法和相关参数设置，得到了一些有意义的计算结果并对此进行了深入分析，最后有以下认识：

（1）结合ALE算法和VOF法的三维有限元方法，可有效捕捉自由液面的变化，模拟液化矿粉受船舶摇摆引起的往复晃荡和舱壁冲击等力学行为，验证减小晃荡的实际措施，并为舱壁结构承压分析打下基础；不同数值模拟软件计算结果可以相互进行比对，以间接验证该问题仿真方法的可行性和正确性。

（2）以高性能计算资源为支撑，可对大规模船载液化矿粉晃荡冲击问题进行细致建模和深入分析，并将问题求解时间的影响控制在有利范围内；开展并行计算时，应该根据求解问题类型和模型力学特点选用合理的区域分解方案和合适的多核并行加速策略，以实现计算资源的优化配置和最佳利用。

[1]雷 海.红土镍矿粉的安全运输[J].航海技术,2011(1):27-28. Lei Hai.Transprotation safety of laterite nickle ore fines[J].Marine technology,2011(1):27-28.

[2]Class N K.Guidelines for the safe carriage of nickel ore(Second edition)[M].Tokyo:Nippon Kaiji Kyokai,2012.

[3]刘桢兵.基于VOF法的液舱晃荡数值模拟及载荷计算[J].重庆理工大学学报（自然科学）,2011,25(3):24-29. Liu Zhenbing.Numerical simulation of liquid tank sloshing and pressure calculation based on volume of fluid method[J]. Journal of Chongqing University of Technology(Natural Science),2011,25(3):24-29.

[4]管延敏,叶恒奎,陈庆任.基于ALE有限元法二维液体晃荡的数值模拟[J].船舶力学,2010,14(10):1094-1099. Guan Yanmin,Ye Hengkui,Chen Qingren.Numerical simulation of 2D sloshing with Arbitrary Lagrangian-Eulerian finite element method[J].Journal of Ship Mechanics,2010,14(10):1094-1099.

[5]陈正云,朱仁庆,祁江涛.基于SPH法的二维液体大幅晃荡数值模拟[J].船海工程,2008,37(2):44-47. Chen Zhengyun,Zhu Renqing,Qi Jiangtao.Numerical simulation of sloshing in two dimensional liquid tank based on SPH method[J].Ship&Ocean Engineering,2008,37(2):44-47.

[6]Yim S C.3-D Wave-structure interaction with coastal sediments-A multi-physics/multi-solution-techniques approach[R]. Corvallis:Oregon State University,2008:1-10.

[7]Karmakar S,Kushwaha R L.Simulation of soil deformation around a tillage tool using computational fluid dynamics[J]. Transactions of the ASAE,2005,48(3):923-932.

[8]王 威.浅谈海运精矿粉的安全运输[J].南通航运职业技术学院学报,2009,8(2):63-65. Wang Wei.On safe transport of ore concentrate powder[J].Journal of Nantong Vocational&Technical Shipping College, 2009,8(2):63-65.

[9]赵月林,孟淑娴.易流态化货物安全运输研究[J].大连海事大学学报,2012,38(4):15-18. Zhao Yuelin,Meng Shuxian.Safety transportation of easily fluidized cargoes[J].Journal of Dalian Maritime University,2012, 38(4):15-18.

[10]Aquelet N,Souli M,Gabrys J,et al.A new ALE formulation for sloshing analysis[J].Structural Engineering and Mechanics,2003,16(4):423-440.

[11]李 政,金先龙,亓文果.流体—结构耦合问题的有限元并行计算研究[J].计算力学学报,2007,24(6):727-732. Li Zheng,Jin Xianlong,Qi Wenguo.Study on the parallel algorithm for finite element analysis of fluid-structure interaction[J].Chinese Journal of Computational Mechanics,2007,24(6):727-732.

[12]LSTC.LS-DYNA keyword user’s manual[M].California:Livermore Software Technology Corporation,2012.

ALE-based parallel numerical simulation for sloshing problem of liquefied ore fines cargo

DING Jun-hong1,JIN Yun-long2,WANG Hui1
(1 Shanghai Supercomputer Center,Shanghai 201203,China;2 Shanghai Ship and Shipping Research Institute,Shanghai 200135,China)

Ore fines cargo with high moisture content is liable to liquefaction,making free surface emerging and intensifying sloshing behavior,and finally pose a threat to marine transportation security.With the benefit of Arbitrary Lagrangian-Eulerian(ALE)finite element method,detailed modeling and refined simulation were implemented to deal with shipping liquefied ore fines sloshing problem,and to investigate three-dimensional sloshing phenomena and characteristics with certain charging ratio and motion state. Two different numerical nodes were used to make comparison and analysis for rationality and accuracy of results.High-performance computing(HPC)resource was utilized to meet the challenge of large scale computing requirement due to time-consuming solving and load cases.In view of model size and fluid-structure intersection property,two domain decomposition plans were proposed to better improve load balance for multicore processors,and compared with parallel performance data to determine more reasonable parallel acceleration strategy.

liquefaction of ore fines;ALE;sloshing;fluid-structure interaction;parallel computing; domain decomposition

TP391.9 U661.3 U674.13+4.1

A

10.3969/j.issn.1007-7294.2015.08.006

1007-7294（2015）08-0927-07

2014-12-18

国家高技术研究发展计划（863计划）课题（2012AA01A308）；国家重点基础研究发展计划资助课题（2012CB723804）

丁峻宏（1977-），男，工学博士，高级工程师，E-mail:Jhding@ssc.net.cn；金允龙（1964-），男，研究员。