面对日益严峻的环境污染、温室气体排放和气候变化等问题, 许多国家都在调整能源发展战略, 煤电投资热度持续降温, 加快可再生能源的开发利用已成为目前世界各国的普遍共识。新世纪以来，风能发电在世界范围内迅速发展，并在可再生能源发电领域中占据重要地位。风能发电可以按照其安装的位置分为陆上风电和海上风电。相对陆上风电, 海上风电在风资源裕度、风电场位置、风力机装机容量等方面具有优势, 在风力机制造安装复杂度、全寿命周期成本、运行维护难度、对电网调峰裕度的要求等方面存在劣势。综合考虑后海上风电的利大于弊，具有很大开发价值，海上风电是世界众多国家推动能源低碳转型的核心路径。

海上风力发电在高纬度寒区海域，例如欧洲波罗的海、亚洲渤海湾、美洲五大湖等区域有很大的发展潜力，对于寒区海上风机而言海冰载荷是重要载荷之一。风能探测的报告显示我国渤海及黄海北部海域风能密度大，然而冬季会有海冰影响风电设备。大量工程项目表明, 海冰会对海洋结构产生巨大的危害。1969年渤海发生特大冰封, 巨大的海冰直接将海上的“海二井”石油平台推倒, 造成了建国以来最大的由于海冰导致的石油平台事故。放眼全球, 阿拉斯加库克湾的采油平台^[1]、北欧的Bothnia湾灯塔^[2]等冰区海洋结构都不同程度地遭到过海冰的破坏。因此可以说, 海冰是关乎海洋结构安全的全球性问题。

1 我国冰载荷的研究现状

目前我国关于冰载荷的研究发展尚不成熟，而且大多集中在海洋平台领域，海冰是高纬度海域海洋结构物设计和安全运营的重要影响因素，研究表明海冰载荷往往是寒区海域影响海洋工程结构物稳定的控制载荷^[3]。早在1966年，我国就开始在渤海海域建设抗冰的海上平台，至今已成功建成了数十座具有抗冰能力的海上平台。

冰载荷的研究只在少数的几个重视寒区海洋能源开发的国家开展，包括美国、北欧、加拿大、俄罗斯、日本、中国等，相对于波浪载荷、风载荷、地震载荷等问题的研究，冰载荷问题的工程需要相比之下没那么广泛，而且各个国家海域的冰情、采用的结构形式等均有着较大的差异，研究重点也根据实际海域有所不同，造成彼此研究的相对独立。我国的冰载荷研究主要围绕渤海展开，目前的研究主要分布在以下几个方面：

1) 渤海的冰情调查，包括历史冰情记录^[4-5]，冰情预报方法^[6]统计分析划分海冰区域以及给定工程设计参数^[7];

2) 海冰的物理力学性能^[8](包括温度、盐度、密度等物理参数和压缩强度、弯曲强度等力学参数);

3) 冰力及其检测方法^[9];

4) 冰与结构的动力相互作用和结构的冰振响应^[10];

5) 结构的冰激疲劳分析以及抗冰设计^[11]。

2 冰载荷数值仿真方法 2.1 数值仿真概况

针对冰载荷研究主要方式有现场实测、原型试验、理论分析和数值模拟。对冰载荷的研究主要针对静冰载荷和动冰载荷。由于海洋环境复杂，影响海冰运动的因素较多，因此如何建立合理有效的数值方法来模拟海冰动力学过程一直都是研究的重点，虽然海冰与结构物相互作用的过程至今没有公认的数学物理模型，不过研究者们还是基于模型实验和现场观测为基础，用数值模拟方法展现冰和结构相互作用的过程，特别是对海冰动态破坏的过程。

在分析冰与一般的结构物作用时，通常所研究的冰尺度较小。对于较小的尺度，若假设其为均匀介质就可以很简单地运用目前最为常用的数值分析方法有限元法^[17-21]。除此以外还有离散元法^[22-23]、非连续变形分析方法^[15-16]等其他数值分析方法在处理不同问题时有着各自优势。

2.2 有限元法

有限元分析是解决工程和数学模型时最为广泛使用的方法。将大型的系统细分为小而简单的元，在满足工程要求的前提下设置元素，将求解域划分为离散域，再根据离散域各自的方程进行求解，就可以用有限数量的未知量代数方程组去求解偏微分方程组。对于不同数学模型的问题有限元的求解方法是大致相同的，其核心思想均为“数值近似”和“离散化”，求解问题的基本步骤^[12]通常是根据问题确定好计算域，然后将计算域离散化，再确定状态变量的方程组和设定边界条件开始进行计算，计算得出结果后再进行分析。处理步骤大概可分成3个阶段，前处理、求解和后处理流程如图 1所示。前处理就是根据研究的问题建立有限元模型并且根据计算需求完成单元网格划分；后处理则是在得到计算结果以后根据计算的设定条件分析结果，以方便其他研究者了解。有限元方法是基于连续介质理论用来描述连续结构体运动应用最广泛的数值计算方法，适合用来分析海洋结构在冰力作用下动力响应。

2.3 离散元法

离散元法是专门用来解决不连续介质问题的数值模拟方法最初往往应用于岩石力学领域，把介质离散为元，相邻的元之间存在某种或几种作用力，元的运动满足牛顿第二定律。近年来，基于离散元方法模拟海冰与船体和海洋平台结构间的相互作用研究取得了较大进展。冰离散元方法的计算单元可采用球体、圆盘、多面体和扩展多面体等不同形态。采用离散单元模型计算更加适合模拟离散颗粒组合体在准静态或是动态条件下的变形和破碎过程，而海冰是属于非连续介质可以离散为具有一定质量的颗粒单元，颗粒之间允许位移重叠或分离，通过这些单元可以灵活构建不同类型的海冰，包括平整冰、浮冰、冰脊和堆积冰等。

2.4 其他数值方法

除了前述2种主流数值分析方法之外，还有许多数值方法用于冰载荷的研究，如差异元素法^[13]、三维颗粒流法^[13]、几何网格法^[14]等。石根华^[15]于1988年提出了一种新型数值分析方法将数学拓扑学和工程实践相结合的非连续变形分析(discontinuous deformation analysis，DDA)。该方法最初是为了解决岩体大变形和大位移问题时介质不连续性导致计算误差大的问题而提出的，平行于有限元法，充分考虑结构介质的不连续性，将结构面切割而成的块体作为分析单元，将动力学与静力学统一起来，用最小势能原理把块体之间的接触问题和块体本身的变形问题统一到矩阵中求解，后来被学者运用到冰与结构的相互作用模拟中。张运良等^[16]就基于DDA算法模拟冰与结构相互作用的过程，将冰破碎的过程可视化，计算了冰撞击结构后每个时间步的结构动态响应。模拟出的冰破碎过程与原型观测到的过程较为吻合，与有限元和边界元方法相比，DDA方法对冰的破碎模拟更为接近实际。不过目前整体而言该方法在应用于冰载荷研究还相对较少。

3 数值仿真方法在抗冰研究中的应用

在研究冰与海洋相互作用时所关注的问题往往是作用在结构上的最大静冰力和动冰力形式、海冰引起结构振动的机理以及结构在冰载荷作用下的失效形式，按照结构与冰作用部分的形状可以分为直立结构基础及非直立结构基础。现在数值仿真在冰载荷方面的应用往往是通过与模型试验数据和现场测试数据对比以验证数值仿真方法的适用性，积累冰载荷领域数值仿真的经验，并在此基础上进一步模拟仿真不同坏境、不同结构时冰与结构相互作用过程，为未来的抗冰结构设计奠定基础。

3.1 在直立结构中的应用

冰与直立结构作用时，可能出现劈裂、弯曲、屈曲和挤压等多种破坏形式，其中以挤压破坏为主。岳前进等^[17-18]、Kärnä等^[19]、李洪升等^[20]对冰与直立结构的作用过程及结构的振动进行了分析，发现冰与直立结构作用主要发生屈曲破坏和挤压破坏，并且不同速度的海冰与结构作用挤压破坏时会形成如图 2示的3种不同形式的冰激振动：准静态振动、稳态振动、随机振动，海冰的挤压破碎也因为冰速的不同而表现为3种不同的破碎模式。当冰速很慢时，出现准静态冰力，对应冰的破碎模式为准静态韧性破坏；随着冰速的增加会出现导致结构稳态振动的冰力，此时的破坏模式为韧脆转换破坏；随着冰速的继续增加会产生导致结构随机振动的冰力，对应破碎模式为脆性破坏。

在静冰力领域，目前已有相应的工程规范，杨耀鹏等^[21]以新建的海上平台为例，通过有限元软件分析其抗冰性能，发现结构设计静力的安全储备比较大，不过冰振对结构造成的潜在影响还有待进一步考究，可能会造成结构的疲劳失效。

在动冰力的研究领域，季顺迎等^[22]通过离散元法模拟计算了海冰与直立海洋平台相互作用的过程, 模拟海冰对不同桩径的直立结构的动冰力过程，确定海冰的破碎规律以及冰载荷的特性，发现平均冰力与桩径呈正比。贾宾等^[23]在前者研究的基础上，采用基于近场动力学理论的数值算法, 数值模拟的主要参数参考了前者的模型，模拟了海冰与柱状结构的作用过程，并考虑了相互作用时产生的结构振动，探究了桩径、冰速等参数对极值静冰力和结构振动位移的影响，冰力部分模拟结果与前者相近，可见这种新型的数值算法也是具有合理性的。

冰激振动是影响海洋结构安全以及导致其损伤和疲劳寿命的重要因素，也是冰载荷研究领域的热点。张大勇等^[24]在分析冰激振动下直立结构的疲劳寿命时，借助有限元数值分析方法，主要考虑了稳态振动和随机振动造成的疲劳损伤影响，选取安全寿命设计方法^[25]进行冰激疲劳估计，明确冰与直立结构相互作用过程，提出一套冰激直立抗冰平台的疲劳寿命计算流程。

在分析冰载荷问题时往往把海洋平台视为刚体结构，而在分析海洋结构的动力特性时则把冰载荷简化成冰力函数处理，随着冰激振动研究的深入，需要进一步研究海冰与海洋结构的动力耦合作用，不少研究者提出了建立离散元-有限元的耦合模型来分析冰激结构振动的内在机理。王帅霖^[26]采用基于区域分解方法的离散元-有限元耦合模型，对单桩直立结构的冰激振动进行了仿真模拟，由梁-壳组合单元构成海洋结构的有限元模型，在耦合界面上发展了相应的接触算法和耦合参数传递算法，考虑了冰速、厚度变化的影响，在冰况相近的前提下，仿真计算出冰载荷、冰激振动加速度和应力分布与现场实测数据趋势相同，初步说明了该耦合模型在分析冰与结构相互作用过程的适用性。

3.2 在非直立结构中的应用

近年来，锥体结构在海洋工程领域不断被采用，不少国家为冰区的海洋平台加装了图 3所示的锥体结构，国内外的研究人员针对冰排与固定倾斜结构作用时受到的冰载荷做了大量的理论、实验和数值研究，通过海冰的物理力学特性发现，冰的弯曲强度远小于挤压强度。锥体取代直立结构被公认为降低冰力的解决办法。

对于锥体结构而言，由于冰排的弯曲强度远小于挤压强度，而在根据现在设计标准，即使是直立结构也有很大的静力安全储备，而锥体则不言而喻，如图 4所示，对于锥体而言，虽然冰排破坏的力大幅减小了，不过会产生冰与斜面的摩擦力所带来的额外载荷，因此研究椎体结构静冰力时不仅需要考虑锥体本身的物理参数还需要考虑冰爬坡力所带来的额外载荷。

Sand等^[27]和Derradji-aouat^[28]，利用有限元软件ANSYS计算了冰排与斜坡结构作用时结构受到的冰力，将冰排理想化成作用在弹塑性基础上的线弹性板，将得到的结果与前人基于实验数据的理论公式进行比对。龚榆峰等^[29]在各项同性海冰失效准则的基础上结合Abaqus有限元软件，开发用于模拟海冰失效的新型单元并且结合了Ralston^[30]给出的塑性极限分析解进行比较, 直观地给出冰损伤过程及对应结构受力的变化过程。Bjørnar等^[27]使用非线性有限元方法模拟冰对固定圆锥体的作用，并分析锥角及摩擦力对冰荷载的影响。利用一种特殊的接触算法模拟了冰原与锥状结构之间的相互作用，使跟踪冰原与结构之间逐渐发展的接触成为可能。在相互作用过程中，由于冰原形态的变化，浮力也随之变化。通过引入一个连续的非线性地基模型来跟踪这一过程，该模型包括了浮力和冰的比重的影响。用2种不同的本构模型来近似冰的力学行为。将材料非线性和冰与结构摩擦的影响也纳入有限元计算当中，精确追踪冰盖与斜坡结构之间的接触，包括摩擦，计算部分或全部被淹没的冰盖上的浮力。探究了冰与结构之间锥角和摩擦系数的变化对结构的影响，并将数值计算的结果与基于塑性极限分析的分析结果进行了比较，发现两者非常接近，可见该有限元模型的可靠性。

锥体结构对减小动冰力以及减小冰激振动方面相比于直立结构表现出色，前文已经提及，冰激振动是目前寒区海洋工程抗冰结构设计的面临的重大课题，不过目前仍没有给出计算冰激结构振动的冰力函数，对海冰引起结构振动的力学机制的解释还没统一，研究者们也在不断探究适合计算海冰动力学的数值方法。王刚等^[31]利用有限元软件，对不同接触宽度的冰与锥体相互作用过程进行了分析，对径向开裂型和环向开裂型2种冰排弯曲断裂形式进行数值仿真计算，研究了接触宽度与冰排弯曲断裂形式的关系。李海等^[32]综合考虑了计算的效率和精确性，发展了一种将拉格朗日坐标与欧拉坐标相耦合的方法，将有限差分法和质点流体力学方法相结合，考虑了海冰的堆积以及涡旋风场的作用，对辽东湾海冰动力过程进行了72 h的数值模拟计算，并且与卫星的遥感图像进行比较，发现模拟的海冰的厚度、密集度分布以及冰速与卫星遥感图像有良好的一致性，可以较好地模拟出渤海海冰的动力学过程。

3.3 数值仿真结果的验证

目前对锥体海洋结构与海冰的相互作用主要通过现场实测、室内试验以及数值模拟展开研究，由于冰-锥相互作用过程的复杂性，很难通过现场实测和试验分析不同参数对该过程的影响，而数值模拟方法可以准确修改参数，再现出实测情况的冰-锥相互作用过程，并在确定数值仿真方法可行性的基础上进一步改变参数模拟各种情况，从而达到冰激振动预警分析、疲劳分析以及动力优化设计的目的。以上数值仿真工作大多基于已有的现场实测或室内试验数据仿真研究以确定数值仿真工作的合理性。关湃等^[33]基于模型试验结果，借助有限元软件对冰锥结构进行优化设计，模拟不同冰况下结构的响应，不断改变结构的各种形状参数以达到结构优化的目的。也有不少研究者专门设计试验以验证数值分析的结果是否合理，汪震宇^[34]用有限元方法分析斜坡与冰的相互过程，并设计了相应的模型试验进行误差分析以及模型的校正。

研究冰载荷的各种研究方法都是相互验证，不断改进的过程，而随着数值分析方法不断发展成熟，而且能够精确地改变各种参数以适应需求，对理解冰与结构相互作用机制具有成本低效率高的优势。

4 在海上风机一体化数值分析的应用 4.1 海上风机一体化分析需求

海上风机也作为一种海洋结构，在冰载荷方面的研究与前述的直立结构和锥面结构有诸多相似之处，在这便不再赘述了。海上风机相比于其他海洋结构，最大的特点在于风机为典型的柔性结构，海冰会引起柔性结构明显的动力放大，冰振响应十分显著，对结构及上部设备造成一定影响^[35]。张大勇等^[36]基于有限元分析法研究动冰载荷对结构的影响，选取单立柱式三桩形式的风电基础并与相同结构形式的导管架平台进行动力特性对比，文章中归纳总结了冰载荷作用下风电基础的各种失效模式，包括静冰力作用下的结构失效和动冰力作用下的疲劳失效、振动失效和土壤失效，发现现有的风电基础设计虽能满足极值冰载荷的要求，但在动冰力作用下对结构造成的影响还有待考量。通过有限元法仿真模拟可以得到验证风电基础结构主要是塔架在载荷下的响应，但并没有模拟出动冰力振动下对内部设备如电机、限速安全机构等的影响，但可以通过结构的冰振位移、速度等指标来间接地考虑动冰力下结构振动对内部设备的影响程度，为寒区风电基础的抗冰设计的改良、相关的实验设计及安全保障提供合理依据。

4.2 常用海上风机一体化分析软件及应用

目前用于海上风机一体化的软件主要为FAST软件，作为风机耦合的软件，不仅能够建立基础结构的模型还能搭建旋转叶片-机舱-塔架组件的详细模型，更加贴近实际的模拟仿真风机在各种载荷作用下的响应。FAST软件^[37]是美国可再生能源实验室(National Renewable Energy Laboratory，NREL)开发，用于计算气动弹性及结构动力学仿真的开源软件。该软件的仿真流程图如图 5所示，FAST作为一款开源的软件，为拓展研究、增加功能提供了基础，部分研究者在软件基础上进行二次开发以适应不同的研究需求。

现代风力机作为最大的旋转机械，扇叶长度甚至可以超过最大的飞机翼展，风载荷对风机的影响是其他海洋结构无法比拟的，所以分析风和海冰的联合载荷与风机作用的过程十分重要。许子非等^[38]考虑了海冰与湍流风载荷的联合作用，利用FAST软件，以近海单桩风力机为研究对象，数值计算在不同厚度及海冰移速的海冰碰撞条件下塔架的动力学响应特性，并对比了有无抗冰锥时响应的不同，发现抗冰结构能够明显减小结构的疲劳载荷。

除了FAST以外，目前用上风机一体化数值分析主流的软件还有HAWC2^[39] (horizontal axis wind turbine simulation code 2nd generation)和Bladed，都是被一些风力机厂商所认证的仿真软件，作为风机仿真软件应用于很多研究项目和工程当中，能够模拟许多现有的风力涡轮机。

HAWC2能够综合考虑空气动力学载荷、水动力载荷和土力学载荷，在时域内使用线性和非线性的方法结合以计算结构动力响应，基于多体公式的非线性气动弹性模型，可以处理复杂的结构。Shi等^[40]利用了HAWC2代码，建立了半经验的全耦合冰载荷模型，研究冰载荷对安装了抗冰结构的海上风电机的影响，同时考虑了风载荷和冰载荷的基础下进行了耦合分析并且与非耦合分析的结果进行了比较，讨论了耦合分析的重要性，并且研究了风速、平均冰厚、漂移速度不同时，冰的物理性质和破冰锥几何形状的对海上风机响应的影响。

除了冰与结构的相互作用之外，海上风力机仍需考虑结冰对风力机的影响，结冰会引起额外的载荷和振动，不同叶片上冰载荷的不平衡也会在一定程度上增加机组部件的疲劳。易贤等^[41]采用数值计算的方法进行结冰计算，模拟结冰形成的动态过程，分析结冰的主要部位以及厚度，还有对风机气动特性的影响，为设计结冰探测方法和除冰方案奠定基础。

风力机作为高耸结构，相比于其他海洋结构柔性特性更为明显，也意味着风力机受冰激振动的影响更为显著，而目前对于冰激振动的研究仍处于探索阶段，离彻底揭示冰激振动机理和抗冰振结构设计的完善还有很长的路要走，随着数值仿真技术的成熟，数值计算方法相比于原型观测和室内试验研究的进程更短，花费小效率高，能规避很多不必要的风险。

5 结束语

目前冰载荷研究的静冰力领域发展相对成熟，各个研究地区均有明确的规范，不过仍未形成统一计算公式或是数学模型，静冰力分析理论还有待进一步发展。新型的海洋结构或是抗冰结构对于静冰力均有较大的安全储备，应当探索降低静冰力工程设计值以提高经济效益的可能性；在动冰力和冰激振动领域，仍需探索可以用于工程的动冰力函数，进一步优化抗冰结构的减震性能，形成相关的工程规范，并且发展数值分析方法以适应动冰力和冰激振动数值计算需求。

理论分析和数值模型的发展必然离不开工程实例与试验，目前我国仅天津大学有冰池实验室以及哈尔滨工程大学建设的室外冰水池实验室，但与德国汉堡水池相比仍有着显著的差距，有待进一步改善实验条件，开展更多冰力实验。

当然，冰载荷领域的研究要想应用到海上风机还需要结合其他领域知识，诸如空气动力学、水动力学、土力学等领域的研究，并结合现代的数值仿真技术，不断发展一体化数值分析技术以贴近实际。

相信随着数值模拟技术与计算机的蓬勃发展，以及目前各国对冰载荷愈加重视，冰载荷的数值仿真方法也会越来越成熟，同时也会有更多的文献方便学习。我国渤海海域能源丰富，不论是传统能源还是新能源的开采都避不开桩式结构，而渤海易受海冰影响，因此我国在抗冰方面的步伐从未停止，特别是近年来各地兴起了海上风电的热潮。通将过数值方法与模型试验、现场观测等手段结合不断发展冰力研究，并在工程实践中不断积累相关经验，相信在不久之后便会形成一套完善海上风机的冰力规范。