液舱晃荡载荷计算与强度校核系统研究

赵 吉，许一敏，程 成，冷文浩

(1.无锡环境科学与工程研究中心，江苏 无锡 214000;2.中国船舶科学研究中心，江苏 无锡 214082;3.江南大学 物联网工程学院，江苏 无锡 214122)

近年来，随着天然气消费量的急剧增长，大型海上运输工具薄膜型LNG 运输船的建造数量和运输频率也在显著增加。目前，该类型船的设计与制造已成为造船领域内的关注热点。当LNG 船在波浪中的运动频率与液舱内液体固有振动频率相近时，液体运动非常剧烈，将会对液舱产生强烈的冲击，严重时还会导致结构的失效和破坏，酿成重大事故，不仅造成严重的生命和财产损失，而且还会导致海洋环境的污染。因此，晃荡载荷已成为LNG 船安全性研究的重要内容之一［1］。

目前关于液舱晃荡的研究可以归结为理论研究［2］、实验研究［3］和数值计算［4］三个方面。随着研究的深入，发现理论研究存在很大的限制，实验研究也遇到很多的困难，而随着计算机技术的高速发展，数值计算已成为现代晃荡研究的重要方法。而且数值计算还具有花费少、速度快、重复性好、风险小等优势。

由于计算液舱晃荡载荷的方法和步骤非常复杂，完全使用商业化的求解器，比如Dytran 等来计算晃荡载荷需要耗费更大量的时间;使用手工计算和加载，需要耗费大量的工作量，而且容易出错。本文针对薄膜型LNG 船，结合中国船级社规范中关于晃荡要求和晃荡载荷的计算部分，整理晃荡载荷计算公式、屈服校核标准，利用MSC.Patran［5］专门提供的二次开发语言PCL (PATRAN Command Language)，设计程序界面［6］，开发了一套集舱室识别、晃荡载荷自动计算与施加、屈服校核、校核报告于一体的系统，并已经开始被船级社应用在实船的校核工作中，如14.7 万方薄膜型LNG 船。

1 流程说明

晃荡载荷计算与校核的过程可以概括为四个步骤，如图1 所示。

1)有限元前处理。实现晃荡载荷计算参数的输入和保存，采用舱室识别功能自动识别并保存载荷施加的作用域，施加边界条件。

2)晃荡载荷计算。读取载荷计算参数进行计算，并施加到相应舱室中生成晃荡载荷，最后与边界条件组合生成晃荡工况。每种工况类型对应多个不同的装载高度，每个装载高度又对应多个水平的晃荡载荷，需要循环计算生成相应的晃荡工况。

3)结构强度计算。调用MSC.Nastran 来计算在施加了载荷后船体结构的响应，并将计算的结果导入到有限元模型中。

4)有限元后处理。读取计算的应力值和构件材料属性，根据构件的校核标准进行屈服校核，并用云图等形式显示校核结果，最后将各个构件的校核结果自动生成到校核报告中。

图1 晃荡载荷计算与校核流程Fig.1 Sloshing loads calculation and assessment process

2 舱室识别

晃荡载荷施加的作用域是舱室，然而手工划分舱室非常复杂，工作量比较大，准确性也不高。舱室识别功能能够根据指定的单元自动搜索相关单元，直至搜索到相应舱室的所有单元，这不仅大大减轻工程人员工作量，也保证了舱室的正确性。本文提出的系统提供了两种舱室识别方法，分别为拓扑算法和切分拼接算法。

2.1 拓扑算法

拓扑算法根据构成舱室的相关构件，通过公共节点、公共边(edges)和网格单元之间的拓扑关系，搜索相临单元，不断循环搜索，直至完成一个舱室的识别。以薄膜型LNG 船液货舱为例，是由内壳板、内底板、内凸型甲板等内部构件组成。搜索舱室时，先要确定舱室的大致范围与起始单元，然后从起始单元出发，找出所有与之相邻的单元，过滤掉舱室范围之外的单元后，将搜索得到的单元保存到舱室单元中，接下来把新搜索到的单元作为起始单元，如此反复，直至找到所有的舱室单元，如图2 所示。

拓扑算法需要在配置文件中定义所有与类型舱室相关联的船体构件，强烈依赖于配置文件，所以此方式不具有一般性和通用性，在船艏或艉部的舱室识别时维护配置文件将耗费极大的工作量，且容易出错。

图2 拓扑算法Fig.2 Topology algorithm

2.2 切分拼接算法

切分拼接算法根据三维立体的几何性质，任何两个相交的面可以确定第三个面，其中第三个面必须满足分别与这两个面相交且不交于同一交线。用该方法搜索舱室先要对整个模型进行预处理，将模型中相交的面都切分为单独的面，然后以这些面为基准，拼成想要的舱室。

程序流程为:首先根据构件之间的相交关系对模型进行切分，将模型切分为若干块没有被分割的面，即初始化。如图3 所示，其中被切分的不同的面被分别标注出来。为了提高程序的效率，保证Patran 的运行稳定性，模型切分产生的“面”输出到数据文件中，并对每个面加上编号，由程序管理。

根据三维立体的几何性质，任何两个相交(有公共edge)的面可以确定第三个面，其中第三个面必须满足如下特征:分别与这两个面相交且不交于同一交线，也即第三个面的两条边分别和已定的两个面是公共边。如图4 所示，面a 和面b 相交，面c 分别与面a 和面b 相交且不交于同一交线，由此得出面c 是属于该三维舱室空间的边界。图4 为实际模型中的三个相关面的拓扑关系。

采用上述的算法原理，根据指定的任一属于舱室空间边界的单元作为起始单元，将切分出的面再拼接成舱室空间。对于选定的起始单元，若起始单元所属的“面”属于唯一的舱室，则识别出的舱室空间较为准确，若起始单元所属的“面”属于两个舱室，则可以再选择相邻相交“面”上的一个单元一起作为起始单元。切分拼接算法可以实现舱室的自动识别，且方法的有效性也得到了验证，是一种可行的技术方案，可大大降低验船师的工作量，缩短审图周期。

图3 模型切分效果Fig.3 Model segmentation effect

图4 拼接示意Fig.4 Stitching schematic diagram

3 晃荡载荷

3.1 晃荡运动水平

液舱的晃荡水平分为以下三种:

水平一:静态或准静态过程，该水平下晃荡运动引起的载荷以静载荷为主，动态运动载荷可忽略。

水平二:晃荡运动的动态不放大过程，该水平下晃荡运动引起的载荷包括静载荷和动载荷，但动载荷还未达到冲击载荷形式。

水平三:晃荡运动的动态放大过程，该水平下晃荡运动引起的载荷包括静载荷和动载荷，其中动载荷是主要载荷，达到冲击载荷的作用形式。

3.2 水平一准静态晃荡载荷

水平一准静态晃荡载荷［7］是由液面线性变化引起的作用于液舱边界上的压力。液货静水压力计算公式为P = ρgh，h 为有限元单元与液面的距离，因此，计算得到晃荡后液面的方程即可计算水平一载荷。图5为计算水平一准静态晃荡载荷流程图。

图5 水平一准静态晃荡载荷计算流程Fig.5 Level-one quasi static sloshing load calculation flow chart

图6 液面纵向运动状态Fig.6 Liquid surface longitudinal motion state

图7 分解示意Fig.7 Decomposition diagram

3.2.1 纵向运动

晃荡运动方向为纵向运动时，需要计算装载液面的总体积V 和临界体积V临。总体积根据舱室的装载高度，结合舱室形状，将舱室划分成长方体和四棱柱进行计算;临界体积划分情况如图6 所示，需将舱室划分成图7 所示的多面体，根据V 与V临的大小关系判断液面线性方程与液舱的交点位置，采用等体积法计算的公式为:

其中:V1为图7(a)类三棱柱的体积，V2为图7(b)类四棱柱的体积。利用积分求得图7(a)类三棱柱体积计算公式为:

其中:a、b 为下底梯形的上底和下底，h1为下底梯形的高，c 为平行于下底梯形的线段，h2为c 与下底梯形的距离。

图7(b)类四棱柱体积计算公式为:

其中:a、b 为下底梯形的上底和下底，h1为下底梯形的高，c、d 为上底梯形的上底和下底，h3为上底梯形的高，h2为上底梯形与下底梯形的距离。

3.2.2 横向运动

晃荡运动为横向运动时，计算方法与纵向运动相似，通过分解计算出总面积与临界面积，计算得到液面线性方程与液舱的交点位置，简化为等面积法来计算舱内液面的方程。

3.3 水平二动态晃荡载荷

水平二动态晃荡载荷是由液货静水压力和晃荡运动引起的作用于液舱边界上的平均晃荡压力，计算公式为:

其中:静水压力P0= ρgh，h 为有限元单元与液面的距离。

平均晃荡压力纵向运动计算公式如下:

横向运动计算如下:

其中:Cl= 0.9x01/［1 + 9* (1 － x0)2］，取值不小于0.25;x0= Tx/Tp;Tp为船体纵摇运动周期;x01=Ct= 0.9y01/［1 +9* (1 － y0)2］，取值不小于0.25;y0= Ty/Tr;Tr为船体横摇运动周期;y01

图8 为水平二动态晃荡载荷的计算流程图。

图8 水平二动态晃荡载荷计算流程Fig.8 Level-two dynamic sloshing load calculation flow chart

3.4 水平三晃荡冲击载荷

水平三中动载荷是主要载荷，为冲击载荷的作用形式，动载荷的冲击作用主要包括顶边舱斜板冲击压力的计算，舱顶冲击压力的计算，液舱下部冲击压力的计算［8］。图9 为水平三晃荡冲击载荷的计算流程。

4 晃荡工况

表1 为研究晃荡载荷典型的4 个运动工况，其中工况LCA_L、LCB_L 为纵摇运动，舱室内液体沿着船长方向运动，引起纵向的液面变化，所以在计算晃荡载荷时主要考虑纵向的变化;而工况LCA_T、LCB_T 为横摇运动，与LCA_L、LCB_L 情况相反，需要考虑船宽方向的液面变化。

4.1 装载高度选取

根据规范的要求，每种运动工况至少要对液舱5%h，10%h，70%h，95%h 这4 个装载高度进行计算。对于每个运动工况来说，在特定装载高度下，还需要根据液舱尺寸、装载液面宽度、液舱频率与船体固有频率比值的不同，施加上述一种或多种水平载荷，并生成对应水平的载荷工况。需要校核所有这些工况，来查看构件的强度是否符合规范要求。

图9 水平三晃荡冲击载荷计算流程Fig.9 Level-three sloshing impact load calculation flow chart

表1 晃荡载荷工况表Tab.1 Sloshing load conditions

4.2 设计思路

从晃荡载荷工况表中，整理计算晃荡载荷垂荡、纵摇、横摇等参数，在Patran 的基础上利用二次开发语言PCL，编写晃荡工况计算的界面，实现用户对晃荡载荷计算参数的输入和修改，程序获取晃荡参数后，根据图10 晃荡工况载荷计算的总体流程，判断晃荡工况所需施加哪几个晃荡水平载荷，按照图5、图8、图9 的计算流程自动计算晃荡载荷，并生成对应的晃荡工况。图11 为本系统晃荡载荷计算前的工况设定界面。

图10 晃荡工况载荷计算流程Fig.10 Sloshing conditions load calculation flow chart

5 有限元后处理

屈服强度校核是判别船体构件在载荷作用下是否满足规范要求的一个标准。本系统通过MSC.Nastran 计算构件结构响应的结果并导入到Patran 后，通过PCL 中的函数res_utl_extract_elem_results2 获取指定的板单元相当应力(σe)和杆单元轴向应力或梁单元正应力(σrod)，与材料屈服极限(Reh)进行相比，即得到屈服利用因子:

规范中对每个构件的屈服利用因子都进行了规定，如表2所示，当屈服利用因子超出要求时表示不合格［8，9］。

图11 晃荡工况计算系数设定Fig.11 Sloshing conditions calculation factor setting

表2 屈服利用因子Tab.2 The yield using factor

6 算 例

本文以14.7 万方薄膜型LNG 船体三舱段模型为例，来验证本系统的有效性。边界条件为后端中心线中和轴处的独立点沿船纵向轴的线位移约束，施加总体垂向和水平弯矩，所有端面纵向单元节点与独立点刚性连接，所有甲板、内底板和外板施加横向弹簧，舷侧、内壳施加垂向弹簧。工况为晃荡载荷工况LCA_L，装载高度为70%h。由于液舱的尺寸、装载液面宽度、液舱频率与船体固有频率的关系，该工况、该装载高度下只需考虑晃荡水平一、晃荡水平二载荷，见图12 晃荡载荷云图。

为了验证实现系统计算结果的正确性，在Mathcad中编写了载荷计算过程，在舱室周界上任意抽样了几个位置上的单元，将程序计算结果与Mathcad 计算结果进行对比，如表3 所示。从表中可以看出本系统计算结果与Mathcad 对比误差率都为0，说明本系统计算结果的准确性。

图12 液舱晃荡载荷云图Fig.12 Sloshing load cloud picture of liquid Tank

本系统根据规范规定的屈服校核准则开发了后处理功能，可以方便、快捷、高效地计算和显示船体各结构的屈服利用因子，界面如图13 所示。图14 为所有结构的屈服利用因子云图，屈服极限为235 N/mm2。图15 为所有结构的屈服利用因子与许用屈服利用因子的比值云图，校核结果表明当前的设计可以满足规范规定的要求。为了方便查看和使用校核结果，本系统还提供了将船体信息、晃荡工况的载荷和构件的校核结果输出到计算报告的功能。

表3 晃荡载荷程序测试对比Tab.3 Sloshing load program test

图13 利用因子云图Fig.13 Using factor cloud picture

图14 利用因子比值云图Fig.14 The ratio of using factor cloud picture

7 结 语

对于工程人员来说，舱室划分、晃荡载荷计算和施加、屈服利用因子的计算都非常复杂，手工操作与计算要求非常高，本文结合中国船级社规范内容，在MSC.patran 软件的基础上，首次提出并开发了一套薄膜型LNG 船的晃荡载荷计算与校核系统，实现参数输入、搜索晃荡载荷的作用域、晃荡载荷自动计算与施加、计算并导入结构响应结果、构件屈服校核和校核报告自动生成功能。针对各个水平的晃荡载荷，对计算流程进行了详细的描述，给出了相应的计算公式，保证载荷施加的正确性。该系统可以大大降低工程人员的工作量，已经被应用于实船审图工作中。

［1］刘新立.船载液体晃荡载荷特性研究［D］.武汉:武汉理工大学，2009.(LIU Xinli.The sloshing characteristics research of the liquid in ship tank［D］.Wuhan:Wuhan University of Technology，2009.(in Chinese))

［2］蔡忠华.液货船液舱晃荡问题研究［D］.上海:上海交通大学，2012.(CAI Zhonghua.Study on the sloshing problems of liquid cargo tanks［D］.Shanghai:Shanghai Jiao Tong University，2012.(in Chinese))

［3］蒋梅荣，任 冰，温鸿杰，等.弹性液舱内液体晃荡实验研究［J］.海洋工程，2013，31(5):1-10.(JIANG Meirong，REN Bing，WEN Hongjie，et al.Experimental study on sloshing in elastic tanks［J］.The Ocean Engineering，2013，31(5):1-10.(in Chinese))

［4］赵炳雄，薛鸿祥，唐文勇.LNG 船液舱晃荡载荷滤波算法研究［J］.船舶工程，2013，35(4):15-18，39.(ZHAO Bingxiong，XUE Hongxiang，TANG Wenyong.Study on filter method of sloshing load for LNG Tank［J］.Ship Engineering，2013，35(4):15-18，39.(in Chinese))

［5］刘兵山，黄聪.Patran 从入门到精通［M］.北京:中国水利水电出版社，2003.(LIU Bingshang，HUANG Cong.Patran from novice to professional［M］.Beijing:Water Power Press，2003.(in Chinese))

［6］张军彦，李昌华，李晓辉，基于MSC.Patran 界面的PCL 二次开发［J］.现代电子技术，2010，16:17-19，22.(ZHANG Junyan，LI Changhua，LI Xiaohui.PCL secondary development based on MSC Patran interface ［J］.Modern Electronics Technique，2010，16:17-19，22(in Chinese))

［7］中国船级社.薄膜型液化天然气运输船检验指南［M］.北京:人民交通出版社，2011.(China Classification Society.Guidelines for survey of membrane tank lng carriers［M］.Beijing:China Communications Press，2011.(in Chinese))

［8］中国船级社.钢质海船入级规范［M］.北京:人民交通出版社，2012.(China Classification Society.Rules for classification of sea-going steel ships［M］.Beijing:China Communications Press，2012.(in Chinese))

［9］伍亮.CSR 散货船结构强度直接计算研究［D］.武汉:武汉理工大学，2012.(WU Liang.Research for direct calculation of bulk carrier’s structural strength based on CSR［D］.Wuhan:Wuhan University of Technology，2009.(in Chinese))