基于原型测量的极地航行船舶船体冰载荷分析

刘瀛昊， 佟福山， 高良田

(哈尔滨工程大学 船舶工程学院, 哈尔滨 150001)

基于原型测量的极地航行船舶船体冰载荷分析

刘瀛昊， 佟福山， 高良田

(哈尔滨工程大学 船舶工程学院, 哈尔滨 150001)

船体与海冰碰撞引起的冰载荷一直是破冰船领域研究的重点。根据S.A. AgulhasⅡ号极地科考补给船在南极海域的原型测量数据，利用影响系数矩阵法和反演法对船体艉肩部的冰载荷进行分析。通过肋骨上不同应变传感器测得的剪应变，在MATLAB中利用系数矩阵转换，得到了不同肋骨上的冰载荷。通过建立艉肩部有限元模型，采用两种不同的离散方式对载荷区域进行离散，并利用Tikhonov正则化求解反演方程，得到外板的冰载荷分布。克服了数据处理过程中解的不适定性，使得两种方法得到的计算结果也极为相近。

冰载荷；原型测量；影响系数矩阵法；反演法；Tikhonov正则化

近年来，随着全球气候变暖，人们在极地区域的活动越来越频繁，对极地航行船舶的操纵性、安全性和舒适性要求也越来越高。船舶在极地航行中会受到来自于海浪、海冰、机械以及螺旋桨等多方面的动态激励[1]。以往研究表明，当船体与海冰发生碰撞时产生的振动等级最高[2]。因此，船体冰载荷的研究尤为重要。测量冰载荷的方法通常为模型试验和原型测量。与模型试验相比，原型测量虽然成本高昂，却是得到实际冰载荷最为精确的方法，对于冰载荷特性研究及规范的补充完善都具有重大意义[3]。但在原型测量中，船舶与海冰碰撞产生的冲击载荷并不能通过直接测量得到，只能测得船体其它部位的结构响应。

根据已知结构响应，建立系统反演模型，重构该环境下的输入激励属于结构动力学中的载荷识别问题[4-5]。动载荷识别法主要分为频域法和时域法[6]。频域法是通过系统响应函数频率与测量得到的结构响应之间的关系来确定输入激励的频谱，或利用模态坐标的变换来计算模态力的特性[7]。频域法主要适用于平稳随机载荷及动载荷的识别。时域法是根据卷积形式下载荷与系统响应的关系，利用运动方程得到载荷的时间历程，可用于非平稳性及冲击性载荷的识别，其结果直观性好，更适合于工程应用[8]。近年来，时域法在车辆交通等领域的研究应用已取得一定成果：LIU等[9-10]采用形函数法和时域内的格林函数法，快速有效地解决了车门动载荷识别问题；朱涛[11]提出了一种基于动态规划的非迭代载荷反演时域法，避免了反演过程中的误差累积，并运用于高速列车轮轨接触力的识别；ROMPANEN[12]研究了一种传感反演法，通过的系统模型的离散变换，成功求得了造纸机滚筒的线性载荷分布。

然而，受极地环境，船舶特殊性等因素影响，时域法在极地船舶领域的应用并不常见。由于人们对冰载荷的特性了解十分有限，船级社规范对于极地航行船舶设计冰载荷的计算也仅局限于简单的近似方法[13]。本文针对极地科考补给船AgulhasⅡ号在南极海域的原型测量数据，应用时域法中的影响系数矩阵法和Tikhonov正则化法对船体艉肩部肋骨处和外板处受到的冰载荷进行了识别分析，避免了反演问题的病态性，得到了极地船舶航行时船体受到的冰载荷分布变化形式及转移过程。

1 AgulhasⅡ号原型测量数据

1.1 应变仪安装位置

本次原型测试在S.A. AgulhasⅡ号破冰船上完成。该破冰船在芬兰Rauma船厂建造完成，是芬兰STX公司专门为南极海域设计的具有PC-5级破冰能力的极地科考补给船，其主尺度见表1。

表1 S.A. AgulhasⅡ号的主尺度

为了监测和评估S.A. AgulhasⅡ号在冰区航行时船体受到的冰载荷，在船艏、艏肩以及艉肩部位安装了应变传感器来测量船体形变(见图1)。本文仅研究作用在船体艉肩部的冰载荷。如图2所示，艉肩部区域共有14个应变传感器，8个位于肋骨上，6个位于板壳上，主要信息见表2。

图1 应变传感器安装位置Fig.1 The strain sensors location

安装应变传感器的肋骨分别为#39+400、#40、#40+400以及#41，除了位于#40+400肋骨上的传感器朝向船体左舷侧，其余传感器均朝向右舷侧。每个肋骨上的传感器均由两个与水平方向呈45°的应变仪构成，用来测量剪应变。用来测量正应变的传感器位于两个肋骨中间的船壳内侧。应变仪的测量范围为±4 000 μstrain，误差为±2 μstrain，相对于南极海域的航行测试来说并未超出传感器的测量范围。所有传感器的样本频率fs为200 Hz，为了消除噪音和动力振荡等短波的影响，需要使用100 Hz的低频过滤器作为测量数据的单位过滤器。艉肩部各个应变仪测得的应变数据经低频过滤后，即可用于后续载荷计算。

图2 艉肩部应变传感器位置(右手边为船艉方向)Fig.2 The strain sensors location on stern shoulder

应变仪名称位置与水平方向夹角/(°)朝向SS16.1#41肋骨+45船艏SS16.2#41肋骨-45船艏SS17.1#41肋骨+45船艏SS17.2#41肋骨-45船艏SS18.1#40+400肋骨+45船艉SS18.2#40+400肋骨-45船艉SS19.1#40+400肋骨+45船艉SS19.2#40+400肋骨-45船艉SS20.1#40肋骨+45船艏SS20.2#40肋骨-45船艏SS21.1#40肋骨+45船艏SS21.2#40肋骨-45船艏SS22.1#39肋骨+45船艏SS22.2#39肋骨-45船艏SS23.1#39肋骨+45船艏SS23.2#39肋骨-45船艏SS24.0#40+400与#41间船壳0内侧SS25.0#40与#40+400间船壳0内侧SS26.0#40与#40+400间船壳0内侧SS27.0#40与#40+400间船壳0内侧SS28.0#40与#40+400间船壳0内侧SS29.0#39+400与#40间船壳0内侧

1.2 应变实测数据

S.A. AgulhasⅡ号于2013-11-28从南非开普敦起航，在南极海域进行了为期78d的航行测试和极地科考，其中从2013-12-07/2014-02-01这段时间内可以观测到明显的海冰。当船只在层冰中转弯时，船的艉部会与海冰边缘发生碰撞，即艉肩部会受到明显的冰载荷作用；当船只开辟航道时，为了使船身前部能够破冰，船舵的操纵力要比转弯时还大，亦即作用在艉肩部的冰载荷也就更大。本文选取了两种海冰与船体右舷发生碰撞时，艉肩部承受的冰载荷进行分析，并将两次撞击分别命名为碰撞1和碰撞2，具体参数见表3，应变数据见图3和图4。

表3 船冰碰撞时的冰况数据

图3 碰撞1的应变数据 Fig. 3 Strain data of ramming 1

图4 碰撞2 的应变数据 Fig. 4 Strain data of ramming 2

2 冰载荷计算方法

2.1 肋骨处冰载荷计算

根据船体上安装的V型应变传感器测得的数据，应用影响系数矩阵法即可得到肋骨受到的载荷[14]。每个肋骨上装有两个应变传感器，用来测量肋骨上部和肋骨下部的剪应变，同一肋骨上不同应变仪测得的剪应变之差与影响系数矩阵相乘，就可以得到肋骨上的剪力(即冰载荷)。

{Fi}=[a]{Δγj) }

(1)

式中：Δγ为肋骨上部和下部剪应变之差；a为影响系数矩阵。

Δγ=γupper-γlower

(2)

[aij]=ΔγjFi

(3)

式中：Δγj(j=1…n)为肋骨j上测得的剪应变之差，n为安装应变传感器的肋骨的数目；Fi为肋骨i上产生的力。a的系数由每个时间点上肋骨产生的力F，以及这一时间点上作用在结构上的Δγ确定。通过艉肩部的肋骨矫正测量，即在肋骨上同时加载拉力并测量感应器的应变数据，就可得到影响系数矩阵a的值为

随后将其代入式(1)即可得到作用在肋骨处的冰载荷。

2.2 外板处冰载荷的计算

在结构动力学中，当输出响应(应变)和结构模型(有限元模型)已知，输入激励(冰载荷)未知时，求解该类问题的方法即为反演法。为了避免反演问题解的不适定性，可以使用单值分解法、交叉法、迭代法或Tikhonov正则化等方法。本文选取Tikhonov正则化来求解外板处冰载荷方程，流程见图5。

图5 反演法确定冰载荷流程图Fig.5 Procedure of the inverse ice-induced load determination

2.2.1 载荷离散

船体与海冰碰撞产生的冰载荷分布是不均匀，可能呈线性分布或散点状分布。多种研究表明，作用在两个肋骨间的板中央区域的压力通常要小于或等于肋骨区域的压力[15-16]。而极地规范只是简单的指出载荷作用区域为一个矩形。

鉴于测量区域应变传感器数量的限制，为了提高冰载荷计算结果的精度，需要将冰载荷进行离散化。离散后的载荷分布应具有其本身的自然特性。此外，还需增加一些前期假设以避免反演计算的病态结果。首先假设载荷是垂直加载在结构上。当冰与船体法向滑动时产生的剪应力忽略不计。其次假设载荷的方向始终是非负的，正方向为进入船体的方向。最后假设冰载荷的作用区域为0.8 m高的条状压力带，分布在两块板的中间位置(如图2)。之所以选择条状压力带，是因为设计水线(Design Waterline, DWL)位于板的上半部分，船舶在航行时水线的位置要低于设计水线；此外，两块板之间肋骨上的剪应变传感器也是对称放置的。

本文采用了两种不同的载荷离散方式。第一种离散方式如图6，固定加载区被分为8个独立的载荷区域，所有单位均为mm。这种离散方式虽然粗糙，且会降低冰载荷分布自然特性的模拟，但却能简单快速地解决Tikhonov方程最小化问题。安装应变传感器的区域在图4中分别标记为左侧(#39+400与#40之间)，中间(#40与#40+400之间)，右侧(#40+400与#41之间)。

图6 离散一Fig.6 Discretization 1

第二种离散方式(图7)共有20个独立的固定载荷区域。与离散一相比，离散二的网格更为精细和复杂，因此也更为符合冰载荷分布的自然特性。

图7 离散二Fig.7 Discretization 2

由于船壳上的正应变传感器均位于中央区域，所以离散一和离散二在中央区域的网格都更为密集。根据实际观测和摄像记录，大部分的冰载荷都作用在0.622 m的区域内，因此可将离散区域的纵向高度由0.8 m缩减到0.622 m，可有效地缩减离散区域高度，降低边界载荷对结果的影响。

2.2.2 有限元模型

根据芬兰STX船厂提供的S.A. AgulhasⅡ号图纸，利用FEMAP有限元软件进行几何建模及网格划分，由253 700个四边形和三角形单元组成，船壳、肋骨及舷侧纵桁等简单几何区域为四边形单元；肋骨和船壳连接处等复杂的几何区域使用三角形单元。图8和图9分别模型的内部和外部视图。中央深色区域为安装应变传感器的监控区域，在网格划分上更为细致，以提高计算精度。

图8 船体内侧的有限元模型 Fig.8 Finite element model from the inner side of the hull

图9 船体外侧的有限元模型 Fig.9 Finite element model from the outer side of the hull

图10中用三角符号标示出了有限元模型的边界条件。在内侧中央位置的两条纵向边界是用来连接舷侧纵桁和船壳内侧的。有限元模型的边界均设为固定边界。由于边界距离应变仪安装区域较远，因而对离散区域的结果影响很小。

图10 有限元模型的边界条件Fig.10 Boundary conditions of the finite element model

艉肩部结构材料为NV AH-50、NV DH-50及NV AH36，船壳板厚度范围为7～21 mm，材料属性见表4。

表4 艉肩部材料属性

2.2.3 系数影响矩阵

利用有限元模型的计算结果，可得到离散一和离散二的系数影响矩阵。Romppanen在2008年通过试验研究得到了固定边界离散的系数影响矩阵Z为

式中：m为传感器的数量；n为离散区域数量。两种离散方式的m和n值见表5。矩阵Z中的元素为固定值，可根据有限元模型中施加在每个压力区域的单位载荷punit确定。元素Zij根据传感器的不同类型的值由式(4)或式(5)确定

(4)

(5)

式中，γij或εij为单位载荷punit施加在压力区域j上时，传感器i产生的剪应变或正应变。即影响矩阵Z中的元素值为在区域j中，传感器i的响应与受到的载荷的比。

表5 传感器与离散区域的数量

2.2.4 Tikhonov正则化

选取Tikhonov正则化方法求解冰载荷方程。作为反演法的一种，Tikhonov正则化最早由前苏联院士 在1977年提出，可以有效地克服方程的不适定性并极为简便地得到方程的优化解[17]。冰载荷的求解方程为pice(f,t)=argminf{‖Z(f)2-e(t)‖2+λD(f)22

(6)

式中：pice(f,t)为离散的冰载荷，f为方程的解向量(2.2.1中前期假设冰载荷方向始终为正，因此式(6)中用(f)2代替f以保证解的非负性);Z为系数影响矩阵；e(t)为随时间变化的应变响应；λ为正则化参数；D为载荷区域的独立系数矩阵。正则化参数λ用来保证方程式(6)等号右侧的两部分范数在同一数量级上。在取值上，λ过小会引起解的不稳定性，过大又会导致解不精确[18]。λ初始值的公式为

(7)

将不同的初始值代入式(6)中，比较方程的优化结果，选取合适的正则化参数λ。以碰撞一为例，应用第二种离散方式测试不同的初始值，结果如图11所示，当λ=0及λ=λ0/10时，载荷结果均出现了不合理载荷波动，随着λ值的增大，意外波动消失的同时，载荷的幅值也相应下降，当λ=λ0时，其载荷幅值明显变小。而当λ≤λ0/2时，在0.63 s与0.72 s处，载荷结果十分相近，表明此时λ值的变化对合理峰值处的结果影响并不敏感。因此，综合考虑，离散二中的正则化参数取λ=λ0/2=1.652×10-19。

图11 不同正则化参数下单位载荷的对比Fig. 11 The comprasion of ice pressure with different regularization parameters

而在测试离散一的正则化参数值时发现，λ值的变化对结果并未产生明显影响，这可能是由于离散网格过于稀疏导致，因此对于第一种离散情况，取λ=0。

3 计算结果对比

3.1 肋骨处载荷结果

将碰撞1和碰撞2测量的原型数据代入式(1)，可得到作用在肋骨处的载荷随时间的变化，如图12。对于撞击1来说，#40和#40+400肋骨承受的冰载荷最大，冰载荷的主要作用时间为0.6～0.8 s；在0.63 s时，#40.5肋骨受到最大的冰载荷约为57 kN。作用在两侧肋骨(#39+400和#41)上的冰载荷均小于30 kN。对于碰撞2来说，作用于两侧肋骨的冰载荷明显大于中央肋骨上的；在1 s时，#41肋骨承受最大的冰载荷为267 kN，随后，载荷开始作用于位于中央的两个肋骨上，其最大值约为峰值的2/3；在1.7 s时，#39.5肋骨承受约200 kN的载荷，其它肋骨此时的载荷几乎可以忽略不计。于碰撞1和碰撞2而言，虽然载荷作用的时长和大小均各不相同，但是随着时间推进，冰载荷均首先作用于#41肋骨上，然后向中央两个肋骨推移，最后作用于#39+400肋骨上，即作用在肋骨上的冰载荷随着船只移动，从艏部方向开始逐渐向艉部方向转移，表明该计算结果是合理的。

图12 肋骨上的冰载荷曲线Fig.12 Ice load of the frames

3.2 外板处载荷结果

根据图5的流程，将原型测量数据代入式(6)，即可得到两种离散模式下每个离散区域的冰载荷随时间的变化，并将其代入式(8)和式(9)

(8)

(9)

式中：a=400 mm；b=155.5 mm(如图6所示)；f为冰载荷单位压力，上角标d1、d2分别为离散一和离散二，下角标i为冰载荷作用的离散区域；Qd1和Qd2即为船体与海冰碰撞期间内外板载荷的变化，如图13所示。 由于我们只关心船体承受的最大冰载荷，因此在计算过程中过滤掉了小幅波动的载荷。

从图13可以看出，离散一和离散二得到的结果十分相近，载荷的峰值均发生在同一时刻，但离散一的计算结果比离散二更大一些。在碰撞1发生时，在外板处共产生了2次载荷波动，分别为0.63 s和0.72 s时，离散一得到的最大值为70 kN，离散二计算出的峰值为48 kN。在碰撞2的过程中，离散一在1.47 s时达到峰值321 kN，此时离散二也达到最大值214 kN，随后在1.7 s时离散二出现了第二次波动，其载荷值达到108 kN，而离散一得到的其余载荷值均在50 kN以下。

除了载荷大小的变化外，尚需关注冰载荷在外板处的分布。本文以两次碰撞时最大载荷处为例，在图14中呈现了艉肩部外板的压力分布，图中坐标原点的选取与图6一致。

图13 外板承受冰载荷曲线图Fig.13 Ice load of the hull

从压力分布图可以看出，离散一与离散二得到的结果基本一致，由于离散二的网格划分比离散一更细致，得到的结果也就更为精确。

3.3 肋骨处与外板处载荷对比

表6中呈现了同一时刻船体不同部位承受的冰载荷及其分布。在碰撞1中，冰载荷主要分布在#39+400与#40肋骨及其之间的外板上，肋骨承受的总载荷为58 kN，介于离散一(70 kN)和离散二(48 kN)之间。在碰撞2中，承受冰载荷的肋骨为#39+400、#40和#40+400，总载荷达316 kN，与离散一得到的外板处载荷值(321 kN)十分相近；从载荷分布来看，离散一仅分布在#39+400与#40之间的外板上，而离散二的两个分布区域与肋骨上的载荷分布更为吻合。

图14 冰载荷峰值处外板压力分布Fig.14 Pressure distributions of peak loads on the hull

计算方法离散方法时间/s载荷分布载荷大小/kN碰撞1影响系数矩阵法0.72#39+400肋骨15#40肋骨43反演法离散一0.72#39+400与#40间下部外板70离散二0.72#39+400与#40间右下部外板48碰撞2影响系数矩阵法1.47#39+400125#4079#40+400112反演法离散一1.47#39+400与#40间上部外板321离散二1.47#39+400与#40间左上部外板214#40+400与#41间左上部外板214

4 结 论

本文基于S.A. AgulhasⅡ号破冰船于南极海域航行时的原型测量数据，应用影响系数矩阵法计算出艉肩部肋骨处的冰载荷随时间的变化，又选用两种离散方式对艉肩部外板区域进行划分，利用Tikhonov正则化反演法得到外板处的冰载荷分布。由于离散网格精度不同，离散一和离散二的计算结果存在一定偏差，两种碰撞情况下离散一的结果均偏大，表明离散区域数量的减少，降低了载荷分布的自然特性，其大小也会受到影响。整体来看，计算得到的肋骨处冰载荷与外板处的数量级相符，表明影响系数矩阵法和Tikhonov正则化反演法可以用于船体冰载荷的计算，且可相互验证，计算结果较为可信。

[ 1 ] SOAL K, BEKKER A. Whole-body vibration comfort on the S.A. Agulhas Ⅱ polar supply and research vessel during a voyage to antarctica[C]//48th UK Group Meeting on Human Responses to Vibration. Ascot: United Kingdom, 2013.

[ 2 ] SILLITOE A, UPCRAFT D, RICH K, et al. Supporting human performance in ice and cold conditions[R]. London: Lloyd’s Register, 2010.

[ 3 ] 李勇，史庆增，宋安.海冰静力作用的特点及几种典型结构的冰力模拟实验[J].海洋学报，1994，16(6): 133-141.LI Yong, SHI Qingzeng, SONG An. The static mechanical properties of sea ice and ice force experiment of several typical structures[J]. Acta Oceanologica Sinica, 1994，16(6): 133-141.

[ 4 ] HAN X, LIU J, LI W J. A computational inverse technique for reconstruction of multi-source loads in time domain[J]. Acta Mechanica Sinica, 2009, 41(4):595-602.

[ 5 ] PEZERAT C, GUYADER J L. Two inverse methods for localization of external sources exciting a beam[J]. Acta Acustica, 1995, 3(1): 1-10.

[ 6 ] UHL T. The inverse identification problem and its technical[J]. Archive of Applied Mechanics, 2007, 77(5):325-337.

[ 7 ] 智浩，文祥荣，缪龙秀,等. 动态载荷的频域识别方法[J]. 北方交通大学学报，2000, 24(4): 5-10.

ZHI Hao, WEN Xiangrong, MIAO Longxiu, et al. Dynamic loading identification in frequency domain[J]. Journal of Northern Jiaotong University, 2000, 24(4): 5-10.

[ 8 ] SUN X S, LIU J, HAN X, et al. A new improved regularization method for dynamic load identification[J]. Inverse Problems in Science and Engineering, 2013, 22(7): 1062-1076.

[ 9 ] LIU J, SUN X S, HAN X, et al. A novel computational inverse technique for load identification using the shape function method of moving least square fitting[J]. Computers and Structures, 2014, 144:127-137.

[10] LIU J, MENG X H, JIANG C, et al. Time-domain Galerkin method for dynamic load identification[J]. International Journal for Numerical Methods in Engineering, 2015, 105(8): 620-640.

[11] 朱涛. 高速列车载荷反演技术及其运用研究[D]. 成都:西南交通大学，2012.

[12] ROMPANEN A. Inverse load sensing method for line load determination of beam-like structures[D]. Finland: Tampere University of Technology, 2008.

[13] SUOMINEN M, KUJALA P. Measured ice loads and design ice loads[C]//4th International Conference on Marine Structures, 2013.

[14] BEKKER A, SUOMINEN M, PELTOKORPI O, et al. Full-scale measurements on a polar supply and research vessel during maneuver tests in an ice field in the Baltic Sea[C]//33rd International Conference on Ocean, Offshore and Arctic Engineering, USA ,2014.

[15] RISKA K, KUJALA P, VUORIO J. Ice Load and Pressure measurements onboard IB Sisu[C]//Proceedings of Seventh International Conference on Port and Ocean Engineering Under Arctic Conditions, 1983:1055-1069.

[16] MAATTANEN M, MARJAVAARA P, SAARINEN S, et al. Ice crushing tests with variable structural flexibility[J]. Cold Regions Science and Technology, 2011, 67(3):120-128.

[17] TIKHONOV A N, ARSENIN V Y. Solution of Ill-Posed problems[M]. Washington: Society for Industrial and Applied Mathematics, 1977.

[18] INOUE H, HARRIGAN J, REID S. Review of inverse analysis for indirect measurement of impact force[J]. Applied Mechanics Reviews, 2011, 54(6): 503-524.

Ice-induced load analysis for hull of an ice-going vessel based on full-scale measurement

LIU Yinghao，TONG Fushan, GAO Liangtian

(College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China)

The ice-induced load on a ship hull is a key point in the study field of icebreaker. Based on the full-scale measurement data here, the influence coefficient matrix method and the inverse method were applied to analyze the ice load on the stern shoulder area of the polar supply and research vessel S.A. AgulhasⅡ. The ice loads on different frames were therefore obtained by transforming the shear strain data measured with sensors attached on those frames and using the coefficient matrix in MATLAB. To get the ice load distribution, the finite element model in the stern shoulder area was built and two different discretization methods were used to discrete the load area. Tikhonov regularization was employed to solve the inverse equation, the ice load distribution on the outside plate was gained. The results showed that the two methods proposed can lead to similar solutions and overcome the ill-posed nature of solutions during data processing.

ice-induced load; full-scaled measurement; influence coefficient matrix method; inverse method; Tikhonov regularization

工信部高技术船舶项目(G014613002)；国家自然科学基金(51509060)；黑龙江省自然科学基金(QC2015059)

2015-06-17 修改稿收到日期： 2016-08-29

刘瀛昊 女, 博士生，1988年生

高良田 男，教授，硕士生导师，1964年生

U661.4

A

10.13465/j.cnki.jvs.2017.07.034