基于动力效应的船体远场冰载荷测量与识别

何帅康，陈晓东，孔帅,2，季顺迎*

1 大连理工大学 工业装备结构分析国家重点实验室，辽宁 大连 116024 2 中国船舶科学研究中心 水动力学重点实验室，江苏 无锡 214082

0 引 言

全球变暖导致的北冰洋冰雪融化与北极航道的逐步开发，使商船航行于北极地区成为可能。2013 年中远海运集团所属“永盛”轮在北极东北航道的首航，拉开了我国商船北极航行的序幕[1]。相比于传统航道，北极航道大幅缩短了东亚与北欧地区的海上运输距离，其商业价值不容忽视。海冰是船舶在冰区航行中所面临的主要威胁之一，船−冰相互作用可能导致船体局部破坏甚至整船倾覆，因此准确、有效地衡量海冰对船舶结构的影响是保障航行安全的重要前提。

国内外针对冰载荷的研究主要采用理论分析、数值模拟[2]、模型试验[3]和实船走航测试等方法[4]。其中，实船走航测试能得到最为真实可靠的船体冰载荷特性[5]。早期冰载荷的确定采用布置压力盒的方式对冰载荷进行直接测量，而在冰载荷的作用下，压力盒很容易脱落或损坏，因此目前多采用测量船体冰激响应（如船体应变、加速度）的方法，间接识别冰载荷。Leira 等[6]和Suyuthi 等[7]针对挪威海岸警卫队KV Svalbard 号破冰船在巴伦支海的航行，安装光纤光栅传感器测量了船−冰作用下的肋骨剪切应变；Suominen 等[8-9]和Kotilainen 等[10]在波罗的海对S. A. AgulhasⅡ号极地科考补给船进行了全船冰载荷、冰激振动、噪声等的测量；Frederking 等[11]在纽芬兰东北海岸对加拿大海岸警卫队Terry Fox 号破冰船进行小冰山撞击试验，分析了该船各海冰参数对总冰力的影响；Lee 等[12]和Kwon等[13]在楚科奇海对韩国Araon号破冰科考船进行了多次局部冰压现场试验和破冰性能研究；Uto 等[14]针对日本Teshio 号破冰巡逻船在鄂霍次克海南部冰区的航行，进行了破冰模式与冰载荷关系的研究试验。国内的现场测量主要针对“雪龙”号破冰船的极地科考航行，进行了冰激振动加速度、应变信息的测量及分析[15]。

根据船体冰激响应反求冰载荷是冰载荷测量研究中的难点[16]。Leira 等[6]根据肋骨剪应变，通过应变与应力之间的关系求出剪应力，并将剪应力沿肋骨积分求得冰载荷；Suominen等[8]考虑到冰载荷在相邻肋骨间的传递，提出基于静力平衡的影响系数矩阵法，并进行了理论推导和数值分析；Kwon 等[13]针对IBRV Araon 船的肋骨剪应变，分别采用直接计算和影响系数矩阵法进行局部冰压力的识别，验证了影响系数矩阵法的工程可用性。然而，无论是直接计算还是采用影响系数矩阵法，都仅在静力平衡的条件下[17]进行识别计算，而冰载荷作为典型的随机动载荷，其动力响应不应被忽略。

本文考虑到冰载荷的动力响应，拟采用格林函数法[18-19]建立冰载荷识别模型，对“天恩”号极地运输船于2019 年在北极东北航道航行的冰载荷进行反演识别，并采用正则化方法避免病态问题求解过程中的不适定问题，识别出局部冰载荷的时间历程，以此为冰区航行船舶的船体冰载荷监测提供一定参考。

1 基于“天恩”号冰级船的船体结构响应监测

“天恩”号是中远海运特种运输股份有限公司针对极地气候打造的首艘冰级轮，于2017 年完成建造，2018 年取道北极东北航道顺利完成北欧首航，其船体主要参数见表1。“天恩”号的冰级为 B1 级（中国船级社入级规范），可自由通行0.8 m厚的当年冰冰区。2019 年7 月20 日，“天恩”号由江苏太仓港出发，取道“冰上丝绸之路”，途经长江水域、中国沿海、日本海、白令海、北冰洋、北海、丹麦海峡、波罗的海于8 月17 日到达目的地挪威。其中，在8 月1 日～8 月8 日期间，“天恩”号在俄罗斯破冰船的引航下顺利穿过了海冰覆盖区域，图1 为“天恩”号在破冰船引航下于冰区航行的航线与实况图。

表1 “天恩”号极地运输船船体主要参数Table 1 Main dimensions of polar class carrier vessel Tian'en

图1 “天恩”号极地运输船的航行路线及沿途冰况Fig. 1 The arctic route of Tian'en and ice conditions

船舶在冰区航行时，海冰对船体造成的荷载远大于波浪、风等环境因素，因此在分析时认为船体的结构变形响应主要由海冰与船体碰撞所产生。船舶舷侧肋骨作为主要的承力构件，在船体监测中被重点关注，因此一般将肋骨剪切应变作为船体结构响应监测信号和载荷识别信号。

由于“天恩”号的船−冰作用区所在位置处为抗横倾水舱，无法安装应变信息采集装置，因此将监测区域设置于远场的水手长储物间左舷侧外板区域(肋位号231～242)，监测区域和船−冰作用区大体位置如图2(a)所示。“天恩”号船艏区域为进行冰区加强采用混合骨架式建造，图2(b)所示为监测区域内部实景，椭圆框内为测点布置区域，区域内加密的肋骨(每半个肋位处设置一根肋骨) 均为T 型钢材料，应变传感器均安装在T 型材腹板处。

图2 “天恩”号极地运输船的监测区域及测点信息Fig. 2 The monitoring area and gauge map of Tian'en

为了扩大应变监测范围，选择在每个肋位处布置测点，在11 根肋骨上分结构上、下位置布置了24 个测点，其中测点1～22 号的布置如图2(c)所示。剩余2 个测点与测点17 号，18 号在同一肋位上，但安装位置更靠近于肋骨翼板，如图2(d)所示，用于比较靠近外板与远离外板的应变信号强弱。考虑到电阻应变片传感器具有灵敏度高、尺寸小、粘结牢固、适用条件好等特点，在每根肋骨的上、下2 个测点处均呈 90°直角布置且沿船宽方向对称分布了2 个电阻应变片，应变片的安装方式以及桥路连接方式分别如图2(e)与图2(f)所示。

根据应变状态理论，在yz平面(肋骨腹板所在平面)内任意方向的线应变[20]可表示为

由以上两式可得

因此，肋骨平面内任意一点的剪切应变可由该测点处所布置的两片垂直对称分布的应变片测量得出。从而，肋骨上、下2 个测点的剪切应变差值可表示为

式中： ∆γ为2 个测点间的剪切应变差值，可用于后续的载荷识别计算； γup为肋骨结构上测点处的剪切应变； γlow为肋骨结构下测点处的剪切应变。

将两电阻应变片以图2(f)所示形式接入电路中。其中：R1和R2分 别代表 45°和 −45°布置的应变片；Rc为 固定阻值电阻；AC端接电源E；UBD为输出电压即所测量的电压。应变变化导致R1和R2阻值变化，UBD与 ε(45)−ε(−45)， 即 γyz，呈线性关系，因此通过测量电路两端电压值可相应转化为肋骨处的剪切应变。

2 船体冰载荷的远场识别方法

2.1 基本原理

根据线性时不变系统的叠加原理可知，任意动态载荷引起的应变响应可用一系列单位脉冲载荷的响应叠加而成[21-22]，如图3 所示。因此，作用在船体结构上的冰载荷和冰激应变响应之间可以表示为卷积积分的形式：

图3 脉冲函数表示的载荷时程曲线Fig. 3 The load history represented by the impulse function

在时域内，设置 ∆t为采样时间间隔，m为采样点数，可将式(6) 的结构动力响应卷积关系式离散化[23]为

式中： εi为t=i·∆t时刻下实测的结构响应值；gi为t=i·∆t时刻下的格林函数值；fi为t=i·∆t时刻下的待识别载荷值。

冰与船的撞击通常表现出多源载荷作用的形式，而其中单个载荷引起的响应依然可以用式(7)表示，结构总的响应则是各单源载荷引起响应的线性叠加，因此多源载荷下响应与激励的关系可表示为

远场冰载荷识别中的格林函数矩阵G的条件数cond(G)很大，当输入数据有微小扰动或计算过程中的舍入误差都可能引起输出数据的很大扰动，导致解严重失真[24]。而在数据测量中，由于测量的系统误差和信息采集中不可避免的噪声影响，船体实测结构响应和实际响应存在一定的差异，可能导致方程无解或者有无穷多解，因此需引入正则化方法来克服冰载荷识别问题的病态性。

正则化的基本思想是给不适定问题添加一个附加条件，求得一组可以接受的解，再从中选出最适合方程的解。本文采用共轭梯度最小二乘(CGLS)迭代算法对冰载荷识别反问题进行求解，其基本思想是应用共轭梯度法求解下式的最小值问题[25]：

式中，Fm为m次迭代后的待识别冰载荷的解。

具体的冰载荷识别流程图如图4 所示。其主要思路是通过有限元分析，分别在每个子加载区域加载单位冲击载荷，求出各加载区域与每个测点响应之间的格林函数矩阵；通过简单处理实测的应变信息找到船体冰激响应数据，构建方程(10)并通过迭代正则化方法解出方程的近似稳定解，即为所识别的冰载荷时间历程。

图4 冰载荷识别的流程图Fig. 4 The flow chart of ice load identification

2.2 数值算例验证

为评估基于格林函数和正则化方法的时域反卷积识别方法对“天恩”号部分舷侧结构所受冰载荷的识别能力，建立相应的有限元模型并在模型上进行多源远场载荷识别的数值算例验证。根据“天恩”号的船艏型线图、船艏基本结构图和现场实拍照片，采用ANSYS 软件建立船艏区域部分舷侧外板的有限元模型。有限元模型完全采用SHELL 181 单元进行建立，模型的主要参数和示意图分别如表2 和图5 所示。骨材尺寸信息及板厚通过现场测量得到。模型四周边界条件设置为刚性固定。

图5 “天恩”号测试区域的有限元模型Fig. 5 The FE model of monitored area of Tian'en

数值算例的测点及加载区域示意图如图6(a)所示。由于强肋骨结构与普通肋骨尺寸、厚度信息不同，4 组测点从右至左分别设置在肋位号239，237，236 和235 处的4 根普通肋骨上。图中三角形处为肋骨的测点选取位置，同一肋骨的上、下两个测点组成一个测点组，从右至左命名为第1～4 个测点组，用于定义肋骨处的剪应变差值。矩形框内为设置的子加载区域，子加载区域水平宽度为船艏区域的肋骨间距，高度设置为水线上、下各0.5 m。首先在图中的1 号和2 号子加载区域上分别施加1 Pa 的单位冲击载荷，载荷时间步长设置为0.005 s，得到各测点的响应，其中1 号子加载区域到各测点的格林函数响应如图6(b)所示。将各子加载区域对应各测点的响应代入式(8)求得G中的每个元素，组装后求出监测区域对应加载区域的格林函数矩阵G。

表2 “天恩”号极地运输船的有限元模型主要参数Table 2 Key parameters of the FE model of Tian'en

图6 数值计算的测点、加载区域示意图以及单位冲击载荷下的格林函数响应Fig. 6 The sketches of loading area, gauging matrix and the response of Green's function under specific surface loading in the numerical simulation

求出对应的格林函数矩阵后，在子加载区域分别施加已知的载荷时程，并将测点组的剪切应变响应代入式(10)解出载荷时程。以往的冰载荷实际测量表明，冰载荷呈现为类似三角波的形式[10]，再考虑到船−冰作用的随机性，分别在1 号和2 号子加载区域中施加单三角波和多三角波形式的载荷。图7 所示为无噪声影响下的多源远场载荷识别算例的结果，图中实线表示施加的均布载荷时程，识别载荷的采样频率设置为50 Hz，从图中可以看出，基于格林函数和正则化方法的冰载荷识别模型能较好地反演出真实载荷时程。

图7 数值算例的计算结果Fig. 7 The results of the numerical simulation

3 冰区航行船体冰载荷的识别计算

为验证格林函数和正则化方法对实船冰激响应的识别能力，基于所建立的“天恩”号部分外板有限元模型，对2 个典型的船−冰相互作用工况下的船体冰激应变响应进行冰载荷识别计算，并对不同工况下的船体冰载荷特征进行对比分析。

3.1 计算工况的选取

根据冰区航行过程中数字摄像机所采集的图像信息，本次冰区航行中船−冰相互作用下的海冰破坏模式主要包括挤压破碎、劈裂破坏和翻转滑移。其中挤压破碎主要出现在海冰碰撞点处的局部区域，整体的劈裂破坏或者翻转滑移取决于海冰尺寸及船−冰碰撞角度。以下选取2 种不同海冰破坏模式下的冰激应变响应进行冰载荷识别分析，2 种计算工况的主要参数见表3。

图8 所示为2 种不同工况下“天恩”号极地运输船与浮冰发生相互作用的视频图像。其中工况1 为“天恩”号航行过程中撞上的最大浮冰块，船体左舷与一块小面积的平整冰迎面撞上。其作用过程可以描述为：最初的浮冰与船体碰撞，产生裂纹；裂纹开始扩展且碰撞点处的海冰发生局部挤压现象，破碎成大量小块碎冰；同时，浮冰的局部裂纹扩展至浮冰整体，使之发生劈裂破坏，从而大块浮冰被分割成多块较小的浮冰；最后，船−冰不再碰撞，冰载荷事件结束。工况2 为一块较大浮冰块与艏部左舷区域相撞，局部的挤压破碎现象不明显，碰撞产生的能量主要转化为浮冰块的动能，使之发生翻转滑移，与船体不再发生碰撞。

表3 计算工况主要参数Table 3 The main parameters of navigation conditions

图8 “天恩”号船体与浮冰的相互作用Fig. 8 The interaction between Tian'en and the ice floe

根据船−冰作用图像信息所记录的时间可从应变采集系统中找到对应的冰激应变时程。在识别计算中选取了肋位号235，236，237，239，240 和241处的6 根相同规格的普通肋骨，分析肋骨上、下两个测点之间的剪切应变差值，其时程曲线如图9所示。其中：图9(a)为工况1 下的冰激剪切应变差值时程曲线，应变变化周期约为1.8 s，应变峰值约为9 µm/m；图9(b)为工况2 下的冰激剪切应变差值时程曲线，应变变化周期约为0.3 s，应变峰值约为1.2 µm/m。从中可以明显发现，冰激应变周期和峰值与海冰破坏模式有很明显的关系，局部挤压状态下的海冰，碰撞后发生小块碎冰堆积现象，与船体持续碰撞，产生的应变响应较高且持续周期较长。而发生翻转滑移现象的小块浮冰，碰撞后会“弹开”，冰载荷事件下仅表现为一次短暂的撞击现象，因此产生的应变响应较低且周期较短。

图9 不同工况下的肋骨剪切应变差值测量结果Fig. 9 The measured results of frame shear strain difference under various navigation conditions

以下选取工况1 和工况2 进行局部冰载荷识别计算。2 个工况的测点及子加载区域的选取如图10 所示，三角形处为肋骨的传感器布置位置，上、下两个应变片传感器组成一个测点组，矩形框内为子加载区域，共设置了6 个测点组和6 个子加载区域。6 个测点组响应分别是肋位号235，236，237，239，240 和241 处的6 根普通肋骨上、下两个应变传感器的剪切应变差值；6 个子加载区域根据视频信息中船−冰作用点的位置确定，子加载区域水平宽度为船艏区域的肋骨间距，高度设置为水线上、下各0.5 m。

图10 不同工况的测点及子加载区域设置Fig. 10 The monitored area and gauge map under different conditions

3.2 船体冰载荷的识别结果

在确定了测点及子加载区域后，分别在6 个子加载区域上施加1 Pa 的单位冲击载荷，载荷时间步长设置为0.005 s。根据各子加载区域对应的测点响应信息求得对应的格林函数矩阵，将矩阵G和实测肋骨剪切应变差值代入式(10)计算，通过CGLS 正则化算法解出冰载荷时间历程。工况1、工况2 计算出的冰载荷时程曲线如图11 所示，2 个工况均可定义为一次冰载荷事件，没有明显的二次碰撞现象，其船−冰碰撞区域面积、作用周期及冰载荷峰值均有所不同。对于每个工况，各子加载区域冰载荷的周期相同而幅值不同，均可以等效为三角波的形式。从曲线光滑程度上看，工况1 的冰载荷时程曲线较为平滑，信噪比较高，而工况2 下的局部冰载荷数值受噪声影响较大，导致冰载荷识别在一定程度上失真。此外，将计算出的各工况下子加载区域最大冰载荷和测点的最大应变响应绘制在一幅图中，如图12 所示。从中可以明显看出，船−冰相互作用下的冰载荷的作用周期与船体的应变响应周期一致，对船体局部区域应变的监测可以有效地实现对冰载荷作用周期的监测。

图11 不同工况下的冰载荷识别结果Fig. 11 The identified ice load under different conditions

图12 不同工况下的最大冰载荷时程曲线与测点的最大应变响应时程曲线Fig. 12 The time-history curves of maximum ice load and maximum strain response of measuring point under different conditions

3.3 船体冰载荷特征分析

船−冰相互作用是一个动态的加载过程，根据计算出的各子区域载荷时程，可以看出冰载荷的有效作用面积往往不局限于一个肋位。工况1 的有效作用肋位为3 个，工况2 的有效作用肋位为2 个，且冰载荷随着肋位的变化表现出衰减的趋势。其原因是随着船−冰相互作用过程的进行，浮冰逐渐破碎，导致有效作用在各肋位上的浮冰尺寸总体呈现出减小的趋势。

识别计算出的冰载荷为均布载荷的形式，实际上是冰压力，而船−冰作用下的冰层压力分布并不是均匀的，水线附近会产生一个明显的高压力区。因此在船体结构监测中，通常采用线载荷(kN/m)的形式描述冰载荷[6]。线载荷和冰压力的关系可由式(12)表示：

式中：FL为线载荷；Fice为识别计算中求得的冰压力；B为划分子区域的宽度。

通过计算转换，可将2 个工况下的总冰载荷以线载荷的形式表示，如图13 所示。工况1(冰厚1.2 m)对应的最大总冰载荷为323.5 kN/m，工况2(冰厚1.5 m)对应的最大总冰载荷为69.3 kN/m。可见在不同海冰破坏模式下，冰厚已不再是影响船体冰载荷特征的最主要因素，浮冰尺寸和船−冰作用面积成为影响冰载荷大小、周期的关键因素。结合海冰图像信息分析，工况1 中有局部海冰挤压破碎堆积现象的发生，在这个过程中海冰与船体保持接触，持续加载，从而导致冰载荷在大小和周期上都大于工况2。此外，通过对比发现，不同工况下冰载荷时程曲线在持续增长阶段的斜率较为相近，即在加载过程中，冰载荷的增长速率不受冰厚、海冰破坏模式等海冰参数的影响，只与海冰破坏强度有关。

图13 不同工况下的总冰载荷时程曲线Fig. 13 The time history curve of global ice load under different conditions

4 结 语

本文针对“天恩”号极地运输船在2019 年7 月北极东北航道的航行，安装应变传感器和图像采集装置测量了冰载荷作用下的肋骨剪切应变及海冰参数。基于实测的冰激响应数据，采用格林函数和正则化方法克服了病态方程求解的不适定问题，进行了冰载荷的识别计算。计算选取了2 个不同海冰破坏模式的工况，结果表明格林函数和正则化方法能有效识别出船−冰作用下的冰载荷时程，识别效果受冰激响应信号信噪比的影响较大。对识别结果对比分析发现，冰载荷时程表现出类三角波的形式，冰载荷作用周期等同于应变响应周期。对于不同的海冰破坏模式，浮冰尺寸和船−冰接触面积成为影响冰载荷大小、周期、作用面积等特性的关键因素。加载过程中的冰载荷增长速率取决于海冰破坏强度，与海冰参数无关。今后还需要进行更多实测数据的识别分析来发现、完善冰载荷作用规律，并建立海冰参数与冰载荷之间的对应关系，为冰区船舶的结构设计与安全航行提供参考依据。