基于波面修正设计波的薄膜型LNG船整船屈服强度校核研究

2022-05-31韦喜忠陈鲁愚许一敏

船舶力学 2022年5期

程成，韦喜忠，陈鲁愚，赵吉，郑刚，许一敏

（1.中国船舶科学研究中心，江苏无锡 214082；2.无锡环境科学与工程研究中心，江苏无锡 214063；

3.中国船级社规范与技术中心，上海 200135）

0 引言

随着船舶向着超大型化发展和计算机性能的提升，船舶进行整船有限元结构强度校核的需求也越来越强烈，特别是薄膜型液化天然气（Liquefied Natural Gas，LNG）船，由于装载货物的特殊性，其船体结构需要更高的可靠性和安全性，仅仅采用舱段结构进行强度校核已逐渐不能满足要求，所以各大船级社已经开始将薄膜型LNG 船的全船有限元结构强度校核作为入级的基本要求[1-3]。中国船级社（简称CCS）在《薄膜型液化天然气运输船检验指南》（2015 版）[4]中正式加入了整船直接计算强度评估的要求，并在2018年合并到《钢质海船入级规范》（2018版）[5]（简称“钢规”）。“钢规”中规定整船分析基于动态载荷法开展，波浪载荷使用水动力直接计算得到的实际装载工况下的载荷，以尽可能真实地反映船舶在海上航行时遭遇到的各种工况。同时要求当采用线性水动力分析时，为了更真实地反映水线附近的舷外海水压力，需要进行波面修正。

近些年国内针对各类舰船整船强度校核开展了深入的研究[6-10]，沪东中华造船（集团）有限公司的王佳颖等[11]采用美国船级社（简称ABS）的DLA/SFA 系列软件开展了薄膜型LNG船全船结构屈服和疲劳强度分析，刘敬喜等[12]使用自编的波浪载荷计算程序开展了基于直接计算方法的LNG 船整船强度评估。但是这些研究都不是基于CCS 规范开展的，而且在按照CCS 规范进行校核时存在以下困难：（1）波面修正会带来船体载荷的不平衡问题，载荷调平非常复杂；（2）液货舱部分装载的情况下确定压力计算参考点十分困难，这些问题目前都缺乏有效的解决手段。

本文基于CCS“钢规”中薄膜型LNG 船整船直接计算的要求，对船体载荷自动调平方法和液舱部分装载下的载荷计算方法进行了研究，开发了一套集与挪威船级社（简称DNV）的波浪载荷计算软件Wadam 接口工具、波浪载荷检查、液舱载荷计算、波浪载荷修正、载荷调平、应力计算和屈服校核于一体的系统，并已经开始应用在实船的校核工作中。

1 校核流程说明

根据CCS“钢规”[5]的要求，薄膜型LNG 船整船屈服强度校核的过程可以概括为模型处理、载荷计算和强度校核三个主要步骤，如图1所示，详细说明如下：

图1 薄膜型LNG船整船屈服强度校核流程Fig.1 Checking process of overall membrane-type LNG carrier yield strength

（1）模型处理。在全船有限元结构模型的基础上，确定船体的湿表面和水动力载荷的作用域，并根据当前装载工况的情况使用质量单元模拟重量分布。通过Patran 和Wadam 的接口输出水动力模型、质量模型和结构模型。

（2）载荷计算。首先使用导出的模型文件在Wadam 中进行波浪载荷预报，并确定等效设计波的参数。将载荷预报结果导入到Patran 中，并正确地施加到作用域。然后根据当前的装载设定和等效设计波的参数计算船体外部载荷和舱内载荷，并创建校核工况。最后对工况的外部载荷进行修正并进行载荷调平。

（3）强度校核。调用Nastran 来计算船体结构的响应，并将计算结果导入到有限元模型中。读取计算的应力值和构件的材料属性，根据构件的校核衡准进行屈服校核，并用云图等形式显示校核结果。最后将各个构件的校核结果自动生成到校核报告中。

2 波浪载荷预报

文中线性波浪载荷的预报和长期预报在Wadam软件中进行，而屈服强度的校核在Patran中进行，因此必须开发Patran 与Wadam 软件的接口，一方面将Patran 中建立的模型信息导入到Wadam 中进行载荷预报，另一方面将Wadam预报的波浪载荷和惯性载荷导入到Patran模型中，并正确地施加到作用域上。

2.1 输出Wadam软件的模型文件

使用Wadam 软件进行波浪载荷预报需要输入水动力模型、结构模型和质量模型，随着Wadam 计算能力的提升可以直接将施加了水动力载荷的整船结构模型的外壳输出为水动力模型，进而避免了预报的水动力载荷需要插值映射到结构模型网格上的问题。同时整船模型中使用了质量单元来模拟货物质量，输出的结构模型中包括了所有的单元，也就反映了装载工况下的实际全船质量，因此可以直接使用质量模型。

通过指定一个船体单元搜索其相连单元的方法，可以实现船体水动力模型的确定，如图2 所示；通过对Wadam 软件后缀为.T的输入文件（简称T文件）格式的研究，开发了模型导出功能，实现Patran与Wadam软件的数据交互。

图2 薄膜型LNG船水动力模型Fig.2 Hydrodynamic model of membrane-type LNG carrier

2.2 导入Wadam软件的预报载荷

Wadam 软件预报的工况信息可以输出到.S 为后缀的文件（简称S 文件）中，波浪载荷和惯性载荷可以输出到.L 为后缀的文件（简称L 文件）中，并且在载荷输出时也输出了工况信息，因此通过对S 文件和L文件格式的研究，可以读取载荷数据和工况数据并创建到Patran中。

载荷创建到Patran 后可以通过直接查看每个工况实部和虚部结果的方式进行检查，如图3 所示，同时也可以通过定义剖面位置然后对剖面处进行载荷积分的方式来检查，如图4～5所示。

图3 查看导入的波浪载荷Fig.3 View of wave loading

图4 船体剖面设定Fig.4 Setting of ship hull section

图5 船体剖面载荷积分结果Fig.5 Results of loading integration of each hull section

3 舱内载荷计算

在进行舱内载荷计算前必须要先识别出舱室周界，如图6所示，并设定舱室形心位置等参数[11]。波浪载荷预报软件预报出了全局坐标系原点处的惯性载荷，在进行舱内载荷计算时需要根据预报的惯性载荷计算出每个舱室形心处的加速度再计算舱内载荷。另外，计算每个液货舱舱内载荷的一个难点是计算出载荷计算点到压力计算参考点的距离，本文根据不同的装载状态采用了不同的计算方法。

图6 液货舱示意图Fig.6 Diagram of liquid cargo tank

3.1 计算详细步骤

每个液货舱每个工况下的舱内载荷计算过程如图7所示，详细说明如下：

图7 液货舱舱内载荷计算流程Fig.7 Calculation process of liquid cargo tank internal loading

（3）计算横摇角φ和纵摇角ψ，

（4）计算由于船舶横摇引起的重力加速度在全局坐标系下的横向分量gT和由于船舶纵摇引起的重力加速度在全局坐标系下的纵向分量gL：

式中，g为9.8 m/s2。

（5）根据公式计算舱室形心处的加速度，

（6）由于采用的是线性水动力载荷预报软件，舱室形心处加速度还需加上重力加速度及分量，

（7）载荷计算点（即舱室单元的形心）到计算参考点的距离就是计算点到液货在考虑船体姿态和运动情况下的液面的距离，因此要确定当前工况下舱内液货的液面方程。假设液面方程为

因为（aL,aT,aV）为液面的法向量，所以可以计算参数a和b，

（8）液面方程中c的大小需要根据液货舱不同的装载状态采用不同的方法分别计算，详细见3.2节中液舱满载情况下确定液面方程和3.3节中液舱部分装载情况下确定液面方程的叙述。

（9）计算第i个载荷计算点到液面的距离hi，并根据以下公式计算出该处的压力：

式中：p0为设定的蒸汽压力（kN/m2）；ρ为液体密度，装液化天然气时为0.5 t/m3，装压载水时为1.025 t/m3；hi为单元形心与液面线的距离（m）。

3.2 液舱满载情况下确定液面方程

对于液舱满载的情况，如图8 所示，只需要遍历舱室周界的所有单元，将所有单元的形心坐标（x，y，z）代入到3.1节计算详细步骤中第七步确定的液面方式，计算出c的值，取最大的c值即可确定液面方程。

图8 液舱满载情况下液面示意图Fig.8 Diagram of liquid surface of liquid cargo tank with full loading

3.3 液舱部分装载情况下确定液面方程

对于液舱部分装载的情况，如图9 所示，在3.1 节计算详细步骤中第七步已经确定了参数a和b的前提下，根据液面下液体的体积与原始装载的体积不变的原则，采用折半查找的方法，不断搜索得到c的大小，具体过程如下：

图9 液舱部分装载情况下液面示意图Fig.9 Diagram of liquid surface of liquid cargo tank with part loading

（1）根据设定的舱室装载的液货质量和液货密度计算出液货原始体积；

（2）将舱室周界所有单元的形心代入到液面方程计算出最大的cmax和最小的cmin；

（3）在cmax到cmin范围内采用折半查找方法确定c的大小，如图10所示，即每次取

图10 折半查找方法确定液面方程Fig.10 Determination of liquid level equation with binary search algorithm

构建液面方程，判断液面下方舱室的体积是否与原始装载的液货体积相等，计算方法如下所述，如果不相等就用c的值代替cmax或者cmin后重新计算c的值，然后再次计算液面下的方体积后进行比较，直到液面下的舱室体积与液货原始体积相同。

液面下的舱室体积计算方法如下：

（1）使用液舱周界的单元和节点信息构建开源的四面体网格划分方法TetGen 所要求的输入文件[13]，然后将液舱内部空间全面划分成四面体网格，如图11所示。

图11 四面体网格划分示意图Fig.11 Diagram of tetrahedral mesh generation

（2）计算每个四面体与液面的相交点，如果不存在相交点，则根据液面的法向量和四面体的形心判断四面体是否在液面下方；如果存在相交点，则根据相交的情况，如图12 所示，将分割出的五面体再次划分成两个四面体，然后判断每个四面体是否在液面下方。

图12 液面与四面体相交情况示意图Fig.12 Diagram of intersection of liquid surface and tetrahedral mesh

（3）计算所有液面下方四面体的体积就可以得到液面下的舱室体积。

4 波面修正

根据CCS 规范的要求，当采用线性水动力载荷预报软件进行波浪载荷预报时，需要进行波面修正。波面修正的主要要求就是根据波浪实际的形状对船体外表面的波浪载荷进行修正，如图13所示。

图13 波面修正示意图Fig13. Diagram of wave surface correction

设计波某一时刻在局部坐标系中的方程为

式中，A为设计波的波幅，ω为设计波所处工况预报时对应的波浪圆频率，t为某一时刻，k使用以下公式计算：

局部坐标系说明如下：

（1）局部坐标系原点为（0，0，d），其中d为设计波所处工况的吃水。

（2）局部坐标系Z轴与全局坐标系Z轴一致。

（3）局部坐标系X轴和Y轴是将全局坐标系绕Z轴顺时针旋转设计波所处工况对应的浪向角后的全局坐标系X轴和Y轴的位置。

在确定了波浪表面方程的前提下针对船体外壳的单元进行波浪载荷修正，详细的修正流程如下：

（1）根据设计波相关参数构造出在局部坐标系下的波浪方程。

（2）将每个单元的形心坐标（x，y，z）转换到局部坐标系下（x,，y,，z,）。

（3）将x,和设计波的相位角，即ωt等于设计波的相位角，代入波面的方程，计算出Zbl并与z,比较，根据不同的情况进行波浪载荷的修正，如图14所示，详细描述如下：

图14 波面修正规则示意图Fig.14 Diagram of wave surface correction rule

（i）对于位于当前工况水线以上的单元

如果z,＞Zbl，则单元压力为0；

如果z,＜Zbl，且z,＞0，则计算压力ρg(Zbl-z,)；

如果z,＜0，则不进行修正。

（ii）对于位于当前工况水线以下的单元

如果z,＞Zbl，则单元压力为0；

如果z,＜Zbl，如果静压+动压指向船体外部，则压力为零；

如果静压+动压指向船体内部，保持原压力。

某工况的设计波形状和修正后的波浪载荷如图15所示。

图15 波面修正结果图Fig.15 Results of wave load with wave surface correction

5 载荷调平

进行波面修正后，可能会导致船体的总力系不平衡，其中的垂向不平衡力、绕X轴及Y轴的不平衡力矩，可以通过浮态调整的方法进行平衡；其中的横向和纵向不平衡力及绕Z轴的力矩可以通过调整加速度的方法进行平衡。整个调平的顺序是先调整船体垂向、横向和纵向的力平衡，再调整绕Y轴和X轴的力矩平衡，在调整力矩平衡的过程中可能会导致力的不平衡，因此需要在调整完力矩后再次判断力是否平衡，然后循环以上过程，直到力和力矩都满足允许的差距范围，最后再调整绕Z的力矩平衡，具体的流程图如图16所示。在调整的过程中有时会出现无法都达到允许的差距范围而一直循环调整的情况，程序会设定最大循环次数来强制结束，功能界面如图17所示。

图16 载荷调平流程Fig.16 Loading balance adjustment process

图17 载荷调平功能界面Fig.17 GUI of loading balance adjustment setting

5.1 重力与浮力调平

重力和浮力调平的目标是使重力和浮力的大小相等，方向相反。重力浮力调平步骤如下：

（1）计算波面修正后的重力G和浮力F。

（2）如果 |G|＞|F|且超出允许的差距范围，抬高水线高度，大小为

式中，S为水线处船体剖面面积。

（3）对于位于载荷预报时设定的水线以下的单元，静压修正为P0+ρgΔh，其中，P0为单元原来的压力，ρ为海水密度。

（4）如果 |G|＜|F|且超出允许的差距范围，降低水线高度，大小为

（5）对于位于载荷预报时设定的水线以下的单元，静压修正为P0-ρgΔh。

（6）以下面公式重新计算d后再次进行波浪载荷修正

式中，d为第一次使用时取设计波所处工况的吃水。

（7）重复第（1）步到第（6）步，直至重力和浮力在允许的差距范围内或者达到最大迭代次数。

5.2 沿X轴和Y轴方向力调平

先计算船体所有载荷沿着X轴方向和Y轴方向的不平衡力Fx和Fy，然后用-Fx和-Fy分别除以整船重量得到X轴方向和Y轴方向的平动加速度的调整量，最后修改设计波所处工况惯性载荷的X方向和Y方向的分量。

5.3 绕Y轴弯矩调平

绕Y轴弯矩调平的目标是通过绕局部坐标系的Y轴旋转水线面，如图18所示，使重心和浮心在同一个平行于YZ平面的平面上。其中：局部坐标系原点为（船舯所在的x坐标值，船舯所在的y坐标值，调整后的吃水d），局部坐标系三个轴的方向与全局坐标系三个轴的方向相同。弯矩调平的步骤如下：

图18 绕Y轴弯矩调平示意图Fig.18 Diagram of bending moment around y-axis balance adjustment

（1）计算重力产生的绕Y轴弯矩MGy和浮力产生的绕Y轴弯矩MFy。

（2）如果|MFy|＞|MGy|并超出允许的差距范围，则浮心的x坐标在重心的前面（即：xF＞xG），需要船艏上抬，即水线面绕着局部坐标系Y轴顺时针转Δψ弧度，反之则波面逆时针转。

式中，ψ为设计波所处工况的纵摇角。

（3）根据新的水线面方程计算每个单元的静水压力（即：根据单元距离水线面高度的变化，对静水压力进行相应的增加或者减小）。

（4）重新进行波浪载荷修正，然后再次计算修正后的重力和浮力分别产生的绕Y轴弯矩MGy和MFy。

（5）如果MGy和MFy的大小关系没有发生变化且超出允许的差距范围，则水线面继续沿着上次转动的方向旋转Δψ。如果浮力力矩与重力力矩的大小关系发生变化且超出允许的差距范围，则水线面逆着上次转动的方向旋转Δψ弧度，计算公式为

（6）重复第（3）步到第（5）步，直至绕Y轴弯矩在允许的范围内相等或者达到最大迭代次数。

5.4 绕X轴弯矩调平

绕X轴弯矩调平的目标是通过绕局部坐标系的X轴旋转水线面，如图19 所示，使重心和浮心在与Z轴平行的一条线上。其中：局部坐标系原点为（船舯所在的x坐标值，船舯所在的y坐标值，调整后的吃水d），局部坐标系三个轴的方向与全局坐标系三个轴的方向相同。弯矩调平的步骤如下：

图19 绕X轴弯矩调平示意图Fig.19 Diagram of bending moment around x-axis balance adjustment

（1）计算重力产生的绕X轴弯矩MGx和浮力产生的绕X轴弯矩MFx。

（2）如果|MFx|与|MGx|大小不同并超出允许范围，则根据情况绕着局部坐标系X轴旋转水线面Δφ弧度，

式中，φ为设计波所处工况的横摇角。

（3）根据新水线面方程计算每个单元的静水压力（即：根据单元距离水线面高度的变化，对静水压力进行相应的增加或者减小）。

（4）重新进行波浪载荷修正，然后再次计算修正后的重力和浮力分别产生的绕X轴弯矩。

（5）如果浮力力矩与重力力矩的大小关系没有发生变化且超出允许的差距范围，则水线面继续沿着上次转动的方向旋转Δφ。如果浮力力矩与重力力矩的大小关系发生变化且超出允许的差距范围，则水线面逆着上次转动的方向旋转Δφ弧度，计算公式为

（6）重复第（3）步到第（5）步，直至绕X轴弯矩在允许的范围内相等或者达到最大迭代次数。

5.5 绕Z轴的弯矩调平

绕Z轴弯矩调平的目标是在以上载荷调平的基础上，通过修正设计波所处工况的惯性载荷来使重力和浮力产生的绕Z轴的弯矩平衡。弯矩调平的步骤如下：

（1）计算重力产生的绕Z轴弯矩MGz和浮力产生的绕Z轴弯矩MFz；

（2）计算船体绕Z轴的转动惯量Izz；

（3）计算绕Z轴的转动加速度修正量，

（4）计算沿着Y轴方向的平动加速度的修正量，

（5）将w,z和a,

y增加到设计波所处工况的惯性载荷的相应分量上。

6 实例计算与结果分析

本文以14.7 万方薄膜型LNG 整船模型为例来验证本系统的有效性。该船有4 个货舱，有限元模型共有439 763个单元、218 851个节点，如图20所示。

图20 整船模型图Fig.20 Overall membrane LNG carrier model

6.1 模型参数与边界条件

设置的计算参数如图21所示，设置的边界条件为：（1）尾端节点1，约束横向线位移；（2）首端节点2，约束三个方向线位移；（3）尾封板左右两侧节点3 和4，约束垂向线位移。

图21 计算参数和边界条件示意图Fig.21 Calculation parameters and LBC setting

6.2 计算工况

根据CCS“钢规”[5]的要求，薄膜型LNG 整船的计算工况由装载工况和波浪工况组合而成。装载工况从装载手册中选取，一般包括满载、正常压载、一舱为空和相邻两舱为空四种工况，波浪工况考虑船中弯矩、L/4处剪力、3L/4处剪力、船首垂向加速度、船首横向加速度和横摇角分别达到极值时的波浪工况，一般共计48 个计算工况。由于篇幅所限，本文仅给出了满载工况下各波浪工况的设计波控制参数，如表1所示。针对每个计算工况按照计算流程进行了舱内载荷计算、波面修正和载荷调平，最后计算出了应力响应。