独立C型罐LNG运输船鞍座结构的数值模拟方法研究

2021-03-04娄本强张洪军董国忠

江苏科技大学学报(自然科学版) 2021年6期

娄本强，张洪军，孙立，董国忠

(1.江苏科技大学船舶与海洋工程学院, 镇江 212100) (2.上海海威斯特保温工程有限公司, 上海 201208) (3.中远海运重工有限公司, 扬州 225211)

随着我国LNG需求开始增加，在储气设施集约建设思想的指导下，我国正引领“碳中和LNG”(carbon neutral LNG)的行业实践[1]，在完善沿海省份LNG 产业链布局的同时，逐步增加 LNG 储运体系的建设力度.LNG的物流可以通过专用汽车或内河沿海船舶运输到末端用户，日本、挪威等国已有使用内河和近海船舶将远洋燃气资源由港口转运、分销至沿线LNG接收终端的先例[2].借鉴国外相关经验，我国亦可以采用中小型LNG运输船为纽带，借助长江航运便利，进行支线运输，供内陆城市使用[3].独立式C型储罐LNG运输船具有装载灵活、投资回收周期短等优点，适用于我国东南沿海及长江流域内河的LNG运输，因此，对该船型的研究具有现实工程意义和经济价值.

C型储罐与船体结构之间相互独立，鞍座结构支撑和限制储罐移动起到保护的作用.在对C型罐LNG运输船进行结构强度计算时，除计入波浪载荷外，还应考虑因加速、减速、转弯和升沉等引起的惯性载荷项.《国际散装运输液化气体船舶构造与设备规则》(简称IGC规则)[4]在惯性载荷合成方法上作了重要修改，用三维加速度椭球代替二维椭圆图，虽然增加了计算难度，却更接近实际，因此，如何合理地准确施加IGC规则中的载荷，会影响支承结构的数值仿真结果，以及对安全评估合理性的判断.

文中以某22 000 m3独立双耳筒C型储罐LNG运输船作为研究对象，利用自编软件工具AWQLoading[5]，将静水载荷、航行中的等效波浪载荷及惯性载荷等加载到典型工况的全船有限元模型上，进行直接计算.文中摒弃线性迭代、转向非线性有限元计算，讨论接触非线性边界以及动/静摩擦系数设置对鞍座结构的仿真效果，比较不同模型在模拟层压木时对层压木支反力与鞍座结构应力分布的影响.在现有的文献中，鲜有文献对如何计入接触面间动/静摩擦力的介绍，甚至直接忽略摩擦影响.文中的工作是对现有C型罐LNG运输船结构强度计算的细化与补充，具有一定的参考价值.

1 LNG运输船及C型储罐支撑系统

1.1 LNG储罐

国际海事组织(IMO)对储罐进行了分类，包括整体型、内部绝热型、薄膜型、半薄膜型、独立型等.其中独立型为自持式储罐，采用全压式或半冷半压式，与船体结构非直接相连.通常将此类独立储罐划分为A型、MOSS型(属B型)和C型3种型式，如图1.

图1 常见的独立型LNG储罐及其所对应的全船有限元模型

C型储罐为卧式货舱容器，储罐形状特殊(图2)，独立放置于船体舱内，储罐材料采用5%～9%的镍钢，使用聚氨酯泡沫材料包覆，确保其在运输过程中，液货保持-163℃的载运温度.双耳筒型储罐较圆柱单筒型储罐能提供较大的仓容，可以更充分地利用内部空间；较三耳型储罐，双耳筒型储罐的技术和工艺相对成熟，因此，双耳型储罐的应用更加广泛.

图2 3种常见的C型储罐型式

1.2 目标船型的舱容与结构布置

LNG运输船是布置型船舶，舱容要求高，主尺度的选择要从储罐总布置及船舶性能等方面来考虑，结合吃水高、内河船舶航行视野宽阔等要求，选择储罐型式、储罐尺寸、压力要求、重量及鞍座结构等内容，与主尺度相互影响[6].文中目标船采用球鼻首型船艏，机舱设置在尾部.机舱前设置4个C型独立液货储罐.货舱卧式容器采用2个鞍座支承，其中1个是固定鞍座，另外1个是滑动鞍座，鞍座与容器间垫有层压木[7].2～4号储罐是常见的标准双筒圆柱形，舱容均为5 570 m3；为配合船艏线型，1号储罐通常采用的双筒锥形，储罐舱容约为5 300 m3.在满载情况下，双筒圆柱储罐的液货总质量约是3 066.5 t，双筒锥形储罐的液货总质量约为2 680.4 t.在中拱和中垂工况下，船体梁将产生较大位移变形.C型独立储罐长度和体积都较大，距离舯剖面较近的储罐的鞍座系统,将承受来自液货和储罐的惯性力，引起局部应力升高[8]，文中选取2号储罐作为目标.

1.3 基于IGC规则的附加惯性加速度

双耳筒型储罐的支撑鞍座呈“碗”状布置(图3)，当面临环境载荷时，由于惯性，储罐会产生运动滞后，船体与储罐之间发生耦合作用.在已发表的研究中有使用线性叠加载荷系数的子工况[9]或者用近似正余弦作用力分布[7]代替整体载荷.

图3 C形独立储罐

按照IGC规则和相关压力容器规范要求，设计中不仅要考虑低温环境，同时还应考虑船体运动对货罐的影响.在更新的IGC规则中，使用加速度椭球代替之前的二维椭圆.惯性加速度沿纵向、横向和垂向的最大分量Ax，Ay，Az是根据北大西洋中在超越概率水平为10-8波浪载荷下得到：

(1)

(2)

(3)

式中：x为液舱重心至船舯的纵向水平距离；z为实际水线面至液舱中心的垂向距离；B为是船舯最大宽度；HGM为横倾高度;Cb为方形系数;V为船速;K为系数，K=13×HGM/B，当K≥1.0时，取计算值，当K<1.0时取1.0.

对液舱内某参考点而言，内部压力Pgd计算为：

(4)

式中：aβ为动载荷引起的加速度；Zβ为沿β方向至液货舱的液柱高度，m；ρ为货物密度，kg/m3.

之前，计算该点所对应的最大压力是采用迭代法，间隔一定的倾斜角进行计算[10]，迭代法对双耳圆筒型C型罐而言稍显臃繁.文中采用解析法，用球坐标(φ,θ)表示aβ，见式(5)，利用矢量数量积理论，得到aβ·Zβ极值的解析解——式(6) (文献[11]中应用于薄膜型LNG液舱).解析解具有精度好、运算效率高等优点.

aβ=(0,0,1)+(Axcosφ,Aysinφcosθ,Azsinφsinθ)

(5)

aβ·Zβ=(zA-zi)+

(6)

值得注意的是纵舱壁通常设置为水密，左右两侧的压力需要分开计算.使用PTC.MathCAD中编程实现解析解的图形化结果，输入2号储罐右侧外型，将数据和计算集成在一个工作表中，得到椭球无因次加速度矢量与重力夹角β= 35.6°，并得到可视化结果：液货高度Zβ=13.74 m，对应的加速度aβ=1.50 m/s2(图4).

图4 第2号液货舱右舱的PTC.MathCAD软件的可视化结果

1.4 基于设计波法的全船有限元直接计算

采用Patran对全船进行有限元建模，使用Nastran求解.总体坐标系原点设置于船体的形心，根据X、Y、Z左手坐标系统建立模型，其中X轴沿船长方向船艏为正；Y轴沿船宽方向右舷为正；Z轴由基线开始，以沿型深方向向上为正.船体的壳板、船体强框架腹板、底部肋板和液货舱内强框架腹板等结构使用板壳单元(shell)；强框架面板、扶强材、纵骨及船底、舷侧上的加强筋和液货舱支撑构件使用梁单元(beam)模拟.采用等效设计波方法将波浪附加弯矩和静水弯矩施加在船体上的，其中波浪载荷是准静态载荷，采用水动力计算软件PRECAL®计算目标船在规则波作用下的船体运动响应、波浪载荷，得到以船舯弯矩为控制参数的等效设计波，波幅为垂向波浪弯矩最大值与典型规则波中船舯弯矩的传递函数比值.

英国劳氏船级社 (LR)的指南[12]要求对规定吃水下满载和隔舱装载等状态进行直接计算，满载状态是实际营运中的主要状态，通常为进一步确保结构强度安全，需要校核隔舱装载状态.全船有限元模型的加载过程复杂，手工操作耗时费力，自编AWQLoading辅助加载程序，通过编辑Patran输出*.bdf文件，实现调整空船质量、加载液货舱和压载水质量等功能，设置二次开发接口，将等效设计波载荷以法向力的形式加载在典型工况下目标船的船壳湿表面上，与静水浮力叠加.和船级社的商业程序(如NAUTICUS®，POSEIDON®和SHIPRIGHT®)进行对比，性能相近.

在AWQLoading中设置参数交互控件，读入合成加速度，将惯性载荷以场的形式加载，得到目标船在“准静态隔舱装载动横摇”的结果(图5)，该等效设计波主导目标参数，包括波浪垂向中拱弯矩和指定点最大压力对应的横摇加速度.

图5 准静态隔舱装载动横摇工况

2 C型储罐支承鞍座结构的模型化

2.1 C型储罐的支承系统

储罐通过支承系统与船体非直接相连，货物及储罐的重量及其动载荷作用在鞍座上.为防止热胀冷缩对货罐产生的附加作用，通常在船底设置双联鞍座形式支撑(图6)，固定鞍座位于艉部，滑动鞍座于艏部.随着C型储罐设计容积的增加，也有工程师提出采用多联鞍座支撑的方案[6].支承系统沿C型罐的周向布置，形成150°左右的包角，提供有效支撑的同时，还隔断了储罐和船体结构之间的热量交换.层压木通常选用北欧地域生长的木材，浸渍树脂，削成薄片，经高温和高压加工成型.成型后的层压木质地坚实、不易开裂[13]，还具有隔热效果好、老化速度慢等优点.

图6 双耳C型LNG船双联鞍座布置

固定鞍座的面板上设有两道挡板形成的槽，由伸入层压木凹槽内的止移键板固定层压木，在层压木与鞍座面板之间和液货舱壳体之间填充环氧树脂胶泥，起到调平和固定储罐的作用；滑动鞍座处层压木分为上下两部分，由位于鞍座面板及储罐筒体上的小挡板固定(图6).在储罐外侧沿鞍座周围布置不锈钢片，缓解剪切力.

2.2 鞍座结构的刚度系数

计算鞍座结构的整体等效刚度以两根串联的弹簧为例(图7(a))，弹簧A的刚度为K1，弹簧B的为K2，串联后弹簧组的总刚度为K.在外力F作用下，弹簧伸长量分别是x1、x2，弹簧组的总伸长量为x.假设不考虑弹簧组自重，变形均在弹性限度内，弹簧A与弹簧B受力相等，即F=K1·x1，F=K2·x2，那么K1x1=K2x2，则弹簧组的总伸长量x=x1+x2，其中x1=F/K1，x2=F/K2，所以，x=F/K1+F/K2，用公式表示为：

K=K1·K2/(K1+K2)

(7)

图7 串联弹簧与固定鞍座结构示意图

结合上述串联弹簧关系推广至多个弹簧串联，如包含不锈钢片、层压木、环氧树脂等的鞍座结构(图7(b))的等效刚度：

(8)

式中：H抗压方向的高度；A构件横截面积；E材料弹性模量;下标1，2，3分别代表不锈钢片、层压木、环氧树脂[14].

2.3 鞍座面板与层压木的面间摩擦力

当鞍座面板与层压木相互接触并发生相对运动或有运动趋势时，接触面上会产生阻碍运动的摩擦力f，按摩擦条件分为静摩擦(static friction)与动摩擦(sliding friction)两种.当只存在运动趋势而未发生相对运动时，面间存在与切向力互为反力的静摩擦力；当静摩擦力达到与法向力W和静摩擦系数μs的乘积相当程度时，会发生相对运动，摩擦力将“由静转动”，大小是法向力W与动摩擦系数μd的乘积，F为剪力，且μs>μd.

(9)

3 数值模拟方法

虽然是否考虑支撑系统内的接触和摩擦非线性对船体总纵强度的影响不大[12]，然而储罐自身及液货重量通过鞍座集中传递给船体梁，会在鞍座位置处产生较大剪力和集中载荷.在校核鞍座局部强度时，其对应力分布的影响还是值得探讨的.无论IGC规则还是各船级社并没有对层压木的模拟给出强制性要求，文中将对数值仿真方案进行介绍和讨论.

3.1 临时杆单元法

临时杆单元法是模拟层压木的传统模型方法，该方法将层压木等效为弹簧或杆单元[15].使用临时杆单元法进行线性迭代：计算并判断单元的拉/压状态，修改模型，删除受拉的层压木单元，更新模型并再次计算直至无层压木单元受拉为止.迭代次数视工况而定，通常经过2～5轮即可得到收敛结果，若计算模型太复杂将导致迭代步骤增多，人工操作容易出错.推荐借助软件工具，设置自动迭代功能，对Patran底层数据文件进行二次开发：读取*.bdf、判断拉/压状态、修改模型、调用Nastran求解器、提交任务等步骤[16].

3.2 线性间隙单元法

MSC®还提供一种可以在求解器内部完成迭代计算的单元——线性间隙单元[17].对该单元的工作原理进行如下简述：从随机选取一组向量开始，一些节点处于接触状态，保证接触部位处于受压状态，且没有侵彻；若检测到受限循环，即迭代返回到前一个状态，重启计算直至收敛.在文献[16]中详细介绍该单元的理论、控制方程及工作原理，并验证了线性间隙单元的计算与使用临时杆单元法的计算结果，发现二者等效.线性间隙单元仿真的优点在于不需要人工介入，设置简单；问题在于对某些复杂工况会产生计算不收敛的现象.

3.3 计及摩擦力的非线性间隙单元法

最理想的方案是既能保证计算的收敛性，又可以考虑层压木的“开/关”状态以及端面的动/静摩擦效果，非线性间隙单元(GAP)正好满足上述要求[16].在单元卡片CGAP中定义(图8)几何信息，初始间隙用U0来表示，载荷通过层压木传递，端点A、B的位移分别是UA、UB，刚度取决于端点的相对偏移量(即UA-UB).单元处于自由状态时(UA-UB≥U0)，分离刚度为KB(一般取KB= 10-14KA)；当偏移量小于初始间隙时，刚度恢复到接触刚度值KA[18]，在单元属性卡片PGAP中除了赋值接触时的法向刚度KA、切向刚度KT和分离时的刚度KB，还可以定义动/静摩擦系数(μs和μd).

图8 GAP单元模拟层压木

(10)

在NLPARM控制卡中对收敛条件的参数进行设置，通过“更新刚度”和“子步内迭代次数”控制计算效率，避免计算发散或歧义.提交计算后，求解器会在不删除单元的前提下，自动判断，更新刚阵，得到收敛的计算结果.

4 计算结果与分析

4.1 层压木的支反力结果

结合加载完成的C型LNG运输船的全船模型，共对5种层压木的数值仿真方案进行对比：只考虑层压木接触开/关，使用方案①“传统临时杆单元法”和方案②“非线性间隙单元法模拟层压木”；方案③“在非线性间隙单元中赋值动/静摩擦系数，μs=0.5和μd=0.25[9]”；方案④”假设无穷大静摩擦系数 (μs=∞)”；方案⑤“假设较大静摩擦系数(μs=50)的情况，忽略动摩擦阻力(μd=0)”.将满载垂荡工况下，支撑第2号储罐的滑动与固定鞍座层压木单元的支反力F结果绘制成图表，见图9，其中x轴为层压木沿周向布置距船舯的距离d，单位为m.

图9 满载动垂荡工况，鞍座支反力结果

同理，隔舱动横摇工况，2号储罐下的滑动与固定鞍座层压木的支反力结果见图10.

图10 隔舱装载动横摇工况，鞍座支反力结果

观察不同数值模拟方案下层压木支反力的异同，发现方案①和②的结果相近，可知当仅考虑接触非线性时，使用GAP仿真与传统方法的结果是一致的，验证了GAP单元的可靠性.方案③是最接近工程实际的模拟，最大支反力小于方案①和②，且分布缓和.方案④、⑤是假想模型：方案④与不计入摩擦力方案的结果相近，究其原因可能是由于无穷大的静摩擦力限制了节点的位移，使计算回归到仅考虑接触非线性的模型；当设置较大静摩擦系数(方案⑤，μs=50)时，接触面间将只存在相对运动的趋势，不存在相对滑动，可以看出周向静摩擦力将分担很大一部分的径向力，此假设模型的支反力最小.因此，对所研究的C型LNG船而言，鞍座接触面之间的摩擦对保护储罐和船体结构是有益的，静摩擦力抑制储罐移动，动摩擦力将减缓储罐的相对滑动.