船冰相互作用研究方法综述

徐 莹 ， 胡志强 ， 陈 刚 ， 徐业峻

（1.上海交通大学海洋工程国家重点实验室，上海200240；2.高新船舶与深海开发装备协同创新中心，上海200240；3.纽卡斯尔大学 船舶科学与技术学院，英国 NE1 7RU；4.中国船舶及海洋工程设计研究院，上海 200011；5.中海油能源发展采油服务公司，天津300452）

0 引 言

极地资源的开发是国际热点之一。有研究表明北极区油气储量占世界新探明油气资源的25%[1]，科技的发展也使得夏季的北极航线成为可能。如图1所示，北极航线相对于目前的国际贸易航线缩短许多，能够有效提高经济效益，并压缩成本。但是，船舶及海洋结构物与冰的相互作用，严重情况下（例如碰撞）可能导致结构物破损甚至沉没，以及原油天然气泄露污染等严重问题，如图2。随着各种海洋平台和船只在北极地区的不断出现，船舶及海洋结构物与冰的相互作用问题成为研究热点。北欧、俄罗斯、美国和加拿大等环北极国家对海冰的研究较早，在船冰相互作用问题的研究中处于国际领先地位。近年来，随着北极开发的不断重视，我国也涌现出一些杰出的研究学者和有意义的科研成果。

海冰从产生到发展经历了复杂的过程，冰的生长环境、寿命和形态，以及结构物与其作用时的速率、结构形式和强度等因素都会影响到冰的力学性能。目前各国学者对冰的材料特性以及船冰相互作用问题的了解有限，因此，研究冰的本构关系，以及开发有效的数值模型是研究船冰相互作用问题的关键。

图1 北极航线与目前航线 Fig.1 The arctic routes and current routes

图2 Reduta Ordona号与冰山相撞[2]Fig.2 Reduta Ordona collided with an iceberg

最早对冰载荷的预测采用经验或半经验公式，但是经验公式法受限于适用范围。试验法是研究海冰特性以及船冰作用的最好方法之一，通过试验可以观察并研究真实冰的失效过程，获得直接的数据资料，但对试验条件要求比较严格，成本较高。随着计算机技术的不断提高，数值模拟已经成为一种经济有效的研究方法。船冰模拟的数值方法主要包括有限元法和离散元法。有限元法具有成熟的理论基础，能够处理比较复杂的力学模型、边界条件以及材料和几何非线性问题，而且商业软件较多，方便实现。但由于冰的破碎具有离散的特性，有限元法在模拟介质的间断和离散性方面能力有限，如何完善有限元方法使之能够更好地模拟冰的离散特性是一个重要课题。离散元方法是一种用途广泛的海冰模拟方法，在表现海冰的离散特性方面具有一定优势，能够较好地模拟船冰作用中海冰从完整状态到破碎状态的过程。目前离散元法开发得并不完善，如何合理地确定模型相关参数（如离散单元间的刚度）以及提高计算效率都是离散元法需要解决的重要问题。由于自然水域中冰的材料，几何和分布具有很大的随机性，因此概率法也普遍被认为适于研究冰的某些特征。本文对上述船冰相互作用的主要研究方法进行介绍，包括试验法、有限元法、离散元法和概率法，并对应用上述方法进行的船冰作用研究进行概述。

1 试验法

船冰相互作用试验研究方法包括压痕试验、压缩与拉伸试验、弯曲试验和大尺度碰撞试验等。

1.1 压痕试验

船冰碰撞时接触面上冰破坏的物理过程对冰载荷的测量至关重要，这一过程可通过压痕试验来观察。船冰碰撞时，两者一旦接触会在冰面形成若干高压区[3]，高压区的冰处于三向应力状态，随着压力提升，高压区的冰首先出现融化、重结晶等现象，低压区的冰出现微裂纹、弹性变形，之后高压区的冰被破坏，并将碎冰挤出，生成新的高压区。图3为冰的碰撞破坏示意图。

加拿大Michel等人[4]采用矩形压头对垂直冰板进行压痕试验。Devinder等人[5]对淡水冰层进行小尺度压痕试验，改变压头宽度、冰厚、压入速度和结构强度。冰——结构接触面安装接触传感器，测量压力。结果表明，低压痕速度下冰发生延性变形，高速下发生连续脆性压碎。低应变率或压入速率下，冰的变形主要是蠕变，蠕变应变通常比弹性应变大；高应变率或压入速度下，断裂前的冰变形大多是弹性的，因为蠕变需要时间。Li Chuanke等人[6]采用四种尺寸的压头在三种速率下进行小尺度压痕试验，试验发现尺度效应与位移率有关：在高速率下，尺度效应比较明显。2013年Browne[7]采用直径20mm的压头冲击多晶冰冰块，发现低温下冰易出现微裂纹并导致剥落（spalling）破坏，高温下冰易出现重结晶并以破碎（crushing）形式破坏。还有一些采用球形或半球形压头进行的多晶冰压痕试验见参考文献[8-10]。2017年，Scopper等人[11-12]采用不同形状的压头，对不同覆盖层的冰以不同压入速度进行试验，以研究这些因素对接触力以及冰破碎情况的影响。试验发现水的浸没会抑制冰的破碎，冰的表面覆盖雪或冰颗粒都会提高接触力。O’Rourke[13]对过往的压痕试验进行了回顾，分析和研究了试验中发生的冰载荷规律性振动（冰激振动）。压头附近的冰有明显的损伤层，接触压力的规律性振动主要来源于该处的损伤，中间夹杂的压力突变一般来源于大块冰的断裂和脱落，并伴随着接触面积的显著减小。

图3 冰的碰撞破坏[3]Fig.3 Ice collision damage

1.2 压缩与拉伸试验

压缩与拉伸试验一般可分为单轴、双轴和三轴加载试验。通过这些试验可得到冰的抗压与抗拉强度，复杂应力和限制条件下的冰强度测试可得到冰的完整失效面（failure envelope），失效面决定了冰在各种压力或张力状态下的失效状态。船冰碰撞中，压缩破坏是冰破坏的一种主要形式，相关试验较多。由于对试验设备要求较高，真正的三轴试验（σ1＞σ2＞σ3）和拉伸试验较少[14]，许多三轴试验为围压三轴试验（σ1＞ σ2= σ3）。

1986年，Timco和Frederking[15]进行了真实的三轴压缩试验，用不同方式对柱状冰、多晶冰和柱状冰混合物进行加载，温度为-2℃和-12℃，得到了能够描述各向异性冰的n型多项式屈服函数。Risk和Frederking[16]对多年冰进行单轴、双轴和三轴压缩试验，并将强度结果分别应用‘Tsai-Wu’和Nadreau破坏准则在p-J2平面上拟合，发现‘Tsai-Wu’准则能更好地拟合试验结果。天津大学马键钧[17]在冰室内进行三轴压缩试验，根据试验结果提出了修正的Derradji-Aout破坏准则。

应变率对冰的破坏过程和冰载荷有较大影响，许多学者在试验中就这一问题展开研究。1997年Jones[18]对淡水冰和波罗的海冰进行高应变率单轴压缩试验，应变率范围10-1-10/s，温度-11℃，发现海冰强度随着应变率提升而增大，并给出了幂次关系式。在应变率为10/s时，淡水冰强度是海冰的1.3倍。2003年，Stephen和Gagnon等人[19]对冰山冰进行单轴抗压试验，应变率范围扩展到4*10-8/s到6/s，温度-10℃。图4为应力峰值随应变率的变化曲线。随着应变率的提升，抗压强度先升高后下降，当应变率最高时，强度会回升，Stephen等人推测这与粒状冰内的裂纹扩展速度有关。

除了应变率，围压大小和温度也是关键因素。Barrette和Jordaan[20]进行三轴压缩恒定应力蠕变试验，发现围压大于某一数值时冰活化能显著增加，冰的强度降低。Sinha[21]对横向各向同性冰模型进行单轴加压/卸载蠕变试验，发现冰模型在加载之后很短时间内只表现出线弹性及滞后弹性行为。更多多轴加载试验还可参考文献[22-28]。

图4 试验测得应力峰值随应变率的变化[19]Fig.4 Peak stress plotted against strain rate from tests

1.3 弯曲试验

弯曲测试一般会导致非均匀应力区，因此弯曲强度常作为指标测试[29]，试验方法包括悬臂梁试验和简支梁试验。悬臂梁试验一般在室外进行，直接将冰梁从自然冰中切出，其中三边自由一边与自然冰相连。在冰梁的自由端施加一个渐变的载荷，当载荷达到某一值时连接端断裂，从而得到相应的弯曲强度。简支梁试验可在室内进行，一般分为三点弯曲试验和四点弯曲试验，将载荷平行且反向地作用在冰梁上[14]。

在von Bock等人[30]对模型冰进行的悬臂梁试验中，断裂面应变不均匀，上部和下部应变较大而中性轴应变较小，使得冰梁在弯曲时弹性响应持续时间较长。韩红卫和李志军[31]对淡水冰进行悬臂梁试验，测量其弯曲强度和有效弹性模量。Frederking等人[32-34]也进行了冰的悬臂梁试验。

Gagnon和Gammon[35]对冰山冰进行四点梁弯曲试验，发现冰梁的弯曲强度随应变率提高而提高，随温度升高而降低。Barrette等人[36]对冰山冰进行四点梁弯曲试验，测量梁的失效强度，建立内部结构特性（晶粒、缺陷）间的关系，发现平均晶粒尺寸越大，气泡数量越多，强度越小。Barrette认为可将缺陷分为三种级别：（1）内部空气（气泡或其他形状）；（2）平面缺陷，表现为气泡构成的二维网络；（3）形成脉络。

大连理工大学的季顺迎等人[37]在2008到2010年两个冬季，对环渤海沿岸9个测点海冰进行三点弯曲试验，发现海冰弯曲强度与卤水体积的平方根呈负指数关系，与加载速率呈线性关系，并得到了海冰弯曲强度的双参数推算公式。在此之后，季顺迎等人[38]对2011至2012年现场采集的渤海莱州湾海冰试样进行室内低温三点弯曲试验，测量海冰的断裂韧度，确定了卤水体积和加载速率共同影响下海冰断裂韧度关系式。

1.4 大尺度碰撞试验

大尺度碰撞试验以及实船测量能够更准确地体现真实冰与结构物碰撞，为此学者们进行了一些大尺度碰撞试验以及总结工作。1995年，加拿大C-CORE公司进行了真实尺寸冰山与结构物碰撞试验，通过拖拽小型冰山与垂直固定的平板进行碰撞。之后C-CORE公司[39-43]进行了CCGS Terry Fox破冰船撞击冰山试验，冰山尺寸为4-73 m，冰山质量为21-22 000 t，冰山表面温度为0℃，0.6 m深处-15.4℃。Timco[44]总结了2011年以前公开发表的大块浮冰与固定结构物的碰撞载荷，包括大尺度试验、实地试验和监测数据，以及小尺度冰碰撞试验，综合以上结果，Timco认为浮冰碰撞力与动能联系紧密且存在幂次关系，并提出了一个粗略估计浮冰碰撞力的经验公式，即，其中F为碰撞力（kN），M为质量（kg），v为浮冰速度（m/s）。Määttänen 等人[45]采用 1：3 比例对冰排与不同强度平板进行了碰撞试验。

天津大学是目前国内少数拥有海冰试验池的院校之一。该试验冰池进行了大量海洋结构物与海冰相互作用的试验研究，史庆增等人[46]在冰池中模拟了辽东湾抗冰平台与冰排碰撞，黄焱[47]对冰激柔性锥体及直立桩结构振动问题进行了冰池模型试验研究。

2 有限元法

有限元法是常用的数值方法，具有成熟的理论基础，商业软件较多，便于实现。开发广泛有效的冰模型对于数值模拟是非常重要的，但需要遵循冰的物理性质和失效方式。

2.1 本构模型与有限元

复杂的环境和形成过程使冰具有多种形态和内部结构，采用一种本构模型模拟所有冰的特点是比较困难的，目前的本构模型往往针对某种冰或在某种情形下有效。本部分将冰的本构模型和有限元模拟结合起来进行讨论。

1978年，Sinha等人[48]开发了一种针对单轴受压的粘弹性本构模型，应变包括弹性、延迟弹性和粘性应变。Karr和Choi等人[49]考虑了微裂纹对延迟弹性应变和粘性应变的影响。Ralston[50]将一种应用于岩土力学的塑性屈服函数引入柱状晶体冰。Jorddan等人[51]提出了两种方法完善应用于连续性分析的本构模型，一个是基于流变模型的方法，另一个是基于指数函数的方法，两种方法都能较好地模拟延迟弹性应变。Lade[52-53]将冰看做一种摩擦材料，通过硬化/软化塑性模型和非关联流动法则模拟不可恢复的变形。

多晶冰的自身损伤与孔隙率是两个重要变量。已有的损伤会导致微裂纹、重结晶与晶粒边界的产生，从而提高冰的延性；孔隙也会提高冰的延性。Singh和Jorddan[54]讨论了自身损伤和孔隙率的影响，用损伤函数模拟裂纹和晶粒边界的变形，并模拟了孔隙率的变化，冰的变形包括弹性、延迟弹性和粘性蠕变变形。Keller等人[55]建立了一个针对多晶冰的粘弹性本构模型，使用分布式纤维束模型构建损伤函数，将损伤变量耦合入延迟弹性变形中。

前文提到，应变率对冰的破坏过程有很大影响，本构模型的构造也应考虑这一点。Jorddan等人[56]认为，在低应变率下，蠕变与微裂纹对冰影响较大，此时冰体现粘弹特性，应变率较高时容易发生脆性失效。Derradji-Aouat[57]针对淡水横向各向同性柱状晶S-2冰，建立三维弹塑性本构模型。该模型由非裂纹子模型和裂纹子模型组成，适用于模拟冰在低应变率下的延性行为，不能模拟高应变率下的脆性失效。

试验模型冰反映了一种冰的理想状态（未变形时各向同性），其生长和冷却过程是已知的，可通过测量得到数值模型需要的重要参数。2013年，Polach[58]针对芬兰阿尔托大学冰池的模型冰——细粒冰，建立了模型冰的数值模型。发现该模型冰屈服强度相对较低，不可恢复应变发生较早，模型中未考虑粘性和温度的影响。

近几年内，采用有限元方法的船冰研究还有很多。2017年，Kolari[59]基于wing-crack模型开发了一个三维各向异性损伤模型，用以描述冰的脆性失效。将该模型嵌入ABAQUS中对冰的单轴压缩和拉伸进行数值模拟并得到了相应的失效模式和强度。Wong[60]针对层冰的弯曲断裂建立了三维本构模型，并采用有限元对圆锥形桥墩与层冰的相互作用进行了分析。Liu Zhenhui[61-62]构建了一种塑性模型模拟冰山与结构物的碰撞，并进行了船侧和船首碰撞分析。Taylor[8]采用单元移除方法，模拟压力软化和破裂过程。Gagnon等人[63]采用一种硬的可压碎泡沫建立了冰破碎的数值模型，该模型能够模拟冰压碎过程的两个重要特性：快速融化和剥离。2012年，Metrikin[64]对相对速度较小的碎冰——结构物相互作用进行模拟，结合了刚体运动与有限元方法。纽芬兰纪念大学[65]和C-CORE[66]采用有限元法对冰山与海底结构的相互作用进行了分析，以保证水下结构的安全性。

哈尔滨工程大学的任慧龙等人[67]探讨了破冰引起的船舯弯矩计算，翟帅帅等人[68]探讨了不同本构模型下的船冰相互作用。大连理工大学武文华等人[69]采用LS-DYNA分析了JZ20-2NW平台的抗冰性能，韩雷等人[70]采用ANASYS分析了冰与锥形结构物的相互作用，张大勇等人[71]对海洋结构物的冰阻力以及冰致疲劳问题进行研究。上海交通大学石础[72]建立了一个考虑温度影响的本构模型，并嵌入LS-DYNA模拟船冰碰撞，分析船体结构的损伤和变形。天津大学的卢安[73]、大连海事大学的宋艳平等人[74]也采用有限元方法开展了冰与结构物相互作用的研究。

2.2 冰的失效方式与准则

失效准则是船冰作用问题的一个重要准则，冰的失效模式包括蠕变失效，辐射裂纹，轴向断裂，屈曲，剥落和压碎，Timco建立了冰的失效图，如图5所示。

中低限制情况下（如雪），摩擦力对断层的形成影响很大，称为库伦断层；高限制情况下，摩擦滑动被抑制，易发生脆性失效，由于塑性变形作用比较明显，称之为塑性断层。冰出现库伦断层或塑性断层也与限制程度有关，Golding[75]和 Schulson[76]对这一问题进行了探讨。Renshaw等人[77]开发了一种符合试验结果的塑性断层失效应力定量模型，并采用该模型对淡水冰的压碎失效图进行完善。

近几年来，加拿大的Derradji-Aouat等人开发的针对冰的多失效面方法得到较为广泛的应用。2000年，Derradji-Aouat[78-80]针对淡水各向同性冰和冰山，根据Mroz’s的多失效面理论（multi-surface failure theory），建立了针对高应变率下的脆性失效准则，之后结合试验数据将失效模型扩展至盐水冰。应用该方法对冰山冰进行模拟和验证，能够捕捉到高应变率（大于10-3/s）下冰的弹性响应和脆性失效。Martonen等人[81]根据Derradji-Aouat的失效准则模拟了冰与结构物的相互作用问题，并与全尺度测量结果进行对比，在某些情况下吻合良好。

Schulson[82]总结了多晶冰的失效机制，尤其是脆性失效。对于冰的拉伸性能采用裂纹成核和扩展描述，对多轴加载的情况采用摩擦裂纹滑动和剪切断裂机制来描述。对于冰在压力作用下由延性到脆性的过渡，可以用裂纹末端蠕变、裂纹扩展进行解释。Kolari等人[83]引入了一种网格更新方法来预测冰的失效，由于该方法不断对网格进行更新，且不预定裂纹扩展路径，容易导致模拟时间过长。Jordaan等人[84]认为损伤力学可以较好地描述冰的软化等特性。Goldstein[85]开发了一种针对大接触区的冰失效模型。

图5 冰的失效模式与失效图Fig.5 Ice failure mode and failure map

2.3 裂纹的模拟——粘合单元法

粘合单元法（cohesive element method）是一种描述裂纹产生和扩展的方法，可在有限元中实现，近年来得到了广泛应用。粘合单元法的主要假设是粘合单元符合某种力-位移关系曲线，可称为软化曲线，该曲线的构造应符合材料的具体特性。基本参数有裂纹初始应力和断裂能。裂纹初始应力是裂纹开始发生的应力，断裂能是完全打开粘结裂纹所需要的外部能量，即曲线下方的面积。这两个参数可以通过试验获得，从而构建出软化函数。该方法模拟裂纹发展的基本过程是，有限元模型中包含普通单元和粘合单元，粘合单元按照软化曲线关系发展，直至到达最大位移时被删除，从而形成真正的裂纹。因此材料分为三个区域：（1）无损伤区域；（2）裂纹过渡区，有应力在裂纹交界面传递；（3）真正的裂纹，材料完全损坏，没有应力在裂纹交界面传递。

Gürtner等人[86]在汉堡水池进行了相关试验校准的冰模型试验，这种基于能量的断裂准则对能量释放率的确定有较高要求。Hilding等人[87]对应用于连续破碎的粘结单元法进行研究，并认为采取该方法可以对全尺度船冰作用下冰的连续破碎进行模拟。Lu Wenjun[88]采用粘合单元法模拟冰与锥形结构的作用，讨论了粘合单元法面临的挑战，例如网格依赖性，并采用两种方法减轻其影响。王峰[89]模拟了倾斜海洋结构物和层冰的相互作用。采用弹塑性软化模型来代表冰的本构关系，通过内聚力单元的失效模拟出层冰的弯曲断裂，并讨论了网格和参数敏感性问题。Gribanov等人[90]建立了一个内聚力模型来模拟柱形冰样本的断裂，通过有限元分析得到应力——应变和损伤累积的结果。

3 离散元法

有限元方法可以处理相对复杂的力学模型和边界条件，但需要明确海冰的本构模型和计算参数。此外，有限元方法对脆性破裂、压碎和冰碎片的大尺度移动模拟能力不足。离散元方法是一种用途广泛的模拟海冰特性的数值方法。离散元法可以在细观尺度下描述海冰材料的结构，能够模拟冰从完整状态到破裂状态的过程，在宏观上又能够合理地模拟海冰与结构作用时的破碎规律，获得冰荷载的冰力时程，具有一定优势。

离散元模型分为两种：硬模型和软模型。硬模型中，粒子遵循非弹性碰撞，认为粒子间的相互作用是瞬间的，遵守动量定理，需满足非渗透约束，通过冲力实现动量交换。软接触中，粒子可以轻微重叠，通过接触关系计算接触力，如弹簧-阻尼模型，分别代表动量交换和能量耗散。接触力的计算基于粒子之间的重叠量，与粒子被赋予的材料性质有关，如弹性、粘性、塑性参数和摩擦系数，典型的接触力模型如图6所示。粒子间的接触可以持续，代价是计算时间较长。离散元法主要有两个方面最为耗时：一是检测粒子单元的接触状态，在计算接触力之前，首先要判断单元间的接触状态，这一过程需要反复的接触探测，非常耗时；二是动力学方程求解。

多晶岩土材料，如冰和石头，由晶体和颗粒构成，内部存在一些缺陷。当样品尺度较小时，只有较少的缺陷，其失效方式主要是新裂纹的发生和扩展，而材料尺寸较大时，已有的关键位置的缺陷会产生较大影响，这是一种尺度效应。Selvadurai等人[91]采用二维离散元方法对船冰相互作用进行模拟，引入了尺寸相关破碎应力（size-dependent fragmentation stress）的概念，并假设应力状态为压力主导的区域为粘塑性失效，应力状态为张力主导的区域为脆性破裂。该模型体现了冰的弹性、粘塑性、脆性断裂等多种特性。

荷兰的Dorival等人[92]采用二维离散元模型模拟冰盖和刚性结构的相互作用。模型划分为六边形网络，如图7所示，轴向弹簧——阻尼器连接位于单元中心的质量点，失效准则为弹簧——阻尼元素的伸长量，断裂后相应元素从模型中移除，最终裂纹表现为沿某曲线的多个元素断裂。图8为模拟出的冰碰撞损坏后的离散元模型。为考虑冰的几何和材料随机性，假设弹簧刚度、阻尼系数、最大伸长量等参数服从高斯分布。实例分析表明该模型能够捕捉到裂纹的产生和扩展，但仍需通过试验得到较准确的相关参数。芬兰的Riikilä等人[93]采用大小随机的二维圆盘粒子，粒子之间通过无质量弹性梁连接，模拟海冰的粘弹性变形和断裂。梁的断裂采用能量法，并考虑阻尼。该模型不仅具备颗粒模型的优势，如对接触和断裂的处理比较直接，同时也能够描述冰的连续粘弹性。Tuhkuri等人[94-98]采用离散元法对层冰的弯曲、拉伸和相应的冰载荷和冰阻力进行了很多研究，并于2018年发表了关于离散元方法在船冰相互作用中应用进展的文章。

图6 弹簧、阻尼器和库伦摩擦力构成的接触模型[100]Fig.6 Contact model consisted of spring，dash-pot and Coulomb friction

图7 六边形点阵模型[92] Fig.7 Hexagonal lattice model

图8 离散元法模拟冰撞击结构物后的破碎[92]Fig.8 Ice sheet damaged after crashing against the structure

2012年，Sun Shanshan等人[99]采用三维盘单元模拟荷叶冰，研究了荷叶冰对于固定圆柱桩的载荷，考虑了冰块受到的拖曳力、附加质量力、浮力等。由于冰载荷的振动特性，参照有义波高的定义，将碰撞力和弯矩的最大三分之一的平均值定义为有义力和有义弯矩。模拟过程中认为荷叶冰尺寸一致，并忽略相邻浮冰之间的冻结。

Hansen等人[100]对锚泊定位于碎冰区的海洋结构物进行响应分析，每块浮冰用圆盘单元表示，单元——单元和单元——结构物的碰撞模拟为线粘——弹性，考虑库伦摩擦接触以防止切向相对滑动。Lau等人[101]采用基于离散元的软件DECICE进行多项分析，并采用新的接触检测序列和平行算法提高软件计算性能，缩短计算时间。Yulmetov等人[102]采用非光滑离散元方法模拟了碎冰区的冰山自由漂移和拖曳。此外还开发了一种数值冰区生成方法，能够根据现场测量的分布规律模拟出随机形状的多边形冰块。Hopskins和Løset等人[103-105]也对离散元方法进行了研究。

大连理工大学采用离散元方法对海冰进行了多年研究，并应用在船舶与海洋结构物方面。季顺迎[106]总结讨论了不同尺度下海冰的离散分布特性以及块体、圆盘、颗粒等不同形态的离散单元的应用。2013年，季顺迎等人[107]采用离散元方法模拟碎冰区海冰与船舶结构的相互作用，碎冰由三维圆盘单元构成，考虑了碎冰受到的浮力、拖曳力和附加质量。同年，对海冰与直立海洋结构作用的破碎过程进行了离散元分析[108]，将海冰单元间的粘接强度设为温度和盐度的函数，确定了不同卤水体积下海冰发生脆性破坏的单轴压缩强度，得到了冰力时程与结构冰振响应。此外，季顺迎等人[109]还模拟了冰层与锥形结构物之间的相互作用，采用球形离散单元，建立法向力、切向力与粒子位移和速度的关系，切向力极限为库伦摩擦力；通过应力修正，考虑温度和盐度对强度造成的影响。

4 概率法

冰的失效往往是多种失效现象同时发生，且由于环境和生长过程不同，导致冰在材料和几何方面存在不均匀性，因此还可以从统计学角度进行分析和预测。近年来船冰问题的概率研究也得到了很大发展。

2006年，Jorddan等人[110-111]的研究表明局部压力值具有随机性，但具备空间与时间上的关联性。CCORE的Li Chuanke和Jordaan等人[112]建立了概率模型，模拟局部压力和整体压力之间的概率平均效应，并采用向量自回归模型展示局部压力之间的空间和时间关系，通过得到的压力——面积关系和压力——方面比关系估算总的冰载荷。Kujala等人[113]认为冰载荷仅在少数情况下服从单参数指数模型。NTNU的Suyuthi[114]建立三参数指数模型描述碎冰区作用于船体的冰载荷，该模型是两种指数分布按照一定比例因子的结合体，并且更符合真实的数据分布规律。

C-CORE的Stuckey等人[115]搭建了一个基于概率算法的模型对冰山设计载荷进行评估。该模型采用贝叶斯技术确定冰山尺寸和碰撞速度的分布，并改良了面积——渗透深度模型和离心模型。Jordann等人[116]建立了一个概率断裂模型以研究冰的破坏过程。首先采用ABAQUS进行应力分析，将结果导入Matlab的PFM程序中，PFM程序能够区分不同的应力状态，根据应力状态选择恰当的裂纹模型，然后对不同的裂纹模型建立失效概率计算方程，并计及裂纹的尺寸、方向和密度等因素的影响。Thijssen等人[117]采用软件SILS（Sea Ice Loads Software），模拟结构物与亚北极区的浮冰、冰脊间的相互作用，进行冰载荷估算和敏感性分析。该软件采用的是ISO19906建议的概率方法，并结合蒙特卡罗算法。

5 其他方法

上述数值方法还可以互相结合以适用于更加广泛的问题。Paavilainen等人[118-121]采用有限元——离散元结合的方法（FEM-DEM）模拟较宽倾斜结构物周围发生的冰块堆积过程。邵帅等人[122-123]采用FEM-DEM方法对海洋结构物的冰激振动进行了研究。Timoshenko梁单元两端连接离散单元，梁间碰撞接触力也采用离散元的方法。数值方法分析船冰问题一般工作量较大且比较耗时，Ehlers等人[124]曾采用数值和解析相结合的方法降低计算成本。此外实地监测也是研究的重要手段，可以获得最直接的数据。

薛彦卓等人[125-127]采用近场动力学（peridynamics）对冰的三点弯曲试验、冰与加筋板的接触以及船舶在冰区的航行进行模拟。叶礼裕等人[128]采用近场动力学对冰——桨相互作用进行了模拟。Robb等人[129]结合平滑粒子流体动力学模型（SPH）与离散元方法模拟了含有漂浮固体物质的自由表面液体流动。

6 结 语

试验法是冰研究的重要方法，也是数值模拟的基础。有限元方法具有深厚的理论基础，应用广泛，相应的商业软件较多，方便实现；离散元方法的优点在于更适合模拟冰的破碎和离散特性；概率法则从统计学角度解释冰的随机性。船冰相互作用主要有两个问题：一是本构模型的开发与完善，二是数值方法的改进。对本构模型的开发一般需要试验数据的支持，也可参考相似材料和领域中的方法，例如混凝土和某些岩石。对数值方法的改进，例如离散元法，已有很多应用但目前开发得并不完善，此外，与冰相关的某些子问题可以采用统计学的方法，并与有限元或离散元方法相结合。